拓端tecdat|R語言非線性混合效應 NLME模型(固定效應&隨機效應)對抗哮喘葯物茶鹼動力學研究

原文鏈接：http://tecdat.cn/?p=24074

原文出處：拓端數據部落公眾號

簡介

茶鹼數據

茶鹼數據文件報告來自抗哮喘葯物茶鹼動力學研究的數據。給 12 名受試者口服茶鹼，然后在接下來的 25 小時內在 11 個時間點測量血清濃度。

head(thdat)

此處，時間是從抽取樣品時開始給葯的時間（h），濃度是測得的茶鹼濃度（mg/L），體重是受試者的體重（kg）。

12 名受試者在時間 0 時接受了 320 mg 茶鹼。

讓我們繪制數據，即濃度與時間的關系：

plot(data=theo.data2) +eo_ine(oaes(group=id))

數據的個體差異

我們還可以在 12 個單獨的圖上繪制 12 個單獨的濃度分布圖，

pl + geom_line() + facet_wrap(~id)

這12個人的模式是相似的：濃度首先在吸收階段增加，然后在消除階段減少。然而，我們清楚地看到這些曲線之間的一些差異，這不僅僅是由於殘差造成的。我們看到病人吸收和消除葯物的速度或多或少。

一方面，每個單獨的特征將通過非線性 葯代動力學 (PK) 模型正確描述 。

另一方面，人口方法和混合效應模型的使用將使我們能夠考慮這種 個體間的變異性。

將非線性模型擬合到數據

將非線性模型擬合到單個受試者

讓我們考慮本研究的第一個受試者（id=1）

