原文链接：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）

 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 模型：

pk.od <- function(pi, t){
  D  <- 320
  ka 
  V  
  ke 
  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 的最小二乘预计定义为

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

每个个体预测浓度 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))

最受欢迎的见解

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 增长回归的 HAR-RV 模型预测 GDP 增长 ”)

9. 应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM