如何绘制混合模型的结果答案

【问题标题】：How to plot the results of a mixed model如何绘制混合模型的结果
【发布时间】：2019-06-10 07:07:09
【问题描述】：

我在 R 中使用 lme4 来拟合混合模型

lmer(value~status+(1|experiment)))

其中值是连续的，状态（N/D/R）和实验是因素，我得到

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

我想以图形方式表示固定效果评估。然而，这些对象似乎没有绘图功能。有什么方法可以图形化地描绘固定效果？

【问题讨论】：

查看 CRAN 上的 coefplot 或 coefplot2 包。并使用data= 参数来构建您的模型拟合过程...
不要认为 coefplot 适用于混合模型。
哎呀。 arm::coefplot 也不适用于 mer 对象。但是coefplot2::coefplot2 可以。
CRAN 上有 coefplot2 吗？我似乎找不到它。我在平台上使用 R 2.14.1：x86_64-pc-linux-gnu（64 位）。
哎呀哎呀哎呀。它在 r-forge 上。我应该把它推到 CRAN。

标签： r

【解决方案1】：

使用coefplot2（在 r-forge 上）：

从@Thierry 窃取模拟代码：

set.seed(101)
dataset <- expand.grid(experiment = factor(seq_len(10)), 
    status = factor(c("N", "D", "R"), levels = c("N", "D", "R")), 
    reps = seq_len(10))
X <- model.matrix(~status,dataset)
dataset <- transform(dataset, 
    value=rnorm(nrow(dataset), sd = 0.23) +   ## residual
    rnorm(length(levels(experiment)), sd = 0.256)[experiment] +  ## block effects
    X %*% c(2.78,0.205,0.887))  ## fixed effects

拟合模型：

library(lme4)
model <- lmer(value~status+(1|experiment), data = dataset)

剧情：

install.packages("coefplot2",repos="http://r-forge.r-project.org")
library(coefplot2)
coefplot2(model)

编辑：

我经常在使用 R-Forge 构建时遇到问题。如果 R-Forge 构建不工作，这个后备应该可以工作：

install.packages("coefplot2",
  repos="http://www.math.mcmaster.ca/bolker/R",
  type="source")

请注意，coda 依赖项必须已安装。

【讨论】：

感谢您的贡献，本。我看到您模拟数据集，构建模型并使用 coefplot2。但是我似乎无法在 CRAN 中找到 coefplot2。这个包还有其他存储库吗？
是的——见我上面的评论，以及上面代码中引用 r-forge 的（编辑的）install.packages() 命令（我今天是个傻瓜）。它在 r-forge ...
我可以看到 coefplot 2 包的当前状态是：“构建失败”并且无法在 R 2.15.2 上安装它。是否放弃了进一步的开发？而对于哪个 R 版本。可以用吗？
coefplot2 的当前 (2018) 状态是什么？这仍然是您的首选，还是现在 CRAN 中有可靠的软件包，性能同样好或更好？
我使用dotwhisker + broom (CRAN) 或 broom.mixed (Github)

【解决方案2】：

我喜欢系数置信区间图，但考虑一些额外的图来理解固定效应可能会很有用..

从@Thierry 窃取模拟代码：

library(ggplot2)
library(lme4)
library(multcomp)
dataset <- expand.grid(experiment = factor(seq_len(10)), status = factor(c("N", "D", "R"), levels = c("N", "D", "R")), reps = seq_len(10))
dataset$value <- rnorm(nrow(dataset), sd = 0.23) + with(dataset, rnorm(length(levels(experiment)), sd = 0.256)[experiment] + ifelse(status == "D", 0.205, ifelse(status == "R", 0.887, 0))) + 2.78
model <- lmer(value~status+(1|experiment), data = dataset)

看看数据的结构...看起来很平衡..

library(plotrix); sizetree(dataset[,c(1,2)])

跟踪固定效应之间的相关性可能会很有趣，尤其是当您拟合不同的相关性结构时。以下链接提供了一些很酷的代码...

http://hlplab.wordpress.com/2012/03/20/correlation-plot-matrices-using-the-ellipse-library/

my.plotcorr(
matrix(c(1,     .891,   .891,
        .891,   1,      .891,
        .891,   .891,   1), nrow=3)
)

最后，查看 10 个实验的变异性以及实验中“状态”的变异性似乎很有意义。我仍在为此编写代码，因为我在不平衡数据上打破了它，但这个想法是......

My2Boxes(m=4,f1=dataset$experiment,f2=dataset$status,x=dataset$value,color=c("red","yellow","green"))

最后，已经提到的 Piniero 和 Bates (2000) 的书从我略读的少量内容中强烈支持格子。所以你可以试一试。也许像绘制原始数据一样......

lattice::xyplot(value~status | experiment, groups=experiment, data=dataset, type=c('p','r'), auto.key=F)

然后绘制拟合值...

lattice::xyplot(fitted(model)~status | experiment, groups=experiment, data=dataset, type=c('p','r'), auto.key=F)

【讨论】：

【解决方案3】：

这里有一些建议。

library(ggplot2)
library(lme4)
library(multcomp)
# Creating datasets to get same results as question
dataset <- expand.grid(experiment = factor(seq_len(10)),
                       status = factor(c("N", "D", "R"),
                       levels = c("N", "D", "R")),
                       reps = seq_len(10))
dataset$value <- rnorm(nrow(dataset), sd = 0.23) + 
                   with(dataset, rnorm(length(levels(experiment)),
                        sd = 0.256)[experiment] +
                   ifelse(status == "D", 0.205,
                          ifelse(status == "R", 0.887, 0))) +
                   2.78

# Fitting model
model <- lmer(value~status+(1|experiment), data = dataset)

# First possibility
tmp <- as.data.frame(confint(glht(model, mcp(status = "Tukey")))$confint)
tmp$Comparison <- rownames(tmp)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) +
  geom_errorbar() + geom_point()

# Second possibility
tmp <- as.data.frame(confint(glht(model))$confint)
tmp$Comparison <- rownames(tmp)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) +
  geom_errorbar() + geom_point()

# Third possibility
model <- lmer(value ~ 0 + status + (1|experiment), data = dataset)
tmp <- as.data.frame(confint(glht(model))$confint)
tmp$Comparison <- rownames(tmp)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) +
  geom_errorbar() + geom_point()

【讨论】：

您能解释一下glht 在您的代码中的用法吗？我读到它是一个用于测试General Linear Hypotheses 的函数，在这里感觉有点不必要......我也只是用一个更复杂的数据集/模型组合和多个固定效果来测试它，它没有给我Comparison 不再命名...有没有办法让您的代码更通用？

【解决方案4】：

此答案说明了较新的 dotwhisker::dwplot + broom.mixed 解决方案。

在模拟中再添加一个变量：

dataset <- transform(dataset, 
    value=rnorm(nrow(dataset), sd = 0.23) +   ## residual
    rnorm(length(levels(experiment)), sd = 0.256)[experiment] +  ## block effects
        X %*% c(2.78,0.205,0.887),
    var2=rnorm(nrow(dataset)))  ## fixed effects

拟合两个不同的模型：

library(lme4)
model <- lmer(value~status+var2 + (1|experiment), data = dataset)
model2 <- update(model, . ~ . -var2)

绘图：

library(broom.mixed)
library(dotwhisker)
dwplot(list(first=model,second=model2), effects="fixed")+
    geom_vline(xintercept=0, lty=2)

（使用effects="fixed" 只获取固定效果参数，默认删除截距）。

broom.mixed 有许多其他选项。当我想做一些复杂的事情时，我可以使用ggplot + ggstance::geom_pointrangeh (+ position="position_dodgev") 来制作我自己的自定义情节，而不是依赖dotwhisker::dwplot()。

【讨论】：