3.3随机截距模型分析结果解释及可视化
3.4随机斜率模型构建、检验及可视化
四、线性混合效应模型选择
4.1多模型比较
4.2模型最优子集筛选
各个数据模型都包含一些假设，当数据满足这些假设时，使用这些模型得到的结果才会更准确。比如广义线性模型要满足“线性”和“独立性”等条件。为了适应更多数据类型，开发出了各种模型，可根据数据情况，选择合适的模型。比如，如果数据不满足“线性”，则可考虑广义可加模型；如果不满足“独立性”，则可以考虑混合效应模型(mixedeffectmodel)。
多水平模型同时包含多个分类因子数据，最主要特点就是非独立性，在多个水平上都存在残差，因此适合混合效应模型，总的来说，其思想就是把高水平上的差异估计出来，这就使得残差变小，估计的结果更为可靠。混合效应模型又称随机效应模型(RandomEffectmodel)，分层线性模型，随机系数模型，方差成分模型等。
比如，一份包含不同施肥处理的多土层微生物数据，其中不同施肥处理与不同土壤深度即是嵌套数据(多水平数据)。
如果用线性模型分析不同施肥处理间微生物群落的差异，则没有考虑到不同土壤深度的微生物群落的差异，暗含了假设：不同土壤深度的微生物群落随不同施肥处理变化的截距与斜率都是相同的。不同土壤深度土壤的微生物群落差异就被线性模型归类到误差中去了。
如何解决这一问题呢？其中一种方法就是使用固定效应模型(FixedEffectModel)，即将不同土壤深度的微生物群落结构随不同施肥处理变化的截距差异和斜率差异都估计出来。在土壤深度分类较少时可以这么做，将次级分类水平作为虚拟变量回归，3分类，则估计2个截距差异和斜率差异。当次级分类水平过多时，需要估计的参数太多，用虚拟变量就会消耗太多的自由度，估计结果不可靠，此时即可考虑使用混合效应模型。
线性混合效应模型假设：
2.残差的方差是同质的-非常重要的假设；
3.随机效应的观测是相互独立的；
4.自变量与因变量存在线性关系；
当然，对数据中异常点的检测也是很有必要的。
此套数据是嵌套数据，包含分类因子depth和tillage。嵌套数据中包含多个分类自变量，数据是非独立性的，模型构建时考虑使用混合效应模型。选择pH为因变量。
#设置工作路径#knitr::opts_knit$set(root.dir="D:\\EcoEvoPhylo\\MEM\\LMM")#使用Rmarkdown进行程序运行Sys.setlocale('LC_ALL','C')#Rmarkdown全局设置#options(stringsAsFactors=F)#R中环境变量设置，防止字符型变量转换为因子library(tidyverse)library(rstatix)library(car)library(lme4)library(lmerTest)#模型固定效应系数检验library(multilevelTools)#模型诊断检验library(extraoperators)library(JWileymisc)#模型诊断及APAStyler以表格形式输出模型结果library(effectsize)#计算标准化回归系数library(influence.ME)#检测高影响观测点library(GGally)#2.1导入数据env<-read.csv("env.csv",row.names=1)env%>%head#2.2计算均值与标准误env%>%group_by(condition)%>%get_summary_stats(names(env[4:ncol(env)]),type="mean_se")%>%left_join(env[1:3],by="condition")->datadata%>%sample_n_by(tillage,size=2)
图1|原始数据表，env.csv。
图2|均值-标准误差计算结果，data。
##2.3.2计算有效样本大小nEffective(n=length(unique(env$depth)),k=nrow(env)/length(unique(env$depth)),dv="pH",id="depth",data=env)
图4|effectivesamplesize。对数据提供的有效独立样本数量的估计。N为独立的样本数量，k为每个独立样本的重复测量次数。
##2.3.3因变量分布、极端值探索###将因变量按照分类因子进行均值分解####depthtmp1<-meanDecompose(pH~depth,data=env)tmp1tmp1$`pHbydepth`tmp1$`pHbyresidual`plot(testDistribution(tmp1[["pHbydepth"]]$X,extremevalues="theoretical",ev.perc=.001),varlab="BetweenDepth")plot(testDistribution(tmp1[["pHbyresidual"]]$X,extremevalues="theoretical",ev.perc=.001),varlab="WithinDepth")#图中黑色点表示极端值。####tillagetmp2<-meanDecompose(pH~tillage,data=env)tmp2do.call(cowplot::plot_grid,c(lapply(names(tmp2),function(x){plot(JWileymisc::testDistribution(tmp2[[x]]$X),plot=FALSE,varlab=x)$Density}),ncol=1))
图5|均值分解图，WithinDepth。
##2.3.4散点图可视化数据分布趋势###widthformat转为longformat(env%>%gather(key="variable",value="value",-condition,-tillage,-depth)%>%rstatix::reorder_levels(variable,order=c("pH","OM","OC","Ammonia","Nitrate","AHN","TN","TP","TK","AP","AK"))->data.p)###设置颜色col1<-colorRampPalette(colors=c("#56B1F7","#132B43"),space="Lab")###绘制分面-分组散点图(data.p%>%ggplot2::ggplot(aes(x=factor(tillage,levels=c("CK","WL","NT","PT")),y=value,color=factor(depth,levels=unique(data.p$depth))))+geom_point(position=position_dodge(0.8),#调整site之间的宽度距离。width=0.5,size=1)+scale_color_manual(name="Depth",#values=ggsci::pal_d3("category10")(10),values=col1(3),breaks=c(unique(as.character(data.p$depth))))+facet_wrap(.~variable,ncol=4,scales="free")+labs(x=NULL,y=NULL)->p1)#绘图函数加上括号，赋值给变量的同时，会在屏幕上显示。###输出pdf图到本地ggsave("p1.pdf",p1,device="pdf",width=unit(10,"cm"),height=unit(8,"cm"))
图6|longformat格式数据表，data.p。
图7|可视化pH数据分布，p1.pdf。
混合效应模型是指包含固定效应和随机效应的模型，固定效应部分就是线性模型模式，随机效应可以分为随机截距模型和随机斜率模型两大类。固定效应因子与随机效应因子没有很严格的区分，更多是根据自己的研究目的进行设置，最后模型都要满足线性混合模型的假设以及符合科学性。
此处使用lme4执行线性混合效应模型，(1+depth)表示depth是随机因子，截距随机。lme4包提供了进行线性混合效应模型、广义线性混合效应模型和非线性混合效应模型拟合和分析的功能。nlme包也可用于混合效应模型的建模。随机效应项在函数表达式中表现为(expr|factor)，expr是一个线性或广义线性等模型公式，factor是随机因子。lme4中用不同的formula区分交叉或嵌套数据。
此过程是模型矩阵与分组因子的一种特殊交互，用以确保固定因子模型矩阵对于分组因子各水平都有不同的影响。分组因子效应被建模为未观察到的随机变量，而不是未知的固定参数。随机截距模型的截距的平均值和标准差是要估计的参数，模型将随机效应的任何非零均值合并为固定效应参数。函数默认执行限制性极大似然估计(REML)，REML假设固定效应结构是正确的,在比较具有不同固定效应的模型时，应该使用最大似然(ML)。因为ML不依赖于固定效应的系数。但ML对两个具有嵌套随机结构的模型进行比较时，对方差项的估计量存在偏差。因此，具有相同固定效应的模型，使用REML比较随机效应不同的模型。
模型(固定效应系数和模型比较)显著性检验方法可选(效果从弱到强)：
1）WaldZ-tests，固定效应系数显著性检验，适合大数据集；
2）Waldt-tests(butLMMsneedtobebalancedandnested)；
3）Likelihoodratiotests，只能用于ML拟合模型，可比较固定效应或随机效应不同的两个模型，那个能更好的拟合数据；
4)Kenward-RogerandSatterthwaiteapproximationsfordegreesoffreedom，可比较两个固定效应模型的差异显著性;
5）MCMC，不适合随机斜率模型;
6）parametricbootstrapconfidenceintervals。
当样本量足够时p值获取可以通过LRT或Waldt检验模型，但样本量较少时更推荐使用Kenward-Roger(只能用于REML)，Satterthwaiteapproximations(REML或ML)或parametricbootstrap。
随机截距模型：模型中只包含随机截距；随机斜率模型：模型中包含随机斜率。
#3.2.1构建随机截距模型control<-lmerControl(optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD",xtol_abs=1e-12,ftol_abs=1e-12,check.conv.singular=.makeCC(action="ignore",tol=1e-12)))##lmerTest::lmer()输出结果包含模型系数检验的p值，lme4::lmer()输出结果则不包含。fm1<-lmerTest::lmer(pH~tillage+A+(1|depth),data=env,REML=FALSE,control=control)#3.2.2模型诊断##JWileymisc::modelDiagnostics()fmd1<-modelDiagnostics(fm1,ev.perc=0.001#可设置为0-1，表示理论分布的比例，超出该比例的值被视为极值。)###绘制残差诊断图plot(fmd1,ask=FALSE,nrow=3)###提取极端值-如果有极端值，则可使用update()函数删除该数据点且更新模型。fmd1$extremeValues#此数据没有异常值，所以为null。
图9|随机截距模型检验图。拟合模型之后，还要检验模型是否符合线性混合效应模型的假设。模型检验图的数量根据固定效应和随机效应数量变化。此模型的模型检验组合图中的第一幅图是残差的分布，应该要大致符合正态分布。第二幅图是残差-预测值图可以大致看出数据是否符合“同方差性”，残差值应该围绕0均匀分布，没有明显的分布趋势。最后，蓝色虚线表示预测值残差的第10%和第90%分位数（根据分位数回归模型估计），如果残差的方差是同质的，蓝色虚线应该是平坦且彼此平行的。如图所示，图中残差分布有明显的趋势，不太适合线性模型，可以改换非线性混合模型，或者残差异质性不严重，也可忽略。最后一幅图是depth随机截距的单变量分布。模型敏感性分析是一种评估模型结果对于参数变化或假设的变化的稳定性和影响的方法，模型假设检验就是敏感性分析的重要内容。
##sjPlot包可以对模型假设检验结果进行可视化fmd1.sj<-sjPlot::plot_model(fm1,type="diag")#fmd1.sj[[1]]#残差正态分布及离群值检验。#fmd1.sj[[2]]$depth#随机效应的Q-Q图。#fmd1.sj[[3]]#残差正态分布检验。#fmd1.sj[[4]]#残差方差同质性检验。###拼图并输出到本地(fm1.assumptions<-cowplot::plot_grid(fmd1.sj[[1]],fmd1.sj[[2]]$depth,fmd1.sj[[3]],fmd1.sj[[4]],labels=c("A","B","C","D")))cowplot::save_plot("fm1.assumptions.png",fm1.assumptions)##performance::check_model()也可对模型假设检验结果进行可视化theme_set(theme_classic(base_size=6))#windows太小，绘图会报错，先将ggplot2绘图的字体设置小一点。pdf("check_model.pdf",width=unit(12,"cm"),height=unit(12,"cm"),family="Times")performance::check_model(fm1)dev.off()
图10|随机截距模型假设检验图，fm1.assumptions.png。sjPlot::plot_model()可以对很多模型的假设检验和模型参数进行绘图，更多参数细节可以plot_model查看函数帮助信息。R中进行相同统计分析的函数很多，使用函数可视化时，还是要知道函数使用的什么数据以及对数据进行过怎样的处理。
图11|随机截距模型假设检验图，check_model.pdf。performance::check_model()的模型检验结果更全，图也更美观。
##检验模型自变量是否具有共线性-一般阈值为5或10等car::vif(fm1)#此模型的两个自变量共线性不强。
图12|方差膨胀因子计算结果。vif输出结果包括variance-inflation(VIF,模型所有项df=1),generalizedvariance-inflationfactors(GVIF,df>1)和自由度校正GVIF值。
##计算unadjustedICC和adjustedICCperformance::icc(fm1)#有时也称为variancepartitioncoefficient(VPC)。
图14|随机截距模型检验图。第8个观测点删除后，tillagePT系数变化过大，可使用sigtest进行检验删除该观测点之后，固定效应参数的显著性是否发生变化，看是否需要删除该观测点。
图15|删除某观测点后模型截距参数显著性的变化。列内容依次代表删除某一观测之后，模型截距参数值(Altered.Teststat)、参数检验显著性(Altered.Sig)以及删除前后截距参数检验结果是否发生变化(Changed.Sig)。Altered.Sig的结果会随着函数中test值的正负而发生改变，Changed.Sig不会。输出结果中只需要看Changed.Sig列即可。因此使用R包，一定要知道他的计算原理，不能完全相信输出结果。
###绘制观测点模型影响力图plot(alt.est,which="dfbetas")plot(alt.est,which="cook",sort=TRUE,cutoff=4/36,xlab="Cook′sDistance",ylab="Samples")plot(alt.est,which="pchange")plot(alt.est,which="sigtest")###删除第8观测点，更新模型及重新诊断模型-此处没有需要删除的观测点，此步只是展示如何更新模型。fm1.2<-update(fm1,data=env[-8,])fmd1.2<-modelDiagnostics(fm1.2,ev.perc=.001)plot(fmd1.2,ask=FALSE,nrow=3)fmd1.2$extremeValues
##推荐使用car包检测模型的高影响力点car::influencePlot(fm1)
图17|car::influencePlot()输出的高影响观测值。推荐使用car包检测模型假设。更多线性模型假设检验信息，可参考R统计绘图-多重线性回归(最优子集法特征筛选，leaps)。
图18|car::influencePlot()输出的离群点、高杠杆值和强影响点的信息整合图。图中纵坐标超过+2或小于-2的观测点可被认为是离群点，横坐标超过0.2或0.3的观测点有高杠杆值(通常为预测值的组合)。圆圈大小与影响成比例，圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点。没有通过检验的点将会在图后列出。
summary()没有对随机效应的方差和整体模型进行检测，接下来介绍随机效应和模型检验的结果。anova()函数可以通过LRT检验进行不同模型的比较。pbkrtest包提供函数可进行只针对REML模型的Kenward-Roger法检验。car::Anova()也可对线性混合效应模型进行检测，且方差分析方法和参数设置更多。
##模型固定效应F检验###TypeIIIanovatablewithp-valuesforF-testsbasedonSatterthwaite's.anova(fm1)#与anova(fm1,type="marginal")结果相同。###基于Kenward-Roger法的模型检验只针对REMLfm1.1<-lmer(pH~AP+tillage+(1|depth),data=env,REML=TRUE,control=control)if(requireNamespace("pbkrtest",quietly=TRUE))anova(fm1.1,type=2,ddf="Kenward-Roger")###Anova中的方差分析方法和参数设置更多-AnalysisofDevianceTable(TypeIIWaldchisquaretests)car::Anova(fm1)
图20|模型的Satterthwaite's检验结果。anova函数中只有一个模型时，则会输出模型固定因子的F检验结果，结果显示三个固定因子的方差检验结果都存在统计学意义。
图21|模型的Kenward-Roger法检验结果。Kenward-Roger法只能用于REML模型，当样本量较小时，更推荐使用此方法。
##ranova()-模型随机效应LRT检验ranova(fm1)
图23|模型随机效应检验结果。ranova()对模型随机效应的似然比检验(LRT)结果支持(1|depth)的方差不显著不为0，则表示depth间无差异，这与depth对模型较低的方差解释率是对应的，则可以把考虑去除模型中的depth的随机截距。如果方差显著不为零，此随机效应则应该保留。
JWileymisc、modelsummary、export和sjPlot等包可以对模型进行检验、提取模型结果以及以文字、表格或图形展示出版需要统计信息：统计值、置信区间、p值等。
##modelTest()可以对模型进行检验并以表格形式提取出版需要统计信息：统计值、置信区间、p值等fm1.da<-JWileymisc::modelTest(fm1)fm1.da$FixedEffects#固定效应检验fm1.da$RandomEffects#随机效应检验fm1.da$EffectSizes#固定效应统计值，包括marginalR2、conditionalR2、cohen’sF2和卡方检验结果。fm1.da$OverallModel$Performance#整体模型检验，包括AIC、BIC、MarginalR2和ConditionalR2等值。#names(modelTest(fm1))#输出结果中包含的列表名。#APAStyler(modelTest(fm1))#可用于去除极端值前后的模型结果比较:APAStyler(list(Original=mt,`OutliersRemoved`=mt3))。APAStyler(modelTest(fm1),format=list(FixedEffects="%s,%s(%s;%s)",RandomEffects=c("%s","%s(%s,%s)"),EffectSizes="%s,%s;%s"),digits=3,pcontrol=list(digits=3,stars=FALSE,includeP=TRUE,includeSign=TRUE,dropLeadingZero=TRUE))
图24|模型综合检验结果，modelTest。数据结果包括AIC/BIC等模型评估指标，随机效应方差和模型残差，固定效应统计量等。输出结果缩写代表的含义，可通过modelTest查看函数的帮助文件。AIC(AkaikeInformationCriterion)；BIC(BayesianInformationCriterion)；LL(LogLikelihood)；LLDF(degreesoffreedomforloglikelihood)；Sigma(residualstandarddeviation)；R2(R2varianceaccountedforinthesample)；F2(Cohen’sf2effectsize,calculatedas[R2/(1R2)])；AdjR2(SamplesizeadjustedR2,abetterestimateofpopulationvarianceaccountedfor)；F(modelFtest)；FNumDF(numeratordegreesoffreedomformodelFtest)；FDenDF(denominatordegreesoffreedomformodelFtest)；P(pvalueformodelFtest)。
##tab_model()以html形式输出模型参数检验结果。(sjPlot::tab_model(fm1,collapse.ci=TRUE,p.style="numeric_stars",file="fm1.html"))#以html表格形式输出模型参数，，使用的是Nakagawa'sR2。
图25|模型参数检验结果，fm1.html。数据结果包括AIC/BIC等模型评估指标，随机效应方差和模型残差，固定效应统计量和MarginalR2/ConditionalR2等。不同函数的主要差异在于R2的计算方法，其他参数只是小数点位数设置差异。sjPlot::tab_model()使用performance::r2()为混合效应模型计算Nakagawa'sR2。
##modelsummary()-能直接将模型参数表格输出为word文档,可以方便的进行多个模型的参数提取及比较，还可以并行提取。(modelsummary::modelsummary(list("fm1"=fm1),stars=TRUE,estimate="{estimate}({std.error}){stars}",#展示模型固定效应统计量、标准误和显著性标记。statistic="conf.int",#展示置信区间。output="fm1.docx"))#使用的performance::r2()计算的R2,会与其他函数的提取的MarginalR2/ConditionalR2值不一样。
图26|模型参数检验结果，fm1.docx。数据结果包括AIC/BIC等模型评估指标，随机效应方差和模型残差，固定效应统计量和边际效应等。modelsummary()使用performance::r2()为混合效应模型计算Nakagawa'sR2,会与其他函数的提取的MarginalR2/ConditionalR2值有些微差别。其他参数只是小数点位数设置差异。
#3.3.1模型系数解释coef(fm1)#未标准化模型系数。effectsize(fm1)#标准化回归系数，不是很有必要，实际预测时使用的是未标准化回归系数。fixef(fm1)#提取未标准化固定效应系数。ranef(fm1)#提取随机效应系数。
图27|未标准化模型系数，coef(fm1)。提取系数可以看到，随机截距模型中，随机效应depth有三个不同的截距系数，各水平对应的其它自变量的斜率系数都是一样的。系数具体含义已在fm1.res中提过，这里不再赘述。
图28|标准化固定因子模型系数，effectsize(fm1)。对自变量和因变量同时经过标准化处理后消除了量纲、数量级等差异的影响，然后构建模型得到的就是标准化回归系数，因此，标准化回归系数可以直观比较不同自变量对因变量的作用大小，但不一定具有统计学意义。未标准化回归系数的意义则是，在其它自变量保持不变的情况下，该自变量变化一个单位导致的因变量的变化量。当自变量中存在分类变量时，使用标准化回归系数那一解释系数意义，且在实际预测时，模型使用的是未标准化系数，因此标准化回归系数可用可不用。
#3.3.2Estimatedmarginalmeans(Least-squaresmeans)-在其它因素保持不变的情况下的固定因子水平的均值。emmeans::emmeans(fm1,~tillage)#tillage各水平的pH最小二乘均值。emmeans::emmeans(fm1,~AP)emmeans::emmeans(fm1,~AP|tillage)#两者交互#env%>%group_by(`tillage`)%>%select(pH)%>%rstatix::get_summary_stats()%>%data.frame()
#3.3.3Post-hoc/ContrastAnalysislibrary(knitr)##tillage水平间估计边际均值的两两比较emmeans::emmeans(fm1,pairwise~tillage,adjust="bonferroni")$contrasts%>%tidy()%>%mutate_if(is.numeric,~round(.,4))%>%kable()##tillage与AP交互效应的估计边际均值的两两比较emmeans::emmeans(fm1,pairwise~tillage*AP,adjust="bonferroni")$contrasts%>%tidy()%>%mutate_if(is.numeric,~round(.,4))##估计边际均值事后检验结果绘图emmeans::pwpp(emmeans::emmeans(fm1,pairwise~tillage*AP,adjust="bonferroni"),type="response")+theme_minimal()
图30|tillage水平间估计边际均值的两两比较。
图31|tillage水平间估计边际均值的两两比较。
图32|估计边际均值事后检验结果绘图。
图34|方差分解结果绘图。柱子中的星号表示固定效应参数检验显著。
#3.3.5绘制效应图(effectplot)##绘制随机效应水平间变化：绘制随机截距和斜率与模型总体截距和斜率的差异。lattice::dotplot(ranef(fm1,condVar=T))##绘制固定效应预测效应图。plot(effects::allEffects(fm1))##绘制模型预测效果图。performance::check_predictions(fm1)
图35|随机效应系数散点图。随机截距模型只有截距的随机效应参数。
图36|模型固定效应预测效应图。AP效应图横轴上的小竖线表示数据点。
图37|模型预测效果图。
可参照此输出描述自己的结果。3.4随机斜率模型构建、检验及可视化
理论知识在随机截距模型部分已提及，此处用随机斜率模型演示如何快速得到上述结果。
此处因为fm1使用的ML法，为了方便后续比较，后续模型都使用的ML法，但是ML法可能会低估随机效应的方差，所以模型构建可优先选用REML法。特别当混合效应模型存在嵌套随机效应时，则不应该选用ML算法。
#3.4.1构建随机斜率模型control<-lmerControl(optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD",xtol_abs=1e-12,ftol_abs=1e-12,check.conv.singular=.makeCC(action="ignore",tol=1e-12)))fm3<-lmerTest::lmer(pH~tillage+AP+(1+tillage|depth),data=env,REML=FALSE)fm3.res<-fm3%>%summaryfm3.res#不是所有回归系数都显著。
#3.4.2模型诊断-检验模型是否符合混合效应模型的假设及是否存在异常值。theme_set(theme_classic(base_size=6))#windows太小，绘图会报错，先将ggplot2绘图的字体设置小一点。pdf("fm3.check.model.pdf",width=unit(12,"cm"),height=unit(12,"cm"),family="Times")(fmd3<-performance::check_model(fm3))dev.off()fmd3$INFLUENTIAL#查看是否存在异常点。
图39|随机斜率模型假设检验综合图。混合线性效应模型需要满足的假设，在图中标题处有注明。与随机截距模型不同的时，此处随机斜率模型增加了三个随机斜率。
#3.4.3模型检验##固定效应显著性检验-图表中有关MarginalR2/ConditionalR2的计算，推荐替换为MuMIn::r.squaredGLMM()的计算结果。###随机效应输出的是标准误options(modelsummary_get="all")(modelsummary::modelsummary(list("fm3"=fm3),stars=TRUE,estimate="{estimate}({std.error}){stars}",#展示模型固定效应统计量、标准误和显著性标记。statistic="conf.int",#展示置信区间。output="fm3.docx"))###随机效应输出的Variance，两个函数的输出可以综合起来，选择自己需要的。(sjPlot::tab_model(fm3,collapse.ci=TRUE,p.style="numeric_stars",file="fm3.html"))##随机效应检验ranova(fm3)##模型整体固定效应检验car::Anova(fm3)##模型间比较anova(fm1,fm3)
图40|随机斜率模型检验结果。随机效应检验结果显示，(1+tillage|depth)项的LRT检验仍然不显著，即在此模型中pH在depth的初始值以及对tillage各水平的响应都没显著差异。不用执行随机斜率模型。模型anova检验结果显示tillage和AP固定因子对pH的影响还是显著的。
#3.4.4方差分解##推荐使用glmm.hp()进行方差分解glmm.hp::glmm.hp(fm3)glmm.hp::glmm.hp(fm3)$hierarchical.partitioning##方差分解结果绘图-详细绘图方式随机截距模型已经画过了，这里使用函数的默认绘图函数进行绘制。plot(glmm.hp::glmm.hp(fm3))
图41|随机斜率模型方差分解结果。模型固定效应的边际R2从~0.85到~0.88，提高的并不多，没必要进行随机斜率模型。
图42|绘制随机效应水平间变化。
##绘制固定因子预测效应图。#plot(effects::allEffects(fm3))plot_model(fm3,type="pred")#tillage与AP各一幅图。
图43|固定因子预测效应图。
##绘制固定效应系数图plot_model(fm3,sort.est=TRUE,#按系数值大小排序。show.values=TRUE,#显示系数值和显著性标签。value.offset=.3)
图44|固定效应系数图。
##绘制模型预测效果图。performance::check_predictions(fm3)
图45|模型预测效果图。
##绘制随机斜率模型及置信区间library(ggeffects)###提取estimatedmarginalmeans(predictedvalues)-与predict()的预测是不一样的。fm3.pred<-ggpredict(fm3,terms=c("AP","tillage","depth"),#当连续自变量多于两个，就不好绘图了,就需要考虑分别绘图。type="re",ci.lvl=0.95)fm3.pred###绘制随机效应fm3.pred%>%plot()###plot_model也可选择自变量提取预测值并绘图#plot_model(fm3,type="pred",#terms=c("AP","tillage","depth"))#terms可以选择预测用自变量。
图46|随机斜率模型及置信区间。如图所示，depth个水平pH在tillage各水平的变化速度相似，没必要构建随机斜率模型。sjPlot是一个很强大的模型参数提取和绘图的包，可以参考帮助文档好好学习一下。
##ggplot2自定义绘图###固定效应+模型预测值绘图，查看模型预测效果。(fm3_plot<-ggplot(env,aes(x=AP,y=pH,colour=tillage,shape=depth))+geom_point(alpha=0.8,size=3.5)+theme_bw()+scale_color_manual(values=ggsci::pal_aaas()(4))+scale_shape_manual(values=c(16,17,15))+scale_linetype_manual(values=c(1,2,4))+#设置depth对应不同线型。geom_line(data=cbind(env,pred=predict(fm3)),#合并env与模型预测pH值。aes(y=pred,linetype=depth#可设置根据depth变化线型。),linewidth=0.75)+#addingpredictedlinefrommixedmodel#facet_wrap(~depth,nrow=1)+#也可根据depth绘制分面图。theme(legend.position="right"))ggsave("fm3.pred.pdf",fm3_plot,device="pdf",family="Times")###绘制混合效应模型预测图-可视化效果一般。#(ggplot(fm3.pred)+#geom_line(aes(x=x,y=predicted))+#slope#geom_ribbon(aes(x=x,ymin=predicted-std.error,ymax=predicted+std.error),#fill="lightgrey",alpha=0.5)+#errorband#geom_point(data=env,#addingtherawdata(scaledvalues)#aes(x=AP,y=pH,colour=depth))+#theme_minimal()#)
图47|随机斜率模型预测效果。可根据随机效应depth绘制分面图，图会更易读。
输出结果部分的marginalR2部分可以自己根据glmm.hp的结果修改。
我们在3.1部分构建了多个线性混合效应模型，可以使用AIC、BIC、R2等模型评价参数对模型进行比较和差异显著性检验。
#4.1.1多个模型参数输出-直接根据模型的AIC、BIC和R2等值和固定效应显著性结果选择模型。modelsummary::modelsummary(list("fm1"=fm1,"fm2"=fm2,"fm3"=fm3,"fm4"=fm4,"fm5"=fm5,"fm6"=fm6,"fm7"=fm7,"fm8"=fm8),#stars=TRUE,estimate="{estimate}[{conf.low},{conf.high}]{stars}",statistic=NULL,output="multi-model.docx")
图48|多模型评估参数输出结果，multi-model.docx。图中只展示了部分内容，用于不同模型间比较。一般来说，AIC、BIC和RMSE值小些比较好，R2值大些比较好。但是所有模型都要考虑偏差-方差平衡，模型不能欠拟合也不能过拟合。根据自己的数据，选择合适的即可。
#4.1.2模型差异比较-复杂模型与简单模型相比。##anova()根据AIC的LikelihoodRatioTest进行模型比较。##固定效应或随机效应不同的模型anova都可比较。##当多于2个模型时，模型顺序放置不同，输出的比较结果可能不一样。所以最好模型间都进行两两比较。anova(fm3,fm1)anova(fm8,fm1)
图49|LRT法模型比较结果。随机截距模型(fm1)与随机斜率模型(fm8),以及与随机斜率+随机截距模型(fm3)的差异都不显著。当多于2个模型时，模型顺序放置不同，输出的比较结果可能不一样，所以最好模型间都进行两两比较。或者使用MarkovChainMonteCarlo(MCMC)或parametricbootstrapconfidenceintervals进行模型比较更好。当使用REML拟合混合效应模型时，REML标准中的一项会随着固定效应结构的变化而变化，从而使得固定效应结构不同的模型之间的比较毫无意义。因此，如果要使用LRT来评估固定效应不同的模型的显著性，则必须使用ML来拟合模型。
##4.1.3pbkrtest包进行Kenward-Rogerapproximation检验##比较的模型必须是REML拟合模型，函数会对ML拟合模型重新进行REML拟合。##实际使用时，发现仅能比较固定效应不同的模型。library(pbkrtest)fm9<-update(fm1.1,.~.-AP)(kr<-KRmodcomp(fm1.1,fm9))#模型顺序依次为复杂模型，简单模型。kr%>%summary##4.1.4Satterthwaiteapproach##实际使用时，发现仅能比较固定效应不同的模型。(sa<-SATmodcomp(fm1.1,fm9))sa%>%summary
图50|固定效应不同模型的比较结果。Kenward-Rogerapproximation和Satterthwaiteapproach的检验结果都显示，复杂模型能更好的拟合数据。
##4.1.5Parametricbootstrap##针对此数据集，对仅随机效应不同的模型的比较效果不好。#PBmodcomp(fm1,fm3,nsim=500,cl=2)#fm1与fm3仅随机效应不同。(pb<-PBmodcomp(fm1.1,fm9,nsim=500,cl=2))pb%>%summary
图51|固定效应不同的模型的Parametricbootstrap法比较结果。基于LRT的Parametricbootstrap检验结果也显示复杂模型能更好的拟合数据。
最优子集回归：算法使用所有可能的特征组合来拟合模型。分析者需要应用自己的判断和统计分析来选择最优模型。当特征数多于观测数时，这个方法的效果就不会好。
此处介绍如何使用MuMIn包进行最优子集回归，筛选特征，根据模型筛选准则筛选模型以及构建平均加权模型。
#4.2.1全模型构建full<-lmer(pH~tillage+AP+TN+TP+TK+Nitrate+Ammonia+AHN+AK+OM+OC+(1|depth),data=env,REML=FALSE,na.action="na.fail")#此参数必须设置，否则后续会报错。full%>%summary#4.2.2生成具有全局模型中固定效应项的组合(子集)的模型的模型选择表，具有可选的模型包含规则。dd<-MuMIn::dredge(#全模型full,#模型子集的排序基准。rank="BIC",#也可按照AIC排序#回归系数标准化方法。beta="partial.sd",#设置输出的模型筛选指数。extra=alist(AIC,BIC,"R^2","adjR^2"),#设置某变量包含在所有模型中。fixed="AP",#设置模型筛选方式,只输出包含AP的模型子集。subset=expression(AP))dd%>%summary()
#4.2.3模型筛选方法一：提取最佳拟合模型，BIC值最低的模型。fm_best<-MuMIn::get.models(dd,1)[[1]]fm_best%>%summary()
图52|根据BIC最小值筛选的最优子集随机截距模型。基于LRT的Parametricbootstrap检验结果也显示复杂模型能更好的拟合数据。BIC最小的模型就是随机截距模型fm1。如果设置REML=TRUE，输出的最佳模型仍然是随机截距模型fm1，但是就是固定效应和随机效应的检验值会有些差异，但是趋势是一致的。输出结果Estimate，Std.Error，tvalue,Pr(>|t|):依次表示对应固定效应的的回归系数及其标准误，回归系数检验的t检验统计量及其p值。分类因子作为自变量时，因子中的每个分类水平都会有一个回归系数，第一个分类水平作为哑变量，其回归系数为0。此模型中tillageCK的回归系数为0。
AICc指标可以纠正估计AIC时因样本量较小而产生的偏差。MuMIn包会计算deltaAICc，还可以使用AICcmodavg包中的AICc函数来比较模型。一般来说，如果模型彼此相差在2个AICc单位(deltaAICc)以内，则它们非常相似。在5个单位内，它们非常相似，相差超过10个单位，则AICc较低的模型拟合更好。利用MuMIn包dredge()生成所有自变量可能组合出的模型，并根据BIC排序。然后使用MuMIn::model.avg()利用模型平均的方法获得所有满足筛选条件的模型加权后的最优模型参数估计。
#4.2.4模型筛选方法二：平均加权模型，##生成所有模型统计量表-包含模型F检验结果d2<-MuMIn::dredge(full,#设置模型包含的terms(回归系数)数目(截距除外)#(1,NA)表示模型中至少要有一个回归系数。m.lim=c(1,NA),rank="BIC",extra=list("BIC","AIC","R^2","*"=function(x){s<-summary(x)c(Rsq=s$r.squared,F=s$fstatistic[[1]])}))d2%>%summary()
图53|多模型参数结果描述统计。
##模型筛选方法二：模型平均加权获得最优模型##根据deltaAICc筛选模型:ΔAICc<2的模型被认为具有相同的优良性。###提取deltaAICc<2的所有模型fm_d2<-MuMIn::get.models(d2,delta<2)#fm_d2#输出结果是一个模型列表数据。###根据信息准则进行模型平均##运行dredge()，未设置rank时，最好先使用get.models()提取模型，再传递给model.avg()，否则predict()等函数可能会对平均后模型不起作用。##平均模型涉及因子变量时，需要注意，模型子集可能包含因子的不同分类，目前对此进行平均会产生毫无意义的结果。fm_avg1<-MuMIn::model.avg(fm_d2,beta="partial.sd")fm_avg1%>%summary()
图54|将所有优良模型加权后的最优模型参数估计。平均加权模型中包含的模型子集以及"full"和““subset””平均模型的回归系数检验结果。平均模型系数：包含“full"和”subset”两种形式，‘subset’(‘conditional’)平均表示只平均该变量回归系数存在的模型。
模型比较最好在同类模型中进行，如，不要使用模型评价参数直接比较lm模型与lmm模型。当模型自变量较多需要进行筛选或想进行模型平均时，推荐使用MuMIn包进行最优子集筛选和模型平均。更详细教程可查看R统计绘图-多重线性回归(平均加权模型/最优子集筛选，MuMIn)。更多输出内容见LMM.html文件。
数据及代码文件：
提取码：zkbr
参考资料：
《白话统计》
资料地址：
推荐阅读
机器学习-分类随机森林分析(randomForest模型构建、参数调优、特征变量筛选、模型评估和基础理论等)
R统计绘图-(O)PLS-DA分析(ropls)
R统计绘图-(O)PLS-DA分析(mixOmics)
机器学习-多元分类/回归决策树模型(rpart包)
R统计绘图-物种丰度差异检验及biomarker鉴定(LEfSe分析、ANOVA和K-W检验)
R统计绘图-LEfSe分析及可视化(microbiomeMarker)
R统计绘图-判别分析(LDA/QDA)
R统计绘图-多重线性回归(平均加权模型/最优子集筛选，MuMIn)
R统计绘图-多重线性回归(最优子集法特征筛选，leaps)
R统计绘图-多重线性回归(最优子集筛选/交叉验证，bestglm)
R统计绘图-多元线性分类分析及绘图(mvpart)
R统计绘图|物种组成堆叠柱形图(绝对/相对丰度)
R统计绘图|物种组成冲积图(绝对/相对丰度，ggalluvial)
R统计绘图|物种组成桑基图(不同分类单元微生物丰度及其从属关系，ggalluvial)-免费送书
R统计绘图|物种组成堆叠面积图(绝对/相对丰度，ggalluvial)
R统计绘图|物种组成气泡图和聚类热图(ggplot2/ComplexHeatmap)
R统计绘图-物种组成分析图(饼图、热图、venn图、花瓣图和UpSet图)
R统计绘图-随机森林分类分析及物种丰度差异检验组合图
机器学习-多元分类/回归决策树模型(tree包)
R统计绘图-NMDS、环境因子拟合(线性和非线性)、多元统计(adonis2和ANOSIM)及绘图(双因素自定义图例)
R统计绘图-RDA分析、Mantel检验及绘图
R绘图-RDA排序分析
R统计绘图-VPA(方差分解分析)
R统计绘图-PCA详解1(princomp/principal/rcomp/rda等)
R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程
R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)
R统计-微生物群落结构差异分析及结果解读
R统计绘图-群落生态学指数计算(α/β多样性、群落稳定性和生态位指数)
R统计绘图-PCA分析及绘制双坐标轴双序图
R中进行单因素方差分析并绘图
R统计-多变量单因素参数、非参数检验及多重比较
R统计绘图-带误差棒径向条形图(RadialColumnChart)
R统计绘图-corrplot绘制热图及颜色、字体等细节修改
R数据可视化之美-节点链接图
R统计绘图-rgbif包下载GBIF数据及绘制分布图
R统计-单因素ANOVA/Kruskal-Wallis置换检验
R统计-正态性分布检验[Translation]
R统计-数据正态分布转换[Translation]
R统计-方差齐性检验[Translation]
R统计-Mauchly球形检验[Translation]
R统计绘图-单、双、三因素重复测量方差分析[Translation]
R统计绘图-混合方差分析[Translation]
R统计绘图-协方差分析[Translation]
R统计绘图-One-WayMANOVA
文献管理|Endnote使用小技巧(修改引文格式、合并文献库、删除重复文献和更新文献信息等)
Alignment--本地blast使用详解1-数据库序列检索下载及比对