Как построить результаты смешанной модели

Я использую lme4 в R для соответствия смешанной модели

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

Я хотел бы графически представить оценку фиксированных эффектов. Однако, похоже, для этих объектов нет функции построения графика. Есть ли способ графически изобразить фиксированные эффекты?

Ответ 1

Вот несколько советов.

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()

Ответ 2

Используя 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

Fit model:

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

Plot:

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 уже установлена.

Ответ 3

Мне нравятся графики доверительного интервала коэффициентов, но может быть полезно рассмотреть некоторые дополнительные графики, чтобы понять фиксированные эффекты.

Кража кода симуляции из @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"))

Наконец, уже упоминавшаяся книга Пиньеро и Бейтса (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)

Ответ 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)

(использование effect effects="fixed" возвращает нам только параметры с фиксированным эффектом, по умолчанию отбрасывая перехват).

broom.mixed есть много других вариантов. Когда я хочу сделать что-то сложное, я могу использовать ggplot + ggstance::geom_pointrangeh (+ position="position_dodgev"), чтобы создать свой собственный график, а не полагаться на dotwhisker::dwplot().