R语言用Rshiny探索lme4广义线性混合模型(GLMM)和线性混合模型(LMM)

随着软件包的进步,使用广义线性混合模型(GLMM)和线性混合模型(LMM)变得越来越容易。

由Kaizong Ye,Liao Bao撰写

最近我们被客户要求撰写关于混合模型的研究报告。由于我们发现自己在工作中越来越多地使用这些模型,我们开发了一套R shiny工具来简化和加速与对象交互的lme4常见任务。 

shiny的应用程序和演示

演示此应用程序功能的最简单方法是使用Shiny应用程序,在此处启动一些指标以帮助探索模型。 

×

广义线性混合模型GLMMGeneralized Linear Mixed Model),是广义线性模型GLM 和线性混淆模型LMM 的扩展形式,于二十世纪九十年代被提出。GLMM因其借鉴了混合模型的思想,其在处理纵向数据重复测量资料)时,被认为具有独特的优势。GLMM不仅擅长处理重复测量资料,还可以用于任何层次结构的数据(因为本质上又是多水平模型)。

提到GLMM,有必要先介绍几个容易混淆的概念:GLM、LMM、MLM、GMM 和GEE。


相关模型简介

广义线性模型 GLM

广义线性模型GLM,是大家经常接触的概念了,比如经典的Logistic模型。GLM是普通线性模型的扩展形式,由于普通线性回归的因变量必须服从正态分布,而实际问题中经常会遇到分类问题或计数问题的建模,GLM采用连接函数(Link Function),将因变量的分布进行了扩展,使得因变量只要服从指数分布族即可(如正态分布,二项分布,泊松分布,多项分布等)。

GLM 可以分解为 Random Component、System Component 和 Link Function 三个部分。Random Component 为残差部分,取决于因变量的分布;System Component 为预测部分,又称 linear predictor,是拟合的关键;Link Function 为连接变化函数,用于将指数分布族转化成正态分布,或者说,对预测结果进行非线性映射(建立 linear predictor与 label 之间的变换关系),是LM成长为GLM的关键环节。

需要强调的是,link function 是从 label 映射到 linear predictor的过程,link function的反函数称为响应函数 response function。响应函数 把 linear predictor 直接映射到了预测目标 label。较常见的link function如 logit函数(又称log-odds);较常用的响应函数如 logistic(又称sigmoid,是二分类中的相应函数)和 softmax(是sigmoid的扩展形式,用于多分类问题),这两个都是 logit 的反函数。

以 Logistic为例:

​因变量为Bernoulli Distribution也就是对二分类问题建模,因变量为Binomial Distribution也就是对多分类问题建模,因变量为Poisson Distribution也就是对计数问题建模(注意区分计数问题和多分类问题)。



在第一个选项卡上,该函数显示用户选择的数据的预测区间。该函数通过从固定效应和随机效应项的模拟分布中抽样并组合这些模拟估计来快速计算预测区间,以产生每个观察的预测分布。

在下一个选项卡上,固定效应和组级效果的分布在置信区间图上显示。这些对于诊断非常有用,并提供了检查各种参数的相对大小的方法。


视频

线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例

探索见解

去bilibili观看

探索更多视频

在第三个标签上有一些方便的方法,显示效果的影响或程度predictInterval。对于每种情况,最多12个,在所选数据类型中,用户可以查看更改固定效应的影响。这允许用户比较变量之间的效果大小,以及相同数据之间的模型之间的效果大小。

预测

预测像这样。

predict(m1, newdata = InstEval[1:10, ])
#>        1        2        3        4        5        6        7        8
#> 3.146336 3.165211 3.398499 3.114248 3.320686 3.252670 4.180896 3.845218
#>        9       10
#> 3.779336 3.331012

预测lmglm

predInte(m1, newdata = Eval[1:10, ], n.sims = 500, level = 0.9,
#>         fit      lwr      upr
#> 1  3.074148 1.112255 4.903116
#> 2  3.243587 1.271725 5.200187
#> 3  3.529055 1.409372 5.304214
#> 4  3.072788 1.079944 5.142912
#> 5  3.395598 1.268169 5.327549
#> 6  3.262092 1.333713 5.304931

预测区间较慢,因为它是模拟计算。

可视化

可视化检查对象的功能。最简单的是得到固定和随机效应参数的后验分布。

 head(Sim)
#>          term        mean      median         sd
#> 1 (Intercept)  3.22673524  3.22793168 0.01798444
#> 2    service1 -0.07331857 -0.07482390 0.01304097
#> 3   lectage.L -0.18419526 -0.18451731 0.01726253
#> 4   lectage.Q  0.02287717  0.02187172 0.01328641
#> 5   lectage.C -0.02282755 -0.02117014 0.01324410

我们还可以快速制作随机效应的图:

head(Sims)
#>   groupFctr groupID        term        mean      median        sd
#> 1         s       1 (Intercept)  0.15317316  0.11665654 0.3255914
#> 2         s       2 (Intercept) -0.08744824 -0.03964493 0.2940082
#> 3         s       3 (Intercept)  0.29063126  0.30065450 0.2882751
#> 4         s       4 (Intercept)  0.26176515  0.26428522 0.2972536
#> 5         s       5 (Intercept)  0.06069458  0.06518977 0.3105805
plotR((m1, n.sims = 100), stat = 'median', sd = TRUE

有时,随机效应可能难以解释

 Rank(m1, groupFctr = "d")
head(ranks)
#>      d (Intercept) (Intercept)_var       ER pctER
#> 1 1866   1.2553613     0.012755634 1123.806   100
#> 2 1258   1.1674852     0.034291228 1115.766    99
#> 3  240   1.0933372     0.008761218 1115.090    99
#> 4   79   1.0998653     0.023095979 1112.315    99
#> 5  676   1.0169070     0.026562174 1101.553    98
#> 6   66   0.9568607     0.008602823 1098.049    97

效果模拟

解释LMM和GLMM模型的结果很困难,尤其是不同参数对预测结果的相对影响。

 
impact(m1, Eval[7, ], groupFctr = "d", breaks = 5,
n.sims = 300, level = 0.9)
 
#>   case bin   AvgFit     AvgFitSE nobs
#> 1    1   1 2.787033 2.801368e-04  193
#> 2    1   2 3.260565 5.389196e-05  240
#> 3    1   3 3.561137 5.976653e-05  254
#> 4    1   4 3.840941 6.266748e-05  265
#> 5    1   5 4.235376 1.881360e-04  176

结果表明yhat根据我们提供的newdata在组因子系数的大小方面,从第一个到第五个分位数的变化。

ggplot(impSim, aes(x = factor(bin), y = AvgFit, ymin = AvgFit - 1.96*AvgFitSE,
ymax = AvgFit + 1.96*AvgFitSE)) +

可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds