In R, tracciando gli effetti casuali da lmer (pacchetto lme4) usando qqmath o dotplot: come farlo sembrare di fantasia?

Question

La funzione qqmath crea grandi diagrammi a bruco di effetti casuali usando l'output del pacchetto lmer. Cioè, qqmath è bravissimo nel tracciare le intercettazioni da un modello gerarchico con i loro errori attorno alla stima puntuale. Un esempio delle funzioni lmer e qqmath sono sotto usando i dati incorporati nel pacchetto lme4 chiamato Dyestuff. Il codice produrrà il modello gerarchico e una trama piacevole usando la funzione ggmath.

library("lme4") 
data(package = "lme4") 

# Dyestuff 
# a balanced one-way classiï¬cation of Yield 
# from samples produced from six Batches 

summary(Dyestuff)    

# Batch is an example of a random effect 
# Fit 1-way random effects linear model 
fit1 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff) 
summary(fit1) 
coef(fit1) #intercept for each level in Batch 

# qqplot of the random effects with their variances 
qqmath(ranef(fit1, postVar = TRUE), strip = FALSE)$Batch 

L'ultima riga di codice produce una bella trama di ogni intercettazione con l'errore attorno ad ogni stima. Ma la formattazione della funzione qqmath sembra essere molto difficile, e ho faticato a formattare la trama. Mi è venuta in mente alcune domande che non posso rispondere, e che penso che altri potrebbero beneficiare anche se stanno usando la combinazione lmer/qqmath:

1. C'è un modo per prendere la funzione qqmath sopra e aggiungi alcune opzioni , ad esempio, rendendo determinati punti vuoti rispetto a quelli riempiti o colori diversi per punti diversi? Ad esempio, puoi fare il pieno dei punti per A, B e C della variabile Batch, ma poi il resto dei punti è vuoto?
2. È possibile aggiungere etichette di assi per ciascun punto (ad esempio, sull'asse superiore o destro, ad esempio)?
3. I miei dati sono più vicini a 45 intercetta, quindi è possibile aggiungere la spaziatura tra le etichette in modo che non si incontrino l'una con l'altra? PRINCIPALMENTE, sono interessato a distinguere/etichettare tra i punti sul grafico , che sembra essere ingombrante/impossibile nella funzione ggmath.

Finora, l'aggiunta di qualsiasi opzione aggiuntiva negli errori producono funzione qqmath dove non otterrebbe gli errori se si trattasse di un complotto di serie, quindi sono in perdita.

ANCHE, se ritieni che ci sia un pacchetto/funzione migliore per tracciare le intercettazioni dall'output di lmer, mi piacerebbe sentirlo! (per esempio, puoi fare punti 1-3 usando dotplot?)

Grazie.

EDIT: Sono anche aperto a un dotplot alternativo se può essere formattato in modo ragionevole. Mi piace l'aspetto di una trama ggmath, quindi sto iniziando con una domanda al riguardo.

Answer 1

Una possibilità è utilizzare la libreria ggplot2 per disegnare un grafico simile e quindi è possibile regolare l'aspetto del grafico.

In primo luogo, ranef oggetto viene salvato come randoms. Quindi le variazioni delle intercettazioni vengono salvate nell'oggetto qq.

randoms<-ranef(fit1, postVar = TRUE) 
qq <- attr(ranef(fit1, postVar = TRUE)[[1]], "postVar") 

oggetto rand.interc contiene solo intercettazioni casuali con nomi di livello.

rand.interc<-randoms$Batch 

Tutti gli oggetti messi in un frame di dati. Per gli intervalli di errore, sd.interc viene calcolato come 2 volte radice quadrata della varianza.

df<-data.frame(Intercepts=randoms$Batch[,1], 
       sd.interc=2*sqrt(qq[,,1:length(qq)]), 
       lev.names=rownames(rand.interc)) 

Se è necessario che le intercettazioni sono ordinate nella trama in base al valore di allora lev.names dovrebbe essere riordinati.Questa linea può essere saltata se le intercettazioni devono essere ordinate per nome di livello.

df$lev.names<-factor(df$lev.names,levels=df$lev.names[order(df$Intercepts)]) 

Questo codice produce trama. Ora i punti saranno diversi per shape in base ai livelli di fattore. risposta

library(ggplot2) 
p <- ggplot(df,aes(lev.names,Intercepts,shape=lev.names)) 

#Added horizontal line at y=0, error bars to points and points with size two 
p <- p + geom_hline(yintercept=0) +geom_errorbar(aes(ymin=Intercepts-sd.interc, ymax=Intercepts+sd.interc), width=0,color="black") + geom_point(aes(size=2)) 

#Removed legends and with scale_shape_manual point shapes set to 1 and 16 
p <- p + guides(size=FALSE,shape=FALSE) + scale_shape_manual(values=c(1,1,1,16,16,16)) 

#Changed appearance of plot (black and white theme) and x and y axis labels 
p <- p + theme_bw() + xlab("Levels") + ylab("") 

#Final adjustments of plot 
p <- p + theme(axis.text.x=element_text(size=rel(1.2)), 
       axis.title.x=element_text(size=rel(1.3)), 
       axis.text.y=element_text(size=rel(1.2)), 
       panel.grid.minor=element_blank(), 
       panel.grid.major.x=element_blank()) 

#To put levels on y axis you just need to use coord_flip() 
p <- p+ coord_flip() 
print(p) 

Answer 2

Didzis' è grande! Giusto per riassumerlo un po ', lo metto in una sua funzione che si comporta molto come qqmath.ranef.mer() e dotplot.ranef.mer(). Oltre alla risposta di Didzis, gestisce anche modelli con più effetti casuali correlati (come qqmath() e dotplot() do). Confronto con qqmath():

require(lme4)       ## for lmer(), sleepstudy 
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 
ggCaterpillar(ranef(fit, postVar=TRUE)) ## using ggplot2 
qqmath(ranef(fit, postVar=TRUE))   ## for comparison 

Confronto con dotplot():

ggCaterpillar(ranef(fit, postVar=TRUE), QQ=FALSE) 
dotplot(ranef(fit, postVar=TRUE)) 

A volte, potrebbe essere utile avere diverse scale per gli effetti casuali - cosa che dotplot() fa rispettare. Quando ho cercato di rilassarmi, ho dovuto cambiare la sfaccettatura (vedi questo answer).

ggCaterpillar(ranef(fit, postVar=TRUE), QQ=FALSE, likeDotplot=FALSE) 

## re = object of class ranef.mer 
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE) { 
    require(ggplot2) 
    f <- function(x) { 
     pv <- attr(x, "postVar") 
     cols <- 1:(dim(pv)[1]) 
     se <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ]))) 
     ord <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x)) 
     pDf <- data.frame(y=unlist(x)[ord], 
          ci=1.96*se[ord], 
          nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)), 
          ID=factor(rep(rownames(x), ncol(x))[ord], levels=rownames(x)[ord]), 
          ind=gl(ncol(x), nrow(x), labels=names(x))) 

     if(QQ) { ## normal QQ-plot 
      p <- ggplot(pDf, aes(nQQ, y)) 
      p <- p + facet_wrap(~ ind, scales="free") 
      p <- p + xlab("Standard normal quantiles") + ylab("Random effect quantiles") 
     } else { ## caterpillar dotplot 
      p <- ggplot(pDf, aes(ID, y)) + coord_flip() 
      if(likeDotplot) { ## imitate dotplot() -> same scales for random effects 
       p <- p + facet_wrap(~ ind) 
      } else {   ## different scales for random effects 
       p <- p + facet_grid(ind ~ ., scales="free_y") 
      } 
      p <- p + xlab("Levels") + ylab("Random effects") 
     } 

     p <- p + theme(legend.position="none") 
     p <- p + geom_hline(yintercept=0) 
     p <- p + geom_errorbar(aes(ymin=y-ci, ymax=y+ci), width=0, colour="black") 
     p <- p + geom_point(aes(size=1.2), colour="blue") 
     return(p) 
    } 

    lapply(re, f) 
} 

Answer 3

Un altro modo per farlo è quello di estrarre i valori simulati dalla distribuzione di ciascuno degli effetti casuali e tracciare quelle. Utilizzando il pacchetto merTools, è possibile ottenere facilmente le simulazioni da un oggetto lmer o glmer e tracciarle.

library(lme4); library(merTools)  ## for lmer(), sleepstudy 
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 
randoms <- REsim(fit, n.sims = 500) 

randoms è ora un oggetto con che assomiglia:

head(randoms) 
groupFctr groupID  term  mean  median  sd 
1 Subject  308 (Intercept) 3.083375 2.214805 14.79050 
2 Subject  309 (Intercept) -39.382557 -38.607697 12.68987 
3 Subject  310 (Intercept) -37.314979 -38.107747 12.53729 
4 Subject  330 (Intercept) 22.234687 21.048882 11.51082 
5 Subject  331 (Intercept) 21.418040 21.122913 13.17926 
6 Subject  332 (Intercept) 11.371621 12.238580 12.65172 

Esso fornisce il nome del fattore di raggruppamento, il livello del fattore stiamo ottenendo una stima per il termine nel modello e la media, la mediana e la deviazione standard dei valori simulati. Possiamo usare questo per generare una trama bruco simili a quelli di cui sopra:

plotREsim(randoms) 

che produce:

Una caratteristica interessante è che i valori che hanno un intervallo di confidenza che non si sovrapponga a zero sono evidenziati in nero. È possibile modificare la larghezza dell'intervallo utilizzando il parametro level su plotREsim rendendo intervalli di confidenza più ampi o più ristretti in base alle proprie esigenze.