R语言对数线性模型loglm函数_R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）...

逃离我推掉我的手 2023-01-05 03:57 139阅读 0赞

## 原文链接： ##

R语言用Rshiny探索lme4广义线性混合模型（GLMM）和线性混合模型（LMM）tecdat.cn

![6bf6181082885556f65ebd1acdd638dc.png][]

##  ##

随着`lme4`包装的进步，使用广义线性混合模型（GLMM）和线性混合模型（LMM）变得越来越容易。由于我们发现自己在工作中越来越多地使用这些模型，我们开发了一套R shiny工具来简化和加速与对象交互的`lme4`常见任务。

##  ##

## shiny的应用程序和演示 ##

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

![eccdd4b02238d0ea43e884b1c023eab7.png][]

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

![52a6a2bcc55025659a6d457e022b6b9a.png][]

在下一个选项卡上，固定效应和组级效果的分布在置信区间图上显示。这些对于诊断非常有用，并提供了检查各种参数的相对大小的方法。这个标签使得利用四个相关的功能`merTools`：`FEsim`，`plotFEsim`，`REsim`和`plotREsim`它们可将自己的使用。

![11e5d7714aa2917bcd6690cef1d92acd.png][]

在第三个标签上有一些方便的方法，通过利用的力量来显示效果的影响或程度`predictInterval`。对于每种情况，最多12个，在所选数据类型中，用户可以查看更改固定效果或其中一个分组级别术语的影响。使用该`REimpact`函数，如果所有其他情况保持相等，则使用模型的预测模拟每个案例，但是通过固定效应或随机效应项的分布来移动观察。这是在因变量的比例上绘制的，这允许用户比较变量之间的效果大小，以及相同数据之间的模型之间的效果大小。

## 预测 ##

标准预测看起来像这样。

<code>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</code>

随着`predictInterval`我们获得更喜欢所生产的标准对象的预测`lm`和`glm`：

<code>\#predictInterval(m1, newdata = InstEval\[1:10, \]) \# all other parameters are optional predictInterval(m1, newdata = InstEval\[1:10, \], n.sims = 500, level = 0.9, stat = 'median') \#> 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 \#> 7 4.215371 2.136654 6.078790 \#> 8 3.816399 1.860071 5.769248 \#> 9 3.811090 1.697161 5.775237 \#> 10 3.337685 1.417322 5.341484</code>

请注意，`predictInterval`它较慢，因为它是计算模拟。它还可以将所有模拟`yhat`值作为属性返回到预测对象本身。

`predictInterval`大量使用包中的`sim`函数`arm`来绘制模型参数的分布。然后，它将这些模拟值组合起来，`yhat`为每个观测值创建分布。

## 绘制 ##

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

<code>feSims <- FEsim(m1, n.sims = 100) head(feSims) \#> 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 \#> 6 lectage^4 -0.01940499 -0.02041036 0.01196718</code>

我们也可以这样绘制：

<code>plotFEsim(FEsim(m1, n.sims = 100), level = 0.9, stat = 'median', intercept = FALSE)</code>

![584b1814851b5b5fa9df78187386bb7f.png][]

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

<code>reSims <- REsim(m1, n.sims = 100) head(reSims) \#> 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 \#> 6 s 6 (Intercept) 0.08055309 0.05872426 0.2182059</code>

<code>plotREsim(REsim(m1, n.sims = 100), stat = 'median', sd = TRUE)</code>

![bfe495729b2052441e7daaa6a678f4b3.png][]

有时，随机效应可能难以解释，并且并非所有这些都与零有意义地不同

<code>ranks <- expectedRank(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</code>

## 效果模拟 ##

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

<code>impSim <- REimpact(m1, InstEval\[7, \], groupFctr = "d", breaks = 5, n.sims = 300, level = 0.9) impSim \#> 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</code>

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

<code>library(ggplot2) ggplot(impSim, aes(x = factor(bin), y = AvgFit, ymin = AvgFit - 1.96\*AvgFitSE, ymax = AvgFit + 1.96\*AvgFitSE)) + geom\_pointrange() + theme\_bw() + labs(x = "Bin of \`d\` term", y = "Predicted Fit")</code>

![0bc9ef0ae5def4980bec03bec76719f5.png][]

## 非常感谢您阅读本文，有任何问题请在下面留言！ ##

[6bf6181082885556f65ebd1acdd638dc.png]: /images/20221119/a7f8120b8dd14b839994c004614f38c1.png
[eccdd4b02238d0ea43e884b1c023eab7.png]: https://img-blog.csdnimg.cn/img_convert/eccdd4b02238d0ea43e884b1c023eab7.png
[52a6a2bcc55025659a6d457e022b6b9a.png]: https://img-blog.csdnimg.cn/img_convert/52a6a2bcc55025659a6d457e022b6b9a.png
[11e5d7714aa2917bcd6690cef1d92acd.png]: https://img-blog.csdnimg.cn/img_convert/11e5d7714aa2917bcd6690cef1d92acd.png
[584b1814851b5b5fa9df78187386bb7f.png]: https://img-blog.csdnimg.cn/img_convert/584b1814851b5b5fa9df78187386bb7f.png
[bfe495729b2052441e7daaa6a678f4b3.png]: https://img-blog.csdnimg.cn/img_convert/bfe495729b2052441e7daaa6a678f4b3.png
[0bc9ef0ae5def4980bec03bec76719f5.png]: https://img-blog.csdnimg.cn/img_convert/0bc9ef0ae5def4980bec03bec76719f5.png