Come per tracciare i risultati di un modello misto

Question

che uso lme4 in R per adattarsi al modello misto

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

dove il valore è continua, lo stato (N/D/R) e l'esperimento sono fattori, e ottengo

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 

Mi piacerebbe rappresentare graficamente la valutazione degli effetti fissi. Tuttavia, sembra che non ci sia alcuna funzione grafica per questi oggetti. C'è un modo in cui posso rappresentare graficamente gli effetti fissi?

Answer 1

Ecco alcuni suggerimenti.

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

Utilizzando coefplot2 (su R-forge):

Rubare il codice di simulazione da @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 

Modello adatto:

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) 

modificare:

Ho spesso avendo problemi con l'accumulo R-Forge. Questo fallback dovrebbe funzionare se la generazione R-Forge non funziona:

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

Nota che la dipendenza coda deve essere già installato.

Answer 3

mi piace la fiducia coefficiente trame di intervallo, ma può essere utile considerare alcuni appezzamenti supplementari per capire gli effetti fissi ..

Rubare il codice di simulazione da @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) 

Ottenere un guardare la struttura dei dati ... sembra equilibrata ..

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

Potrebbe essere interessante tracciare la correlazione tra effetti fissi, soprattutto se si adattano strutture di correlazione diverse. C'è un po 'di codice fresco fornito al seguente link ...

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

Infine sembra rilevante per guardare la variabilità tra i 10 esperimenti, nonché la variabilità tra "Stato "all'interno di esperimenti. Sto ancora lavorando al codice per questo perché lo spezzo su dati sbilanciati, ma l'idea è ...

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

Infine, il già citato Piniero e Bates (2000) libro reticolo fortemente favorita da quel poco che ho scremato .. Così si potrebbe dare che un colpo. Forse qualcosa di simile tracciare i dati grezzi ...

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

E poi riportando i valori adattati ...

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