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

原文出處：拓端數(shù)據(jù)部落公眾號

簡介

茶堿數(shù)據(jù)

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

head(thdat)

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

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

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

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

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

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

pl + geom_line() + facet_wrap(~id)

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

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

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

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

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

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

1. the.dat.dta$id==1 ,c("tme)]

2. plot(data=teo1

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

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

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

?

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

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

我們首先需要實(shí)現(xiàn)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?函數(shù)將此（非線性）模型擬合到數(shù)據(jù)

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

2. coef(km1)

并繪制預(yù)測濃度 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))

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

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

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

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

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

nls(ocetn ~ kme1(ps, tme)

1. nll <- predct(kmll, ewta=n.f)

2. p+geom_line(ewd,astm,=rdal,clu="390" )

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

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

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

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

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

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

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

1. 


2. for (i in (1:N)) {

3. pkmi <- nls(cocetatn ~ pk.mdl1(psi, time)

4. pred <- c(prd, prdit(kmi, neta=ewf))

5. }

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

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)

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

第一個(gè)基本模型

到目前為止，單個(gè)參數(shù) (ψi)被認(rèn)為是固定效應(yīng)：我們沒有對可能的值做出任何假設(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)是完整的觀測集，ψ=(ψi,1≤i≤N) 單個(gè)參數(shù)的 N向量，

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

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

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

如果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ù)測，采取的理念是各個(gè)參數(shù)的最可能值最適合計(jì)算最可能的預(yù)測。

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

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

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

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

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

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

其中

然后對數(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é)方差矩陣的對角元素的平方根。

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

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

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

1. Data(dta ? ? ? = data,

2. grp ? ? ?= id",

3. prditors = "time",

4. repose ? = "con")

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

1. molct <- function(pi,id,x) {

2. D ? <- 320

3. fe <-D*a/(V*(a-e))*(exp(-e*t)-exp(-a*t))

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

1. Model(modl = moelpt,

2. pi ?= c(a=1,V=20,ke=0.5))

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

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

1. list(map=TRUE,seed=632545)

2. mix(model, dat,optns)

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

results

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

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

1. pred(fit1)

2. plot.fit(fit1)

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

pltobsv(fit1,lvl=1)

殘差與時(shí)間和個(gè)人預(yù)測的關(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ù)測因變量成正比： 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ù)值，則可以使用對數(shù)轉(zhuǎn)換。然后我們可以用兩個(gè)等效表示來編寫模型：

Model( ero.dl="exp")

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

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

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

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

或者，等效地，

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

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

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

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

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

logit 函數(shù)定義在 (0,1)上并取其在 RR 中的值：對于 (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。對于(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 的對數(shù)正態(tài)分布，那么 par 應(yīng)該是向量 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的方差。

帶有協(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對人口中的典型個(gè)體來說為零，ψpop就是ψi的典型值。

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

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

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

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

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

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

1. Model(

2. trasf ? = c(1,0,1),

3. covri = c(0,1,0))

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

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

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

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

或者，等效地，

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

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

1. lw70 <- log(weight/70)

2. Data(data,

3. res=c("cerato"),

4. cova=c("lw70"))

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

1. Model(

2. trans.pr ? = c(1,1,1),

3. cor = c(0,1,0))

4. 


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

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

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

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

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語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

4.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例2

5.R語言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

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

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

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

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