全文下載鏈接：http://tecdat.cn/?p=24074

最近我們被客戶要求撰寫關(guān)于非線性混合效應(yīng) NLME模型的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

茶堿數(shù)據(jù)文件報(bào)告來自抗哮喘藥物茶堿動(dòng)力學(xué)研究的數(shù)據(jù)。給 12 名受試者口服茶堿，然后在接下來的 25 小時(shí)內(nèi)在 11 個(gè)時(shí)間點(diǎn)測(cè)量血清濃度

head(thdat)

此處，時(shí)間是從抽取樣品時(shí)開始給藥的時(shí)間（h），濃度是測(cè)得的茶堿濃度（mg/L），體重是受試者的體重（kg）。

12 名受試者在時(shí)間 0 時(shí)接受了 320 mg 茶堿。

讓我們繪制數(shù)據(jù)，即濃度與時(shí)間的關(guān)系：

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

相關(guān)視頻

**

拓端

，贊14

**

拓端

，贊21

**

拓端

，贊7

**

拓端

，贊17

數(shù)據(jù)的個(gè)體差異

我們還可以在 12 個(gè)單獨(dú)的圖上繪制 12 個(gè)單獨(dú)的濃度分布圖，

pl?+?geom_line()?+?facet_wrap(~id)

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

一方面，每個(gè)單獨(dú)的特征將通過_非線性_?藥代動(dòng)力學(xué) (PK) 模型正確描述?。

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

將非線性模型擬合到數(shù)據(jù)

將非線性模型擬合到單個(gè)患者

讓我們考慮本研究的第一個(gè)主題（id=1）

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

?我們可能想為這個(gè)數(shù)據(jù)擬合一個(gè) PK 模型

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

對(duì)該數(shù)據(jù)寫入具有一階吸收和線性消除的單室模型

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

讓我們計(jì)算定義為 ψ?的最小二乘估計(jì)

我們首先需要實(shí)現(xiàn)PK模型：

pk.od?<-?function(pi,?t){??D??<-?320 ??ka? ??V?? ??ke? ??f??<-?D*a/V/(a-k)*(exp(-e*t)-exp(-k*t))

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

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

并繪制預(yù)測(cè)濃度 f(t,ψ^)

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

將獨(dú)特的非線性模型擬合到幾個(gè)患者上

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

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

在該模型中，ψ?的最小二乘估計(jì)定義為

讓我們將該nls?函數(shù)與來自 12 個(gè)受試者的合并數(shù)據(jù)一起使用?。

?nls(ocetn?~?kme1(ps,?tme)

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

這些估計(jì)的 PK 參數(shù)是典型的 PK 參數(shù)，并且該 PK 曲線是該患者樣本的典型 PK 曲線。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語言LME4混合效應(yīng)模型研究教師的受歡迎程度

左右滑動(dòng)查看更多

01

02

03

04

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

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

將多個(gè)非線性模型擬合到多個(gè)患者

相反，我們可以為每個(gè)受試者擬合具有不同參數(shù)的相同 PK 模型，正是我們?cè)谏厦鎸?duì)第一個(gè)患者所做的：

其中 ψi 是患者 ii 的 PK 參數(shù)向量。

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

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

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

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

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

第一個(gè)基本模型

到目前為止，單個(gè)參數(shù) (ψi)被認(rèn)為是固定效應(yīng)：我們沒有對(duì)可能的值做出任何假設(shè)。

在群體方法中，假設(shè) N 受試者是從相同的個(gè)體群體中隨機(jī)抽樣的。然后，每個(gè)單獨(dú)的參數(shù) ψi 被視為一個(gè)隨機(jī)變量。

我們將開始假設(shè) ψi是獨(dú)立且正態(tài)分布的：

其中 ψpop 是總體參數(shù)的 d 向量，Ω是? d×d方差-協(xié)方差矩陣。

備注：?這個(gè)正態(tài)性假設(shè)允許我們將每個(gè)單獨(dú)的參數(shù) ψi 分解為固定效應(yīng) ψpop 和隨機(jī)效應(yīng) ηi：

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

我們還將開始假設(shè)殘差 (eij)是獨(dú)立且正態(tài)分布的：eij～iidN(0,a2)。

總之，我們可以等效地表示一個(gè)（非線性）混合效應(yīng)模型

i)?使用方程：

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

ii)?或使用概率分布：

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

任務(wù)、方法和算法

總體參數(shù)的估計(jì)

模型參數(shù)為θ=(ψpop,Ω,a2)。θ的最大似然估計(jì)包括使_似然函數(shù)_相對(duì)于 θ 最大化，??定義為

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

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

單個(gè)參數(shù)的估計(jì)

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

這個(gè)條件分布的模式被定義為

該估計(jì)稱為 ψi 的最大后驗(yàn) (MAP) 估計(jì)或經(jīng)驗(yàn)貝葉斯估計(jì) (EBE)。

備注：?由于 f 是 ψi的非線性函數(shù)，因此沒有 ψ^i的解析表達(dá)式。然后應(yīng)使用牛頓算法來執(zhí)行此最小化問題。

然后我們可以使用條件模式來計(jì)算預(yù)測(cè)，采取的理念是各個(gè)參數(shù)的最可能值最適合計(jì)算最可能的預(yù)測(cè)。

似然函數(shù)的估計(jì)

對(duì)給定模型執(zhí)行似然比檢驗(yàn)和計(jì)算信息標(biāo)準(zhǔn)需要計(jì)算對(duì)數(shù)似然?

對(duì)于非線性混合效應(yīng)模型，不能以封閉形式計(jì)算對(duì)數(shù)似然。在連續(xù)數(shù)據(jù)的情況下，通過高斯線性模型近似模型允許我們近似對(duì)數(shù)似然。

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

設(shè)?ψf(t,ψ)是f(t,ψ)關(guān)于ψ的導(dǎo)數(shù)的行向量。然后，

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

其中

然后對(duì)數(shù)似然函數(shù)近似為

Fisher信息矩陣的估計(jì)

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

然后可以通過觀察到的 FIM 的逆來估計(jì) θ^?的方差-協(xié)方差矩陣。θ^ 的每個(gè)分量的標(biāo)準(zhǔn)誤差 (se) 是標(biāo)準(zhǔn)偏差，即方差-協(xié)方差矩陣的對(duì)角元素的平方根。

對(duì)茶堿數(shù)據(jù)擬合 NLME 模型

讓我們看看如何將我們的模型擬合到茶堿數(shù)據(jù)。

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

Data(dta???????=?data, ??????????????????????????grp??????=?id", ??????????????????????????prditors?=?"time", ??????????????????????????repose???=?"con")

結(jié)構(gòu)模型是以前使用的一階吸收和線性消除的單室模型。

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

需要人口參數(shù)向量ψpop的結(jié)構(gòu)模型和一些初始值

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

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

種子是用于隨機(jī)數(shù)生成器的整數(shù)：使用相同的種子多次運(yùn)行算法可確保結(jié)果相同。

list(map=TRUE,seed=632545)mix(model,?dat,optns)

可以顯示估計(jì)算法的結(jié)果摘要

results

還可以使用單個(gè)參數(shù)估計(jì)值

這些單獨(dú)的參數(shù)估計(jì)可用于計(jì)算和繪制單獨(dú)的預(yù)測(cè)

pred(fit1) plot.fit(fit1)

可以顯示多個(gè)診斷擬合圖，包括觀察值與單個(gè)預(yù)測(cè)的圖

pltobsv(fit1,lvl=1)

殘差與時(shí)間和個(gè)人預(yù)測(cè)的關(guān)系圖，

pltsateresi(fit1,?levl=1)

模型的一些擴(kuò)展

殘差模型

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

恒定誤差模型：

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

因此， yij 的方差隨時(shí)間保持不變：

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

誤差模型可以定義為Model 的參數(shù)

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

比例誤差模型：

比例誤差模型假設(shè) eij的標(biāo)準(zhǔn)偏差與預(yù)測(cè)因變量成正比：eij= bf(tij,ψi)εij 其中 εij～iidN(0,1)。然后，

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

組合誤差模型：

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

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

指數(shù)誤差模型：

如果已知 y 取非負(fù)值，則可以使用對(duì)數(shù)轉(zhuǎn)換。然后我們可以用兩個(gè)等效表示來編寫模型：

Model(?ero.dl="exp")

單個(gè)參數(shù)的變換

顯然，并非所有分布都是高斯分布。首先，正態(tài)分布有支持度R，與許多在精確區(qū)間取值的參數(shù)不同。例如，有些變量只取正值（如體積和轉(zhuǎn)移率常數(shù)），其他變量則被限制在有界區(qū)間內(nèi)。

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

即假設(shè)存在一個(gè)單調(diào)的函數(shù)h，使得h(ψi)是正態(tài)分布。為了簡單起見，我們?cè)谶@里將考慮一個(gè)標(biāo)量參數(shù)ψi。然后我們假設(shè)

或者，等效地，

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

對(duì)數(shù)正態(tài)分布：

對(duì)數(shù)正態(tài)分布確保非負(fù)值，廣泛用于描述生理參數(shù)的分布。

如果 ψi服從對(duì)數(shù)正態(tài)分布，則以下 3 種表示是等價(jià)的：

對(duì)數(shù)正態(tài)分布：

logit 函數(shù)定義在 (0,1)上并取其在 RR 中的值：對(duì)于 (0,1)中的任何 x，

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

概率正態(tài)分布：

probit函數(shù)是與標(biāo)準(zhǔn)正態(tài)分布N(0,1)相關(guān)的反累積分布函數(shù)（量化函數(shù)）ψ-1。對(duì)于(0,1)中的任何x。

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