1.  
   the.dat.dta$id==1 ,c("tme)]
2.  
   plot(data=teo1

 我們可能想為這個數據擬合一個 PK 模型

其中 (yj,1≤j≤n) 是該受試者的 nn PK 測量值，f 是 PK 模型，ψ是該受試者的 PK 參數向量， (ej,1≤ j≤n)是殘差。

對該數據寫入具有一階吸收和線性消除的單室模型

 

其中 ψ=(ka,V,ke) 是模型的 PK 參數，D 是給予患者的葯物量（此處，D=320mg）。

讓我們計算定義為 ψ 的最小二乘估計

我們首先需要實現PK模型：

1.  
   pk.od <- function(pi, t){
2.  
   D <- 320
3.  
   ka
4.  
   V
5.  
   ke
6.  
   f <- D*a/V/(a-k)*(exp(-e*t)-exp(-k*t))

然后我們可以使用該 nls 函數將此（非線性）模型擬合到數據

1.  
   nls(neatin ~p.me1(psi, time))
2.  
   coef(km1)

並繪制預測濃度 f(t,ψ^)

1.  
   e. <- dafme(tm=sq(0,40,=.2))
2.  
   w.pd1 <- pedct(pk, newaa=wdf)
3.  
   line(da=new., aes(x=tie,y=re1))

將獨特的非線性模型擬合到幾個患者上

與其將這個 PK 模型擬合到單個患者，我們可能希望將相同的模型擬合到所有患者：

其中（yij,1≤j≤ni）是受試者i的ni PK測量值。這里，ψ是N個受試者共享的PK參數的向量。

在該模型中，ψ 的最小二乘估計定義為

讓我們將該nls 函數與來自 12 個受試者的合並數據一起使用 。

 nls(ocetn ~ kme1(ps, tme)

1.  
   nll <- predct(kmll, ewta=n.f)
2.  
   p+geom_line(ewd,astm,=rdal,clu="390" )

這些估計的 PK 參數是典型的 PK 參數，並且該 PK 曲線是該患者樣本的典型 PK 曲線。

根據定義，它們沒有考慮患者之間的變異性，因此不能提供良好的個體預測。

line(data=e.d, aes(x=im,y=pe.al)) + faetap(~ id)

將多個非線性模型擬合到多個受試者

相反，我們可以為每個受試者擬合具有不同參數的相同 PK 模型，正是我們在上面對第一個患者所做的：

其中 ψi 是患者 ii 的 PK 參數向量。

在該模型中，ψi 的最小二乘估計定義為

1.  
    
2.  
   for (i in (1:N)) {
3.  
   pkmi <- nls(cocetatn ~ pk.mdl1(psi, time)
4.  
   pred <- c(prd, prdit(kmi, neta=ewf))
5.  
   }

每個個體預測濃度 f(t,ψ^i)似乎很好地預測了 12 個受試者的觀察濃度：

1.  
   nc <- lengh(nwdtie)
2.  
   tepred <- data.rame(d=rp(1:12),acc),tie=renew.fime12 fpre=pre)
3.  
   line(dta=te.re, aes(x=me,y=frd)) + factrp(id)

非線性混合效應 (NLME) 模型

第一個基本模型

到目前為止，單個參數 (ψi)被認為是固定效應：我們沒有對可能的值做出任何假設。

在群體方法中，假設 N 受試者是從相同的個體群體中隨機抽樣的。然后，每個單獨的參數 ψi 被視為一個隨機變量。

我們將開始假設 ψi是獨立且正態分布的：

其中 ψpop 是總體參數的 d 向量，Ω是  d×d方差-協方差矩陣。

備注： 這個正態性假設允許我們將每個單獨的參數 ψi 分解為固定效應 ψpop 和隨機效應 ηi：

其中 ηi∼iidN(0,Ω)。

我們還將開始假設殘差 (eij)是獨立且正態分布的：eij∼iidN(0,a2)。

總之，我們可以等效地表示一個（非線性）混合效應模型

i) 使用方程：

其中 eij∼iidN(0,a2) 和 ηi∼iidN(0,Ω),

ii) 或使用概率分布：

模型是(y,ψ)的聯合概率分布，其中y=(yij,1≤i≤N,1≤j≤ni)是完整的觀測集，ψ=(ψi,1≤i≤N) 單個參數的 N向量，

任務、方法和算法

總體參數的估計

模型參數為θ=(ψpop,Ω,a2)。θ的最大似然估計包括使似然函數相對於 θ 最大化，  定義為

如果f是ψi的非線性函數，那么yi就不是高斯向量，似然函數L(θ,y)就不能以封閉形式計算。

在非線性混合效應模型中存在幾種最大似然估計的算法。特別是，隨機近似EM算法（SAEM）是一種迭代算法，在一般條件下收斂到似然函數的最大值。
 

單個參數的估計

一旦θ被估計出來，條件分布p(ψi|yi;θ^)就可以用於每個個體i來估計個體參數向量ψi。

這個條件分布的模式被定義為

該估計稱為 ψi 的最大后驗 (MAP) 估計或經驗貝葉斯估計 (EBE)。

備注： 由於 f 是 ψi的非線性函數，因此沒有 ψ^i的解析表達式。然后應使用牛頓算法來執行此最小化問題。

然后我們可以使用條件模式來計算預測，采取的理念是各個參數的最可能值最適合計算最可能的預測。

似然函數的估計

對給定模型執行似然比檢驗和計算信息標准需要計算對數似然 

對於非線性混合效應模型，不能以封閉形式計算對數似然。在連續數據的情況下，通過高斯線性模型近似模型允許我們近似對數似然。

實際上，我們可以將個體 i的觀測值 (yij,1≤j≤ni)的模型線性化，該模型圍繞預測的個體參數 ψ^i 的向量。

設∂ψf(t,ψ)是f(t,ψ)關於ψ的導數的行向量。然后，

在此之后，我們可以通過正態分布來近似向量 yi 的邊緣分布：

其中

然后對數似然函數近似為

Fisher信息矩陣的估計

使用線性化模型，最大似然估計 (MLE) θ^ 的方差以及置信區間可以從觀察到的 Fisher 信息矩陣 (FIM) 中導出，而 FIM 本身是從觀察到的似然導出的：

然后可以通過觀察到的 FIM 的逆來估計 θ^ 的方差-協方差矩陣。θ^ 的每個分量的標准誤差 (se) 是標准偏差，即方差-協方差矩陣的對角元素的平方根。

對茶鹼數據擬合 NLME 模型

讓我們看看如何將我們的模型擬合到茶鹼數據。

我們首先需要定義應該使用數據文件的哪一列以及它們的作用。在我們的示例中，濃度是因變量 yy，時間是解釋變量（或預測變量）t，id 是分組變量。

1.  
   Data(dta = data,
2.  
   grp = id",
3.  
   prditors = "time",
4.  
   repose = "con")

結構模型是以前使用的一階吸收和線性消除的單室模型。

1.  
   molct <- function(pi,id,x) {
2.  
   D <- 320
3.  
   fe <-D*a/(V*(a-e))*(exp(-e*t)-exp(-a*t))

需要人口參數向量ψpop的結構模型和一些初始值

1.  
   Model(modl = moelpt,
2.  
   pi = c(a=1,V=20,ke=0.5))

可以定義幾個選擇和運行算法的選項，包括單個參數的估計 (map=TRUE)、Fisher 信息矩陣的估計和線性化對數似然 (fim=TRUE) 或重要性采樣的對數似然（ll.is=TRUE）。

種子是用於隨機數生成器的整數：使用相同的種子多次運行算法可確保結果相同。

1.  
   list(map=TRUE,seed=632545)
2.  
   mix(model, dat,optns)

可以顯示估計算法的結果摘要

results

還可以使用單個參數估計值

這些單獨的參數估計可用於計算和繪制單獨的預測

1.  
   pred(fit1)
2.  
   plot.fit(fit1)

可以顯示多個診斷擬合圖，包括觀察值與單個預測的圖

pltobsv(fit1,lvl=1)

殘差與時間和個人預測的關系圖，

pltsateresi(fit1, levl=1)

模型的一些擴展

殘差模型

在模型 yij=f(tij,ψi)+eij 中，假設殘差 (eij)是均值為 0 的高斯隨機變量。 (eij)在非線性混合效應模型中的方差。

恆定誤差模型：

殘差 (eij) 是獨立同分布的：

因此， yij 的方差隨時間保持不變：

其中 εij∼iidN(0,1)。

誤差模型可以定義為Model 的參數

Model(mo=md1p, p0=c(ka=1,V=20,ke=0.5), mdl="constant")

比例誤差模型：

比例誤差模型假設 eij的標准偏差與預測因變量成正比： eij= bf(tij,ψi)εij 其中 εij∼iidN(0,1)。然后，

Model(modl=dl1pt,error="prori")

組合誤差模型：

組合誤差模型將常數和比例誤差模型相加組合： eij=(a+ bf(tij,ψi))εij其中 εij∼iidN(0,1)。然后，

Model(moel=d1ct, mde="bined")

指數誤差模型：

如果已知 y 取非負值，則可以使用對數轉換。然后我們可以用兩個等效表示來編寫模型：

Model( ero.dl="exp")

單個參數的變換

顯然，並非所有分布都是高斯分布。首先，正態分布有支持度R，與許多在精確區間取值的參數不同。例如，有些變量只取正值（如體積和轉移率常數），其他變量則被限制在有界區間內。

此外，高斯分布是對稱的，這並不是所有分布都具有的屬性。擴展使用高斯分布的一種方法是考慮我們感興趣的參數的某種變換是高斯的。

即假設存在一個單調的函數h，使得h(ψi)是正態分布。為了簡單起見，我們在這里將考慮一個標量參數ψi。然后我們假設

或者，等效地，

其中 ηi∼N(0,ω2)。

對數正態分布：

對數正態分布確保非負值，廣泛用於描述生理參數的分布。

如果 ψi服從對數正態分布，則以下 3 種表示是等價的：

對數正態分布：

logit 函數定義在 (0,1)上並取其在 RR 中的值：對於 (0,1)中的任何 x，

具有 logit 正態分布的單個參數 ψi 在 (0,1)中取值。ψ 的 logit 服從正態分布，即，

概率正態分布：

probit函數是與標准正態分布N(0,1)相關的反累積分布函數（量化函數）ψ-1。對於(0,1)中的任何x。

具有概率正態分布的單個參數 ψi 在 (0,1) 中取值。ψi的概率呈正態分布：

每個單獨參數的分布可以使用參數 transform.par 定義（0=normal，1=log-normal，2=probit，3=logit）。默認為正態分布，即向量為 0。

例如，如果我們想使用 V 的正態分布和 ka 和 ke 的對數正態分布，那么 par 應該是向量 c(1,0,1)：

1.  
   Model(model ,
2.  
   psi ,
3.  
   trns.par = c(1,0,1))
4.  
    

備注：這里，ω2ka和ω2ke是log（kai）和log（kei）的方差，而ω2V是Vi的方差。

帶有協變量的模型

讓ci=(ci1,ci2,...,ciL)為個體協變量的向量，即數據中可獲得的個體參數的向量。我們可能想用這些協變量來解釋非觀察到的個體參數（ψi）的部分變異性。

我們將只考慮協變量的線性模型。更准確地說，假設 h(ψi) 是正態分布的，我們將 h(ψi)分解為固定效應和隨機效應：

備注:如果協變量ci1, ..., ciL對人口中的典型個體來說為零，ψpop就是ψi的典型值。

讓我們考慮一個模型，其中體積Vi是正態分布，是重量wi的線性函數。

假設人口中一個典型個體的體重是wpop，這個個體的預測體積不是β0，而是β0+βwpop。

如果我們使用中心體重wi-wpop，我們現在可以把模型寫成

事實上，現在對一個典型個體的預測體積是Vpop。

假設我們決定在茶鹼研究中使用70公斤作為典型體重。現在需要包括wi-70。

這里，只有體積 VV 是重量的函數。因此，協變量模型被編碼為向量 (0,1,0)。

1.  
   Model(
2.  
   trasf = c(1,0,1),
3.  
   covri = c(0,1,0))

這里，β^w70=0.33意味着重量增加1kg會導致預測的體積增加0.33l。

檢驗H0：βw70=0與H1：βw70≠0的P值為0.01，那么我們可以拒絕H0，並得出結論：預測的體積隨着重量的增加而顯著增加。

想象一下，我們現在用對數正態分布來表示體積Vi。現在是對數體積，它是轉化后的重量的一個線性函數。

我們可以假設，例如，對數體積是中心對數重量的線性函數。

或者，等效地，

我們看到，使用這個模型，一個典型個體的預測體積是Vpop。

Data對象現在需要包括log（wi/70）這個協變量。

1.  
   lw70 <- log(weight/70)
2.  
   Data(data,
3.  
   res=c("cerato"),
4.  
   cova=c("lw70"))

協變量模型再次編碼為（行）向量 (0,1,0)，但變換現在對於三個參數編碼為 1

1.  
   Model(
2.  
   trans.pr = c(1,1,1),
3.  
   cor = c(0,1,0))
4.  
    

隨機效應之間的相關性

到目前為止，隨機效應被認為是不相關的，即矢量-協方差矩陣Ω是一個對角矩陣。

隨機效應之間的相關性可以通過輸入參數covari引入，這是一個大小等於模型中參數數量的方形矩陣，給出了模型的方差-協方差結構。1s對應於估計的方差（在對角線上）或協方差（非對角線元素）。矩陣Ω的結構應該是塊狀的。

例如，考慮一個模型，其中ka在人群中是固定的，即ωka=0（因此對所有i來說kai=0），而log(V)和log(ke)是相關的，即ηV和ηke)是相關的。

1.  
   Model(
2.  
   covai = t(c(0,1,0)),
3.  
   covain = matrix(c(0,0,0,0,1,1,0,1,1),nrow=3))
4.  
    

---

最受歡迎的見解

1.基於R語言的lmer混合線性回歸模型

2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

3.R語言線性混合效應模型實戰案例

4.R語言線性混合效應模型實戰案例2

5.R語言線性混合效應模型實戰案例

6.線性混合效應模型Linear Mixed-Effects Models的部分折疊Gibbs采樣

7.R語言LME4混合效應模型研究教師的受歡迎程度

8.R語言中基於混合數據抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM