W R, wykreślając losowe efekty z lmer (pakiet lme4) za pomocą qqmath lub dotplot: jak sprawić, by wyglądał fantazyjnie?

Question

Funkcja qqmath tworzy świetne wykresy losowe przy użyciu wyjścia z pakietu lmer. Oznacza to, że qqmath jest świetny w wykreślaniu przechwyceń z modelu hierarchicznego z ich błędami wokół oszacowania punktu. Przykład funkcji lmer i qqmath znajduje się poniżej przy użyciu wbudowanych danych w pakiecie lme4 o nazwie Dyestuff. Kod wygeneruje model hierarchiczny i ładną fabułę za pomocą funkcji 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 

Ostatnia linia kodu tworzy naprawdę ładny wykres każdego punktu przecięcia z błędem wokół każdego oszacowania. Ale formatowanie funkcji qqmath wydaje się być bardzo trudne i starałem się sformatować fabułę. Mam wymyślić kilka pytań, na które nie mogę odpowiedzieć, a myślę, że inni mogą również skorzystać z, jeśli są one za pomocą kombinacji lmer/qqmath:

1. Czy istnieje sposób, aby wziąć powyżej funkcję qqmath i dodać kilka opcji, takich jak: czy niektóre punkty są puste, czy wypełnione, czy różne kolory dla różnych punktów? Na przykład, czy możesz ustawić punkty dla A, B i C zmiennej Batch, ale pozostałe punkty są puste?
2. Czy można dodać etykiety osi dla każdego punktu (np. Wzdłuż górnej lub prawej osi )?
3. Moje dane mają blisko 45 przechwyceń, więc istnieje możliwość dodania odstępów między etykietami, aby nie do siebie pasowały? GŁÓWNIE, jestem zainteresowany rozróżnianiem/etykietowaniem punktów na wykresie , który wydaje się być niewygodny/niemożliwy w funkcji ggmath.

Do tej pory dodanie dowolnej dodatkowej opcji w funkcji qqmath powoduje błędy, które nie będą powodowały błędów, jeśli jest to standardowy wykres, więc jestem w błędzie.

TAKŻE, jeśli uważasz, że jest lepszy pakiet/funkcja do wykreślania przechwyceń z lm wyjściowego, chciałbym to usłyszeć! (na przykład, czy możesz zrobić punkty 1-3 używając dotplot?)

Dzięki.

EDYCJA: Jestem również otwarty na alternatywny dotplot, jeśli można go w rozsądny sposób sformatować. Po prostu podoba mi się wygląd fabuły ggmath, więc zaczynam od pytania na ten temat.

Answer 1

Jedną z możliwości jest użycie biblioteki ggplot2 do narysowania podobnego wykresu, a następnie można dostosować wygląd wydruku.

Po pierwsze, obiekt ranef zostaje zapisany jako randoms. Następnie warianty przechwyceń są zapisywane w obiekcie qq.

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

Obiekt rand.interc zawiera tylko przypadkowe przechwycenia z nazwami poziomów.

rand.interc<-randoms$Batch 

Wszystkie obiekty umieszczane w jednej ramce danych. Dla przedziałów błędów sd.interc jest obliczany jako 2-krotność pierwiastka kwadratowego wariancji.

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

Jeśli trzeba, że ​​przechwytuje są uporządkowane w zależności od wartości działki następnie lev.names powinna być kolejność.Linię tę można pominąć, jeśli przechwyty powinny być uporządkowane według nazw poziomów.

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

Ten kod powoduje utworzenie wykresu. Teraz punkty będą się różnić o shape według poziomów czynników. Odpowiedź

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' jest super! Aby go trochę zawinąć, umieściłem go w swojej własnej funkcji, która zachowuje się podobnie jak qqmath.ranef.mer() i dotplot.ranef.mer(). Oprócz odpowiedzi Didzisa, obsługuje ona również modele z wieloma skorelowanymi efektami losowymi (takimi jak qqmath() i dotplot()). Porównanie do 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 

porównanie dotplot():

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

Czasami może to być przydatne mieć różne wagi dla efektów losowych - coś, co dotplot() wymusza. Kiedy próbowałem to rozluźnić, musiałem zmienić fasetowanie (patrz: 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

Innym sposobem, aby to zrobić, aby wyodrębnić wartości symulowanych z podziału każdej z efektami losowymi i wykreślić tych. Korzystając z pakietu merTools, można łatwo uzyskać symulacje z obiektu lmer lub glmer i narysować je.

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

randoms jest teraz obiekt o który wygląda następująco:

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 

Zapewnia nazwę czynnika grupującego, poziom czynnika jesteśmy uzyskanie szacunkowych na termin w modelu oraz średnią, medianę i odchylenie standardowe symulowanych wartości. Możemy to wykorzystać, aby wygenerować wykres gąsienica podobne do tych powyżej:

plotREsim(randoms) 

która produkuje:

Jedną miłą cechą jest to, że wartości, które mają przedział ufności, które nie pokrywają się zera są podświetlone na czarno. Można zmienić szerokość przedziału, używając parametru level, aby uzyskać szersze lub węższe przedziały ufności w zależności od potrzeb.