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.