簡介
茶堿資料
茶堿資料檔案報告來自抗哮喘藥物茶堿動力學研究的資料。給 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）
- the.dat.dta$id==1 ,c("tme)]
- 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
函數将此（非線性）模型拟合到資料
- nls(neatin ~p.me1(psi, time))
- coef(km1)
并繪制預測濃度 f(t,ψ^)
- e. <- dafme(tm=sq(0,40,=.2))
- w.pd1 <- pedct(pk, newaa=wdf)
- 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) - nll <- predct(kmll, ewta=n.f)
- 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 個受試者的觀察濃度：
- 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)
非線性混合效應 (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 是分組變量。
- Data(dta = data,
- grp = id",
- prditors = "time",
- repose = "con")
結構模型是以前使用的一階吸收和線性消除的單室模型。
- molct <- function(pi,id,x) {
- D <- 320
- fe <-D*a/(V*(a-e))*(exp(-e*t)-exp(-a*t))
需要人口參數向量ψpop的結構模型和一些初始值
- Model(modl = moelpt,
- pi = c(a=1,V=20,ke=0.5))
可以定義幾個選擇和運作算法的選項，包括單個參數的估計 (map=TRUE)、Fisher 資訊矩陣的估計和線性化對數似然 (fim=TRUE) 或重要性采樣的對數似然（ll.is=TRUE）。
種子是用于随機數生成器的整數：使用相同的種子多次運作算法可確定結果相同。
- list(map=TRUE,seed=632545)
- mix(model, dat,optns)
可以顯示估計算法的結果摘要
results 還可以使用單個參數估計值
這些單獨的參數估計可用于計算和繪制單獨的預測
- pred(fit1)
- 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)：
- Model(model ,
- psi ,
- trns.par = c(1,0,1))
備注：這裡，ω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)。
- Model(
- trasf = c(1,0,1),
- 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）這個協變量。
- lw70 <- log(weight/70)
- Data(data,
- res=c("cerato"),
- cova=c("lw70"))
協變量模型再次編碼為（行）向量 (0,1,0)，但變換現在對于三個參數編碼為 1
- Model(
- trans.pr = c(1,1,1),
- cor = c(0,1,0))
随機效應之間的相關性
到目前為止，随機效應被認為是不相關的，即矢量-協方差矩陣Ω是一個對角矩陣。
随機效應之間的相關性可以通過輸入參數covari引入，這是一個大小等于模型中參數數量的方形矩陣，給出了模型的方差-協方差結構。1s對應于估計的方差（在對角線上）或協方差（非對角線元素）。矩陣Ω的結構應該是塊狀的。
例如，考慮一個模型，其中ka在人群中是固定的，即ωka=0（是以對所有i來說kai=0），而log(V)和log(ke)是相關的，即ηV和ηke)是相關的。
- Model(
- covai = t(c(0,1,0)),
- covain = matrix(c(0,0,0,0,1,1,0,1,1),nrow=3))