每個(gè)單獨(dú)參數(shù)的分布可以使用參數(shù) transform.par 定義（0=normal，1=log-normal，2=probit，3=logit）。默認(rèn)為正態(tài)分布，即向量為 0。

例如，如果我們想使用 V 的正態(tài)分布和 ka 和 ke 的對(duì)數(shù)正態(tài)分布，那么 par 應(yīng)該是向量 c(1,0,1)：

Model(model??, ??????????????????????????psi???, ??????????????????????????trns.par?=?c(1,0,1))

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

帶有協(xié)變量的模型

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

我們將只考慮協(xié)變量的線性模型。更準(zhǔn)確地說，假設(shè) h(ψi) 是正態(tài)分布的，我們將 h(ψi)分解為固定效應(yīng)和隨機(jī)效應(yīng)：

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

讓我們考慮一個(gè)模型，其中體積Vi是正態(tài)分布，是重量wi的線性函數(shù)。

假設(shè)人口中一個(gè)典型個(gè)體的體重是wpop，這個(gè)個(gè)體的預(yù)測(cè)體積不是β0，而是β0+βwpop。

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

事實(shí)上，現(xiàn)在對(duì)一個(gè)典型個(gè)體的預(yù)測(cè)體積是Vpop。

假設(shè)我們決定在茶堿研究中使用70公斤作為典型體重。現(xiàn)在需要包括wi-70。

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

Model( ????????????????????????????trasf???=?c(1,0,1), ????????????????????????????covri?=?c(0,1,0))

這里，β^w70=0.33意味著重量增加1kg會(huì)導(dǎo)致預(yù)測(cè)的體積增加0.33l。

檢驗(yàn)H0：βw70=0與H1：βw70≠0的P值為0.01，那么我們可以拒絕H0，并得出結(jié)論：預(yù)測(cè)的體積隨著重量的增加而顯著增加。

想象一下，我們現(xiàn)在用對(duì)數(shù)正態(tài)分布來表示體積Vi。現(xiàn)在是對(duì)數(shù)體積，它是轉(zhuǎn)化后的重量的一個(gè)線性函數(shù)。

我們可以假設(shè)，例如，對(duì)數(shù)體積是中心對(duì)數(shù)重量的線性函數(shù)。

或者，等效地，

我們看到，使用這個(gè)模型，一個(gè)典型個(gè)體的預(yù)測(cè)體積是Vpop。

Data對(duì)象現(xiàn)在需要包括log（wi/70）這個(gè)協(xié)變量。

lw70?<-?log(weight/70) Data(data, ????????????????????????res=c("cerato"), ????????????????????????cova=c("lw70"))

協(xié)變量模型再次編碼為（行）向量 (0,1,0)，但變換現(xiàn)在對(duì)于三個(gè)參數(shù)編碼為 1

Model( ??????????????????????????trans.pr???=?c(1,1,1), ??????????????????????????cor?=?c(0,1,0))

隨機(jī)效應(yīng)之間的相關(guān)性

到目前為止，隨機(jī)效應(yīng)被認(rèn)為是不相關(guān)的，即矢量-協(xié)方差矩陣Ω是一個(gè)對(duì)角矩陣。

隨機(jī)效應(yīng)之間的相關(guān)性可以通過輸入?yún)?shù)covari引入，這是一個(gè)大小等于模型中參數(shù)數(shù)量的方形矩陣，給出了模型的方差-協(xié)方差結(jié)構(gòu)。1s對(duì)應(yīng)于估計(jì)的方差（在對(duì)角線上）或協(xié)方差（非對(duì)角線元素）。矩陣Ω的結(jié)構(gòu)應(yīng)該是塊狀的。

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

Model( ??????????????????????????covai?=?t(c(0,1,0)), ??????????????????????????covain?=?matrix(c(0,0,0,0,1,1,0,1,1),nrow=3))

本文摘選?《?R語言非線性混合效應(yīng) NLME模型(固定效應(yīng)&隨機(jī)效應(yīng))對(duì)抗哮喘藥物茶堿動(dòng)力學(xué)研究?》?，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語言用線性混合效應(yīng)（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關(guān)系
R語言LME4混合效應(yīng)模型研究教師的受歡迎程度R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實(shí)例
R語言混合線性模型、多層次模型、回歸模型分析學(xué)生平均成績GPA和可視化
R語言線性混合效應(yīng)模型（固定效應(yīng)&隨機(jī)效應(yīng)）和交互可視化3案例
R語言用lme4多層次（混合效應(yīng)）廣義線性模型（GLM），邏輯回歸分析教育留級(jí)調(diào)查數(shù)據(jù)R語言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例
R語言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應(yīng)模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言建立和可視化混合效應(yīng)模型mixed effect model
R語言LME4混合效應(yīng)模型研究教師的受歡迎程度
R語言 線性混合效應(yīng)模型實(shí)戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
SPSS中的多層（等級(jí)）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)
用SPSS估計(jì)HLM多層（層次）線性模型模型