cómo trazar los resultados de un modelo mixto

Question

que utilizo lme4 en R para ajustar el modelo mixto

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

donde el valor es continua, el estado (N/D/R) y el experimento son factores, y conseguir

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 

me gustaría representar gráficamente la evaluación de efectos fijos. Sin embargo, parece que no hay función de trazado para estos objetos. ¿Hay alguna forma en que pueda representar gráficamente los efectos fijos?

Answer 1

Aquí hay algunas sugerencias.

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

Answer 2

Usando coefplot2 (en r-forge):

el robo de los código de simulación de @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 

modelo Fit:

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

Terreno:

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

edición:

frecuencia me he estado teniendo problemas con la acumulación de I-Forge. Este repliegue debería funcionar si la acumulación R-Forge no está funcionando:

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

Tenga en cuenta que la dependencia coda debe estar ya instalado.

Answer 3

me gustan las parcelas del intervalo de confianza coeficiente, pero puede ser útil tener en cuenta algunas parcelas adicionales para entender los efectos fijos ..

el robo de los códigos de simulación de @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) 

Obtener una mirar la estructura de los datos se ve equilibrada ... ..

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

Puede ser interesante rastrear la correlación entre los efectos fijos, especialmente si se ajustan a diferentes estructuras de correlación. Hay un cierto código fresco proporcionado en el siguiente enlace ...

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

Finalmente parece relevante para mirar la variabilidad a través de los 10 experimentos, así como la variabilidad a través de "status "dentro de los experimentos. Todavía estoy trabajando en el código para esto ya que lo descompongo en datos desequilibrados, pero la idea es ...

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

Finalmente la ya mencionada Piniero y Bates (2000) libro de celosía fuertemente favorecida por lo poco que he desnatada .. por lo que podría dar que un tiro. Tal vez algo así como el trazado de los datos en bruto ...

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

Y entonces el trazado de los valores ajustados ...

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