关于数据挖掘:非线性混合效应-NLME模型对抗哮喘药物茶碱动力学研究附代码数据

全文下载链接：http://tecdat.cn/?p=24074

最近咱们被客户要求撰写对于非线性混合效应 NLME模型的钻研报告，包含一些图形和统计输入。

茶碱数据文件报告来自抗哮喘药物茶碱动力学钻研的数据。给 12 名受试者口服茶碱，而后在接下来的 25 小时外在 11 个工夫点测量血清浓度

head(thdat)

此处，工夫是从抽取样品时开始给药的工夫（h），浓度是测得的茶碱浓度（mg/L），体重是受试者的体重（kg）。

12 名受试者在工夫 0 时承受了 320 mg 茶碱。

让咱们绘制数据，即浓度与工夫的关系：

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

相干视频

**

拓端

，赞14

**

拓端

，赞21

**

拓端

，赞7

**

拓端

，赞17

数据的个体差异

咱们还能够在 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 曲线。

点击题目查阅往期内容

R语言LME4混合效应模型钻研老师的受欢迎水平

左右滑动查看更多

01

02

03

04

依据定义，它们没有思考患者之间的变异性，因而不能提供良好的个体预测。

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))

本文摘选 《 R语言非线性混合效应 NLME模型(固定效应&随机效应)反抗哮喘药物茶碱动力学钻研 》 ，点击“浏览原文”获取全文残缺材料。

点击题目查阅往期内容

R语言用线性混合效应（多程度/档次/嵌套）模型剖析腔调高下与礼貌态度的关系\
R语言LME4混合效应模型钻研老师的受欢迎水平R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model剖析藻类数据实例\
R语言混合线性模型、多层次模型、回归模型剖析学生均匀问题GPA和可视化\
R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例\
R语言用lme4多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据R语言 线性混合效应模型实战案例\
R语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据\
R语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状\
R语言基于copula的贝叶斯分层混合模型的诊断准确性钻研\
R语言建设和可视化混合效应模型mixed effect model\
R语言LME4混合效应模型钻研老师的受欢迎水平\
R语言 线性混合效应模型实战案例\
R语言用Rshiny摸索lme4狭义线性混合模型（GLMM）和线性混合模型（LMM）\
R语言基于copula的贝叶斯分层混合模型的诊断准确性钻研\
R语言如何解决线性混合模型中畸形拟合(Singular fit)的问题\
基于R语言的lmer混合线性回归模型\
R语言用WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型\
R语言分层线性模型案例\
R语言用WinBUGS 软件对学术能力测验（SAT）建设分层模型\
应用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM\
R语言用WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型\
SPSS中的多层（等级）线性模型Multilevel linear models钻研整容手术数据\
用SPSS预计HLM多层（档次）线性模型模型