问题描述
我有一组线性混合模型，并创建了一个平均模型.我想绘制该模型对一个因子的两个水平的拟合，该模型包含在平均模型中.一个简单的例子:
I have a set of linear mixed models, and have created an average model. I'd like to plot the model fits for two levels of a factor, included in the average model. A simple example:
library(lme4)
library(MuMIn)
mtcars2 <- mtcars
mtcars2$vs <- factor(mtcars2$vs)
gl <- lmer(mpg ~ am + disp + hp + qsec + (1 | cyl), mtcars2,
REML = FALSE, na.action = 'na.fail')
d <- dredge(gl)
av <- model.avg(d, subset = cumsum(weight) <= 0.95)
summary(av)
Call:
model.avg(object = d, subset = cumsum(weight) <= 0.95)
Component model call:
lme4::lmer(formula = mpg ~ <7 unique rhs>, data = mtcars2, REML = FALSE, na.action = na.fail)
Component models:
df logLik AICc delta weight
13 5 -77.81 167.92 0.00 0.37
123 6 -76.34 168.05 0.13 0.35
134 6 -77.54 170.43 2.51 0.11
1234 7 -76.25 171.16 3.24 0.07
23 5 -79.85 172.00 4.08 0.05
2 4 -81.63 172.75 4.83 0.03
124 6 -78.99 173.34 5.42 0.02
Term codes:
am disp hp qsec
1 2 3 4
Model-averaged coefficients:
(full average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 25.457505 6.467643 6.648016 3.829 0.000129 ***
am 4.103425 1.861593 1.898182 2.162 0.030636 *
hp -0.043829 0.017926 0.018265 2.400 0.016415 *
disp -0.009419 0.011834 0.011983 0.786 0.431821
qsec 0.081973 0.284147 0.292015 0.281 0.778929
(conditional average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 25.45751 6.46764 6.64802 3.829 0.000129 ***
am 4.46519 1.46823 1.51835 2.941 0.003273 **
hp -0.04651 0.01471 0.01515 3.070 0.002140 **
disp -0.01793 0.01068 0.01099 1.632 0.102634
qsec 0.40421 0.51757 0.53873 0.750 0.453075
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Relative variable importance:
hp am disp qsec
Importance: 0.94 0.92 0.53 0.20
N containing models: 5 5 5 3
我想绘制完全平均模型估计的am效果.
I want to plot the effect of am as estimated by the full averaged model.
通常我会使用lsmeans::lsmeans(gl, ~am)或lmerTest::lsmeansLT(gl, 'am')并绘制两组的最小二乘均值及其置信区间.
Normally I would use lsmeans::lsmeans(gl, ~am) or lmerTest::lsmeansLT(gl, 'am') and plot the least squares means for the two groups and their confidence intervals.
如何为普通模型做同样的事情?
How can I do the same for the average model?
推荐答案
(经过一些讨论和进一步的发现，这是经过修订的答案.请注意，我是emmeans软件包的作者.)
(This is a revised answer, after some discussion and further findings. Note that I'm the emmeans package author.)
这似乎有用.
首先，定义 emmeans 程序包所需的方法:
First, define methods needed by the emmeans package:
library(emmeans)
terms.averaging = function(x, ...)
terms(x$formula)
recover_data.averaging = emmeans:::recover_data.lm
### NOTE: still have to provide 'data' argument
emm_basis.averaging = function(object, trms, xlev, grid, ...) {
bhat = coef(object, full = TRUE)
m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
X = model.matrix(trms, m, contrasts.arg = object$contrasts)
V = vcov(object, full = TRUE)
dffun = function(k, dfargs) NA
dfargs = list()
list(X=X, bhat=bhat, nbasis=estimability::all.estble, V=V,
dffun=dffun, dfargs=dfargs, misc=list())
}
因为没有一个，所以需要terms方法.其他两个是从lm对象的现有方法改编而成的.现在有一个陷阱:vcov()调用要求对象具有非NULL "modelList"属性.并且您的av对象失败.但是model.avg帮助页面底部的示例显示了操作方法:
The terms method is needed because there isn't one. The other two are adapted from the existing methods for lm objects. Now there is one catch: the vcov() call requires the object to have a non-NULL "modelList" attribute. And your av object fails. But the examples at the bottom of the help page for model.avg shows what to do:
cs95 = get.models(d, cumsum(weight) <= .95)
AV = model.avg(cs95)
现在，AV具有必需的属性.现在我们得到:
Now, AV has the required attribute. Now we get:
em = emmeans(AV, ~ am, at = list(am = c("0", "1")), data = mtcars)
em
## am emmean SE df asymp.LCL asymp.UCL
## 0 15.42665 2.985460 NA 9.575257 21.27805
## 1 19.53008 1.986149 NA 15.637297 23.42286
pairs(em)
## contrast estimate SE df z.ratio p.value
## 0 - 1 -4.103425 1.861593 NA -2.204 0.0275
请注意，对比结果与模型摘要表中av的估计值和未经调整的SE相匹配.
Note that the contrast result matches the estimate and unadjusted SE for av in the model summary table.
注意:使用coef(..., full = FALSE)和vcov(... full = FALSE)会产生一个非正定协方差矩阵，从而导致EMM的负方差估计.
Note: Using coef(..., full = FALSE) and vcov(... full = FALSE) yielded a non-positive-definite covariance matrix, resulting in negative variance estimates for the EMMs.
而且我警告，尽管这似乎在计算上可行，但这并不意味着答案是正确的！
And I caution that while this seems to work computationally, that does not imply that the answers are right!
这篇关于从平均模型中拟合离散变量的绘图模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！