线性混合模型学习

cococo1015

已于 2024-04-18 15:06:36 修改

阅读量960

点赞数 9

文章标签：学习

于 2024-04-18 15:03:38 首次发布

本文链接：https://blog.csdn.net/cococo1015/article/details/137879091

版权

前言

学习文献（doi: 10.1111/nph.19521）中的数据分析方法（线性混合模型看根性状对实验处理的响应这一部分）。（未系统学习统计方法，因此本文内容不包括原理，只是记录总结流程，及别人得出的使用方式。）

文献数据分析流程：

①利用R包：glmmTMB 构建线性混合模型

glmmTMB包的学习可以参考：广义线性模型glmmTMB包_glmmadmb-CSDN博客；以及R语言官网文件

②利用R包：DHARMa 检验模型是否合适

DHARMa包在网页上还未找到大佬们的中文解读，因此最好的参考就是R包解释的网页或GitHub（https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html）

③将全模型和零模型、简化模型对比进行评估（AIC）

④模型显著性的检验方法（Likelihood ratio test）-----（此步不确定）

⑤利用R包：EMMEANS 进行事后检验（post hoc test），并用Benjamini– Hochberg procedure法矫正p值

首先为什么要构建线性混合模型：可参考R中线性混合效应模型及其应用（一），个人理解：假设对不同木本植物的性状a与实验处理（eg. 遮荫或光照）、生长位置（eg. 不同纬度）之间的相关性构建模型（性状a为连续型变量），但是性状a的变异可能不止是由实验处理的不同引起的，因为每个不同的木本植物本身性状a都不一样，故同一个木本植物的性状a对于不同实验处理或者在不同生长位置的表现并不独立，因此可以对每个木本植物设置随机截距，构建线性混合模型。

此处以实验数据为例：

setwd(……)
df<-read_excel(……)
head(df)
#结果
#ID      species               mycorrhizal_type elevation nutrient_type    RD
#  <chr>   <chr>                 <chr>            <chr>     <chr>         <dbl>
#1 C-FX-G1 Liquidambar formosana AM               high      CK            0.259
#2 C-FX-G2 Liquidambar formosana AM               high      CK            0.297
#3 C-FX-G3 Liquidambar formosana AM               high      CK            0.315
#4 C-FX-G4 Liquidambar formosana AM               high      CK            0.326
#5 C-FX-G5 Liquidambar formosana AM               high      CK            0.276
#6 C-FX-G6 Liquidambar formosana AM               high      CK            0.334

①利用R包：glmmTMB 构建线性混合模型

数据集研究了不同树种（species）在不同海拔（elevation：high/low）以及不同养分水平（nutrient_type：CK为背景值/I为添加养分）的根直径（RD）。

这里将根直径RD作为因变量，species作为随机效应，elevation和nutrient_type为固定效应。建立如下：

mod1<-glmmTMB(RD~elevation*nutrient_type +(1|species),data=df,family = gaussian(link = "log"))

此处的(1|species)即为随机效应的部分，其中”|”左边为1代表固定斜率但对不同树种有随机截距，若表达式写成(1+nutrient_type|species)则说明为随机斜率和截距模型，代表不同树种本身的直径存在差异，并且在养分添加的影响下，他们的直径变化趋势不同。

表达式这里的理解可以参考：R语言，GLMM 模型，lme4包中的 lmer（）的使用_r语言lmer-CSDN博客（虽然不是一个包的，但是方法都差不多）

nutrient_type*elevation也可以写成：nutrient_type + elevation + nutrient_type: elevation，表明是养分添加、海拔、以及它们的交互作用一起作为固定因子。

此处可参考：R语言回归分析_r语言回归分析结果解读-CSDN博客

参数“family=”为残差的分布方式？（这里为参考文献写出的，具体待研究）

（此处发现一个点：在用nlme包的lme()函数时或是lme4/lmerTest包的lmer函数时，在formula这一部分写x~a*b+……（出现交互项），对于a和b是什么类型的数据是没有要求的，都可以写（数值变量或是分类变量），但当使用glmmTMB()函数时，若a和b其中有一项是数值型变量时，写x~a*b+……会报错，此时要么选择a和b均为分类变量，要么使用factor函数将a或b变为分类变量，待研究）

②&③步模型检验和评估

利用DHARMa包对建立的模型进行检验：

library(DHARMa)
sim1<-simulateResiduals(fittedModel = mod1)  #simulateResiduals()函数属于DHARMa包
plot(sim1)

得到如下的结果图：

左边qq图的KS检验p值显著，说明该模型的构建并不是很好

（根据R包说明网页，若出现以下几个情形，则说明是overdispersion：

