mixtos - Cómo trazar los resultados de un modelo mixto.

modelos lineales mixtos en r

Question

Uso lme4 en R para ajustar el modelo mixto

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

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

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 manera en que pueda representar gráficamente los efectos fijos?

Answer 1

Aqui 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

Me gustan los gráficos de intervalo de confianza de coeficientes, pero puede ser útil considerar algunos gráficos adicionales para comprender los efectos fijos.

Robando el código 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)

Eche un vistazo a la estructura de los datos ... se ve equilibrado ..

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

Puede ser interesante hacer un seguimiento de la correlación entre los efectos fijos, especialmente si ajusta diferentes estructuras de correlación. Hay un buen código 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 observar la variabilidad en los 10 experimentos, así como la variabilidad en el "estado" dentro de los experimentos. Todavía estoy trabajando en el código para esto mientras lo rompo en datos no balanceados, pero la idea es ...

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

Finalmente, el ya mencionado libro de Piniero y Bates (2000) favoreció fuertemente el enrejado de lo poco que he hojeado. Tal vez algo como trazar los datos en bruto ...

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

Y luego trazando los valores ajustados ...

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

Answer 3

Usando coefplot2 (en r-forge):

Robando el 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 de ajuste:

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

Trama:

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

editar :

Con frecuencia he tenido problemas con la versión R-Forge. Este respaldo debería funcionar si la compilación de R-Forge no funciona:

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

Tenga en cuenta que la dependencia de coda ya debe estar instalada.