R数据分析：用lme4包拟合线性和非线性混合效应模型
快一个月没更新文章啦，今天收到好几个粉丝的催更私信，好的吧，实在对不住大家期待的眼神，看样子不能再拖啦，想想写啥好呢，大家咨询比较多的，混合模型算一个，今天就继续给大家写写混合模型如何做吧。
混合模型一般都可以用lme4这个包解决，lme4既可以做线性混合模型，也可以做广义线性混合模型还可以做非线性混合模型，大家有需要可以只研究这一个包就行。
所谓混合模型就是既有固定效应又有随机效应的模型：
“mixedeffects”, denotes a model that incorporates both fixed- and random-effects terms in a linear predictor expression from which the conditional mean of the response can be evaluated
第一部分 线性混合模型
直接上例子，数据是来自一篇研究睡眠剥夺的文献，整个数据大概长下图这样，其中我们的受试者在day0的时候可以睡到自然醒，在之后的日子里所有的受试者就只能睡3个小时了，我们的响应变量是Reaction，就是对受试者做的测验的响应时间，我现在关心睡眠剥夺后，响应时间的变化情况：
对于这么一个纵向数据，
我们来捋一捋：我们只有18个人受试者，每个受试者随访10次，我们需要明白的是，此时我们的每一次测量是嵌套在人的水平上的，我们可以认为，不同人自己的10次测量是有强烈的相关性的，而不同人之间的这种关系又不一定是相同的。
直观一点，我们可以画出来每一天所有人响应时间和睡眠剥夺的变化，画出来就是下图：
可以看到我们上面的这个大图是由很多个小图组成的，每一个小图中横轴就是睡眠剥夺的时间，纵轴是反应时间。每个小图就代表着我们要研究的睡眠剥夺和反应时间的关系（具体到人），但是我们也应该注意到这种关系在不同的人上是不同的，体现在：关系的斜率不同和截距不同。（这个关系的不同可以很明显的在图中看出来）
所以我们就可以拟合一个带有随机效应的混合模型：
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)运行代码后得到下面的结果：
结果中有随机效应的标准差和固定效应的β估计，我得到的截距是251.4，斜率是10.5，这两个系数就是我们研究的总体关系的表示，通常需要在文献中汇报，就意味着睡眠不剥夺的时候人的反应时间是251.4，而睡眠每剥夺多一天反应时间增长10.5。
上面这个是最简单的混合模型。我们继续看：
lme4包高水平设置介绍
- 混合模型公式
对于一个常见的混合模型，我们可以在lme4包中写出来如下差不多的混合模型公式：
resp ~ FEexpr + (REexpr1 | factor1) + (REexpr2 | factor2) + ...这个公式中FEexpr就是固定效应，(REexpr1 | factor1)and(REexpr2 | factor2)都是随机效应，理论上你可以弄很多个随机效应但是实际操作中我们不会关心那么多。
- 理解混合模型公式
我们看到每一个随机效应在公式中的表达都是(expr | factor)这样的。竖杠前面的expr就是一个常规的回归公式，竖杠后面的factor就是一个常规的因子，你可以把竖杠想象成回归公式和因子的交互：
One way to think about the vertical bar operator is as a special kind of interaction betweenthe model matrix and the grouping factor。This interaction ensures that the columns of themodel matrix have different effects for each level of the grouping factor.
这种交互的意思就是在因子的不同水平，我们的回归是不一样的，这也正好和我们前面的解释相对应，就是在不同的人的水平睡眠剥夺和响应时间的关系不一样。
写到这，希望大家能记住下面这张表：
这个表就给我们展示了常见的随机效应的设置，比如(1 | g)，就是说在因子g的不同水平，我们响应变量的截距都不一样。表中的第二行有个offset，表示没有固定效应。如果我们的数据是一个三层嵌套数据，我们可以用第三行的设定来表示随机截距；如果你的数据没有直接嵌套但是在g1和g2的不同水平上存在相关，那么可以用第四行的设定，这个在项目反应理论中比较常见。
在lme4中，默认认为同一个模型的截距和斜率是存在相关的，如果你确定截距和斜率无关那么设定随机效应的时候就可以用两个竖杠，或者把截距和斜率分开来写，就是说(x || g)和x +(1 | g) + (0 + x | g)表达的随机效应都是一样的。
比如如果我认为睡眠剥夺和反应时间随机效应的截距和斜率无关，我便可以做如下设定：
fm2 <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)#截距和斜率无关的设定有时候我们拟合一个后又想尝试对模型进行改变，但又不想重写，此时就可以直接对相似的模型基础上进行更新：
- 模型的更新
比如我想在fm1的基础上去掉随机斜率只留随机截距，我就可以用updata写出如下代码：
fm3 <- update(fm1, . ~ . - (Days | Subject) + (1 | Subject))#模型的更新到底哪一个模型更好呢？
可以用anova方法进行模型间的比较：
anova(fm1, fm2, fm3)运行代码会输出比较的结果：
其中，从模型比较的结果可以看出，给模型增加一个截距和斜率无关的随机效应相比会使得模型的deviance变小，进一步将随机效应设定为相关，并不能够显著地减小deviance，从而我们就可以知道fm2才是对数据拟合最好的模型。
第二部分 非线性混合模型
非线性混合模型就是通过一个连接函数将线性模型进行拓展，并且同时再考虑随机效应的模型。
The fixed-effects parameters describe the general patterns of the data and random-effects parameters describe specific clusters. If the model is nonlinear in the parameters,it is called a nonlinear mixed-effects model (Davidian &Giltinan, 2003)
非线性混合模型常常在生物制药领域的分析中会用到，因为很多剂量反应并不是线性的，如果这个时候数据再有嵌套结构，那么就需要考虑非线性混合模型了。
看下面这个图，这个图描绘了不同人用了茶碱过后的反应，时间是横轴，残留是纵轴，和开篇线性模型中睡眠剥夺和反应时间的例子一样，我们把每个人的关系都做了图出来，不过从图中可以明显看出这种关系并不是简单线性的。
其实这种不是线性的关系存在的情况很多。
比如渐进回归：
再比如逻辑增长：
此时我们要注意到像这两非线性关系模型的参数都不是简单的一个斜率加个截距了。都有φ1，φ2，φ3三个额外参数。
这儿先给大家写一个逻辑增长的实际例子：我现在有一个关于树木周径的数据集，每棵树随访了7次，每次随访记录数的年龄age，和周径，我现在想研究在所有树木中时间和周径的关系。
很自然，我们可以想到不同的树这个关系应该是不一样的，我们想探求的一定是考虑了树水平的变异之后的总体关系，所以不妨先画出来每个树的关系：
从图中可以看到我们总共有5棵树，基本关系是一致的，但存在些许变异相关（所以考虑混合模型），而且这个关系并不是线性的（时间越大周径基本不改变了），所以我们应该考虑非线性的混合模型。
具体地，我们可以用nlmer方法来拟合非线性混合模型，方法参数包括3部分：首先是响应变量，然后是非线性函数，然后是混合效应公式：
The formula argument fornlmeris in three parts: the response, the nonlinear model function depending on covariates and a set of nonlinear model (nm) parameters, and the mixed-effects formula.
比如对我们的数据我就可以写出如下SSlogis方法的代码：
print(nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym | Tree, Orange, start = c(Asym = 200, xmid = 770, scal = 120)), corr = FALSE)此时我们选择的非线性函数是逻辑增长函数SSlogis，刚刚给大家解释了这个函数是有3个参数的，在上面的代码中，age是我们的预测变量，Asym, xmid, scal分别是额外的三个参数（之前的逻辑增长的式子和Asym/(1+exp((xmid-input)/scal))等同）：
进一步，拟合逻辑增长是要我们给出这些参数的初始值的，然后从初始值通过梯度下降寻找各个参数的最优解：
SSlogis has an attribute called "initial", which is a function that nls can call to compute reasonable starting values for fitting a logistic function to the input data.
所以我们看到代码中都给出了响应参数的初始值。
运行上面代码后输出如下结果：
我们可以看到结果中的固定效应里面有Asym，Xmid，scal参数的估计结果。
那么这些参数如何解释呢？
留个悬念我们下期再更。关注关注关注，嘿嘿。
小结
很久没更新了，今天给大家写了如何用lme4做混合模型，包括线性和非线性的例子，感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，请关注后私信回复“数据链接”获取所有数据和本人收集的学习资料。如果对您有用请先收藏，再点赞转发。
也欢迎大家的意见和建议，大家想了解什么统计方法都可以在文章下留言，说不定我看见了就会给你写教程哦，另咨询代做请私信。
R数据分析：用lme4包拟合线性和非线性混合效应模型相关推荐
- R语言非线性混合效应 NLME模型(固定效应随机效应)对抗哮喘药物茶碱动力学研究
最近我们被客户要求撰写关于非线性混合效应 NLME模型的研究报告,包括一些图形和统计输出. 相关视频:线性混合效应模型(LMM,Linear Mixed Models)和R语言实现 线性混合效应模型( ...
- lme4包中的glmer()函数来拟合混合效应模型，并进行逻辑回归分析。
在R语言中,我们可以使用lme4包中的glmer()函数来拟合混合效应模型,并进行逻辑回归分析.glmer()函数的用法与glm()函数类似,但是它能够处理具有随机效应的数据. 下面是一个示例代码,使 ...
- R语言 面板数据分析 plm包实现(一) ——LSDV和固定效应模型
系列文章 R做面板数据分析:R语言 面板数据分析 plm包实现(一) --LSDV和固定效应模型 如果想看随机效应模型怎么做,参见这篇文章 R语言 面板数据分析 plm包实现(二)--随机效应模型 如 ...
- R语言使用glmnet包拟合lasso-cox回归模型（包含生存时间和结果标签）、使用lasso-cox模型进行特征筛选、使用sapply函数对特征数据进行标准化z-score
R语言使用glmnet包拟合lasso-cox回归模型(包含生存时间和结果标签).使用lasso-cox模型进行特征筛选.使用sapply函数对特征数据进行标准化z-score 目录
- R语言使用glmnet包拟合lasso-cox回归模型（生存时间和结果标签）、lasso-cox模型进行特征筛选、plot函数可视化cv.glmnet模型获得的最佳lambda曲线位置及其1个标准差线
R语言使用glmnet包拟合lasso-cox回归模型(包含生存时间和结果标签).使用lasso-cox模型进行特征筛选.plot函数可视化cv.glmnet模型获得的最佳lambda曲线位置及其1个 ...
- R语言使用rms包拟合cph生存分析模型（包含生存时间和结果标签）、绘制不同生存时间节点的列线图nomogram（例如，6个月生存风险、12个月生存风险等）、使用逐步回归筛选最佳的cox回归模型
R语言使用rms包拟合cph生存分析模型(包含生存时间和结果标签).绘制不同生存时间节点的列线图nomogram(例如,6个月生存风险.12个月生存风险等).使用逐步回归筛选最佳的cox回归模型 目录
- R语言使用xgboost包拟合xgboost回归模型、使用predict函数和训练好的模型进行预测推理、计算回归模型的评估指标MAE、MSE、RMSE
R语言使用xgboost包拟合xgboost回归模型.使用predict函数和训练好的模型进行预测推理.计算回归模型的评估指标MAE.MSE.RMSE 目录
- R语言线性混合效应模型（固定效应随机效应）和交互可视化3案例
最近我们被客户要求撰写关于线性混合效应模型的研究报告,包括一些图形和统计输出. 视频:线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例 线性混合效应模型(LMM,Li ...
- R语言线性混合效应模型不同类型的比较
在R语言中,可以使用lme4包中的lmer函数进行线性混合效应模型的分析.在分析不同类型的线性混合效应模型时,我们可以在模型公式中通过更改因变量.自变量.随机效应等参数来创建不同的模型. 下面是一个简 ...
- R语言使用epiDisplay包的lroc函数可视化logistic回归模型的ROC曲线并输出诊断表、输出灵敏度、1-特异度、AUC值等、设置lwd参数自定义ROC曲线线条的粗细（宽度）
R语言使用epiDisplay包的lroc函数可视化logistic回归模型的ROC曲线并输出诊断表(diagnostic table).输出灵敏度.1-特异度.AUC值等.设置lwd参数自定义ROC ...
最新文章
- 绩效真的重要吗？绩效管理系统有哪些？
- 心理所发表关于神经科学研究可信度的评论文章
- 让jquery中的load不缓存方法
- BUUCTF-[网鼎杯 2020 青龙组]singal——angr学习记录
- R语言聚类算法的应用实例
- Qt Creator开发基于小部件的应用程序
- 前端工程构建工具——Yeoman
- presto领读 查询引擎翻译
- Exp3 免杀原理与实践 20164314
- Username is not in the sudoers file. This incident will be reported
- 单例模式简单示例与优化
- lisp医院化验系统_医院LIS系统解决方案
- 狸窝全能视频转换器功能介绍
- 信息安全管理的效益分析
- python爬虫之获取谷歌浏览器所有cookie
- 华硕电脑连接不上wifi_WiFi6 | ASUS中国
- Android实现一键开启自由窗口、分屏、画中画模式——分屏模式
- 面试经历（某大型机器人少儿编程培训机构）
- 关于Oracle parallel（并行）的几个基本常识
- 工业大数据应用技术国家工程实验室