1) qq图呈现“S”形 & KS test significant；2) Dispersion test significant；3) outlier test significant；4) quantile fit are spread out too far（即右边图中的数值过于分散）

若出现倒S型则说明underdispersion了，此时可以给模型加上”dispformula=”的参数

有时模型还会出现零膨胀现象，可以通过testZeroInflation()检验，通过” ziformula =”参数调节。（具体还不太明白，但相关知识可网页上查询）

调节模型：

mod2<-glmmTMB(RD~elevation*nutrient_type+(1+elevation|species)+(1|species),data=df,ziformula= ~elevation,family = gaussian(link = "log"))

再运行以上代码：

Sim2<-simulateResiduals(fittedModel = mod2)  
plot(sim2)

则没有显著的p值了。

但是用anova()对mod1和mod2进行检验

anova(mod1,mod2)
#结果
#Data: df
#Models:
#mod1: RD ~ elevation * nutrient_type + (1 | species), zi=~0, disp=~1
#mod2: RD ~ elevation * nutrient_type + (1 + elevation | species) + , zi=~elevation, disp=~1
#mod2:     (1 | species), zi=~0, disp=~1
#     Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
#mod1  6 -645.33 -624.15 328.66  -657.33                         
#mod2 11 -639.51 -600.69 330.76  -661.51 4.1877      5     0.5227

p值并不显著，说明两个模型差异并没有很大，且mod2的AIC值反而更大，因此还是选择mod1（不知道对于这部分的解释是否正确）

注：由于此处的mod1和mod2在固定效应部分没有发生变化，而在随机效应部分发生改变，因此这里mod1和mod2应该使用限制性最大似然估计法REML，但这glmmTMB函数中该参数默认TRUE，因此不用额外加。

（相关总结：

北大的R语言教程（37 线性混合模型 | R语言教程）

该教程中提到：

为了比较随机效应相同、固定效应不同的模型，需要两个模型都使用最大似然估计；

；

lme2和lme3的固定效应相同，仅随机效应有差别，当固定效应相同时用anova()函数比较两个模型的随机效应部分，应使用REML估计

lme2和lme3的固定效应相同，方差结构不同且构成嵌套关系，可以用基于REML的似然比检验进行比较：

）

其他关于ML和REML（https://blog.csdn.net/u014182497/article/details/82252456；https://blog.csdn.net/qq_39355550/article/details/81809467；https://zhuanlan.zhihu.com/p/66062778）

④模型显著性的检验方法（Likelihood ratio test）

可以建立零模型：只有随机截距；或建立简化模型，比如将固定效应中的elevation或nutrient_type或它们的交互项elevation*nutrient_type去掉其中一项，进行anova()分析，

例如：

reduced_model1<-glmmTMB(RD~elevation+nutrient_type+(1|species),data=df,family = gaussian(link = "log"))
anova(mod1, reduced_model1)

注：这里只改变了固定效应，所以应该用最大似然法（ML）进行估计，以下稍微修改：

mod1<-glmmTMB(RD~elevation*nutrient_type +(1|species),data=df,family = gaussian(link = "log"),REML=FALSE)
reduced_model1<-glmmTMB(RD~elevation+nutrient_type+(1|species),data=df,family = gaussian(link = "log"),REML=FALSE)
anova(mod1, reduced_model1)
#结果
#Data: df
#Models:
#reduced_model1: RD ~ elevation + nutrient_type + (1 | species), zi=~0, disp=~1
#mod1: RD ~ elevation * nutrient_type + (1 | species), zi=~0, disp=~1
#               Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
#reduced_model1  5 -647.30 -629.65 328.65  -657.30                         
#mod1            6 -645.33 -624.15 328.66  -657.33 0.0319      1     0.8583

p值不显著，说明海拔和养分添加的交互作用对根直径的影响不大。

但如果去掉海拔因素：

reduced_model2<-glmmTMB(RD~nutrient_type+(1|species),data=df,family = gaussian(link = "log"),REML=FALSE)
anova(mod1,reduced_model2)
#结果：
#Data: df
#Models:
#reduced_model2: RD ~ nutrient_type + (1 | species), zi=~0, disp=~1
#mod1: RD ~ elevation * nutrient_type + (1 | species), zi=~0, disp=~1
#               Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)   
#reduced_model2  4 -638.99 -624.88 323.50  -646.99                            
#mod1            6 -645.33 -624.15 328.66  -657.33 10.333      2   0.005705 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

说明海拔对于根直径的影响显著。

关于全模型和零模型、简化模型的对比可能写的不是很详细，只代表对网上其他讲解的理解。看了大部分大佬的解读，一开始构建就要从全模型开始构建，然后通过对AIC或BIC或者残差分析等的评估逐渐简化得到最优模型，因此此步不是特别严谨，希望以后可以深入学习一下。

⑤事后检验（post hoc test）

可以来看看summary(mod1)的结果

关于summary结果的解读可以参考：

R语言|1.重复测量连续数据：混合效应模型分析 - 知乎（网络上还有很多）

解读以上summary的结果，在固定效应部分，只有低海拔的和截距项（即高海拔不添加养分）相比是显著的，但是对于线性模型来说，是会设定一个基准水平，后续所有的水平都会与该基准水平比较，因此线性模型不可以看出不同分类之间的差异是否显著（基线设置可参考：R语言中混合线性模型的实现以及参数解析_the random-effects parameters and the residual var-CSDN博客后部分提到），并且不能知道分类变量总体是否显著，方差分析可以检验变量总体的显著性（理解为对模型summary无法了解分类变量总体是否显著影响结果，如海拔或是否添加养分整体是否影响根直径，只能看到不同项之间的对比，而在构建简化模型少了某个效应后，再用anova方差分析，得到的结果就是分类变量整体的影响了）。例如利用本数据，我想要探究在不添加养分或添加养分的情况下，木本植物的根直径在高海拔和低海拔是否有差异，期望画出参考文献中这样的图：

但是只有summary的结果并不能看出每个具体分组之间的显著性。（这里提一下从别的网页中总结的，summary(mod1)可以看固定效应和随机效应，而方差分析anova(mod1)可以看主效应和交互作用，简单效应：交叉效应中单一变量的作用，参考：如何理解主效应简单主效应和交叉作用？ - 知乎）

灰字部分不是本步骤中想要讲述的重点，但是不知道另外可以写到哪里，就放这儿了。

因此采用emmeans包的emmeans()函数：

emmeans(mod1,pairwise~elevation|nutrient_type)
#结果
#$emmeans
#nutrient_type = CK:
# elevation emmean    SE  df lower.CL upper.CL
# high       -1.20 0.143 246    -1.48   -0.913
# low        -1.25 0.143 246    -1.53   -0.970

#nutrient_type = I:
# elevation emmean    SE  df lower.CL upper.CL
# high       -1.17 0.143 246    -1.45   -0.887
# low        -1.23 0.143 246    -1.52   -0.951

#Results are given on the log (not the response) scale. 
#Confidence level used: 0.95 

#$contrasts
#nutrient_type = CK:
# contrast   estimate     SE  df t.ratio p.value
# high - low   0.0576 0.0263 246   2.190  0.0294

#nutrient_type = I:
# contrast   estimate     SE  df t.ratio p.value
# high - low   0.0643 0.0267 246   2.407  0.0168

根据我想要分析数据的形式：不同养分添加形式下不同海拔之间的根直径差异，即在没有养分添加的条件下，高低海拔的木本植物的根直径有什么差异，以及在养分添加条件下有什么差异，写出如上代码，结果可以看到两个部分，第一部分是关于均值标准误等信息，第二部分是显著性，说明不管是否有养分添加，高低海拔的根直径都是有显著差异的。

参考文献中还提到的Benjamini–Hochberg procedure矫正p值，即在emmeans()函数中有可选参数”adjust=“，可以设置为“adjust=‘BH’”，即：

emmeans(mod1,pairwise~elevation|nutrient_type,adjust="BH")

至此大致过了一遍文献中提到的线性混合模型步骤啦，还是统计小白，一些比较细节的地方还没有搞得很明白，简单记录学习过程和总结，如果有错误希望大家指出！

参考及摘录（附件）

列出比较推荐的几个参考，整个过程很详细：

推荐阅读但很多步骤还未自己尝试过的：

基于R的线性混合效应模型分析 - 简书

R统计绘图-线性混合效应模型详解(理论、模型构建、检验、选择、方差分解及结果可视化)_r 混合效应模型-CSDN博客

R语言中混合线性模型的实现以及参数解析_the random-effects parameters and the residual var-CSDN博客

可以学习的course：

Zuur et al., 2009 Mixed Effects Models and Extensions in Ecology with R（https://link.springer.com/book/10.1007/978-0-387-87458-6）

Linear Mixed-Effects Models Using R（https://link.springer.com/book/10.1007/978-1-4614-3900-4）

以下是在学习过程中记录的参考网页以及摘录总结（顺序较乱，如果以后有时间再重新编辑整理）：见附件。

cococo1015

关注

9
点赞
踩
22

收藏

觉得还不错? 一键收藏
1
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