Ein Raupe-Plot von nur den "signifikanten" Zufallseffekten aus einem Mixed-Effects-Modell

Question

Ich hatte großartige Erfahrungen, um Hilfe hier zu fragen, und ich hoffe, wieder Hilfe zu bekommen.

Ich schätze ein ziemlich großes Mixed-Effekt-Modell, in dem einer der zufälligen Effekte über 150 verschiedene Ebenen hat. Das würde ein Standard-Raupe-Plot ziemlich unlesbar machen.

Ich möchte, wenn alles möglich ist, eine Raupe Handlung von nur den Ebenen des zufälligen Effekts, die für einen Mangel an besseren Begriff "signifikant" sind. Das heißt: Ich möchte ein Raupe-Diagramm, in dem entweder der zufällige Schnittpunkt oder die zufällige Steigung für einen variierenden Koeffizienten ein "Konfidenzintervall" hat (ich weiß, dass es nicht ganz das ist), das null nicht enthält.

Betrachten Sie dieses Standardmodell aus den sleepstudy Daten, die standardmäßig mit lme4 sind.

library(lme4) 
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 
ggCaterpillar(ranef(fit,condVar=TRUE), QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]] 

Ich würde diese Raupe Plot bekommen.

Die Raupe Grundstück Ich benutze kommt aus this code. Beachten Sie, dass ich dazu tendiere, weniger konservative Grenzen für die Intervalle zu verwenden (d. H. 1,645 * se und nicht 1,96 * se).

Grundsätzlich mag ich eine Raupe Handlung, die die Ebene nur einschließen würde für 308, 309, 310, 330, 331, 335, 337, 349, 350, 352 und 370, da diese Werte entweder abfängt hatten oder Pisten wessen Intervalle nicht Null enthalten. Ich frage, weil mein Raupe-Plot von über 150 verschiedenen Ebenen irgendwie unlesbar ist und ich denke, dass dies eine sinnvolle Lösung sein könnte.

Nachführbarer Code folgt. Ich schätze wirklich jede Hilfe.

# https://stackoverflow.com/questions/34120578/how-can-i-sort-random-effects-by-value-of-the-random-effect-not-the-intercept 
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE, reorder=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, ]))) 
if (reorder) { 
    ord <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x)) 
    pDf <- data.frame(y=unlist(x)[ord], 
        ci=1.645*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))) 
} else { 
    pDf <- data.frame(y=unlist(x), 
        ci=1.645*se, 
        nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)), 
        ID=factor(rep(rownames(x), ncol(x)), levels=rownames(x)), 
        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 of the Random Effect") + ylab("Random Effect") 
} 

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


library(lme4) 
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 
ggCaterpillar(ranef(fit,condVar=TRUE), QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]] 
ggsave(file="sleepstudy.png") 

Answer 1

Zunächst vielen Dank für „signifikant“ in Anführungszeichen setzen ... jeder der Lektüre dieses, dass die Bedeutung erinnern sollte keine statistische Bedeutung in diesem Zusammenhang (es könnte besser sein, einen Z- zu verwenden Statistik (Wert/std.error) Kriterium wie | Z |> 1,5 oder | Z |> 1,75 stattdessen nur zu betonen, dass dies nicht eine folger Schwelle ...)

landete ich ein wenig bekommen weggetragen ... Ich entschied, dass es besser wäre, Dinge ein bisschen zu refaktorisieren/zu modulieren, also schrieb ich eine augment Methode (entworfen, um mit derzu arbeitenPaket), die nützliche Datenrahmen von ranef.mer Objekte erstellt ... sobald dies erledigt ist, sind die Manipulationen, die Sie wollen, ziemlich einfach.

Ich habe den Code augment.ranef.mer am Ende meiner Antwort - es ist ein bisschen lang (Sie müssen es Quellcode, bevor Sie den Code hier ausführen können).

library(broom) 
library(reshape2) 
library(plyr) 

Übernehmen Sie die augment Methode auf die RE-Objekt:

rr <- ranef(fit,condVar=TRUE) 
aa <- augment(rr) 

names(aa) 
## [1] "grp"  "variable" "level"  "estimate" "qq"  "std.error" 
## [7] "p"   "lb"  "ub"  

Nun ist die ggplot Code ist ziemlich einfach. Ich benutze geom_errorbarh(height=0) anstatt geom_pointrange()+coord_flip() weil ggplot2coord_flip mit facet_wrap(...,scales="free") nicht verwenden kann ...

## Q-Q plot: 
g0 <- ggplot(aa,aes(estimate,qq,xmin=lb,xmax=ub))+ 
    geom_errorbarh(height=0)+ 
    geom_point()+facet_wrap(~variable,scale="free_x") 

## regular caterpillar plot: 
g1 <- ggplot(aa,aes(estimate,level,xmin=lb,xmax=ub))+ 
    geom_errorbarh(height=0)+ 
    geom_vline(xintercept=0,lty=2)+ 
    geom_point()+facet_wrap(~variable,scale="free_x") 

nun die Ebenen, die Sie behalten möchten, finden:

aa2 <- ddply(aa,c("grp","level"), 
      transform, 
      keep=any(p<0.05)) 
aa3 <- subset(aa2,keep) 

aktualisieren Raupe Grundstück mit nur Ebenen mit „signifikant“ Hängen oder abfängt:

g1 %+% aa3 

Wenn Sie nur markieren wollte "signifikante" Stufen, anstatt "nicht signifikante" Stufen vollständig zu entfernen

ggplot(aa2,aes(estimate,level,xmin=lb,xmax=ub,colour=factor(keep)))+ 
    geom_errorbarh(height=0)+ 
    geom_vline(xintercept=0,lty=2)+ 
    geom_point()+facet_wrap(~variable,scale="free_x")+ 
    scale_colour_manual(values=c("black","red"),guide=FALSE) 

---

##' @importFrom reshape2 melt 
##' @importFrom plyr ldply name_rows 
augment.ranef.mer <- function(x, 
           ci.level=0.9, 
           reorder=TRUE, 
           order.var=1) { 
    tmpf <- function(z) { 
     if (is.character(order.var) && !order.var %in% names(z)) { 
      order.var <- 1 
      warning("order.var not found, resetting to 1") 
     } 
     ## would use plyr::name_rows, but want levels first 
     zz <- data.frame(level=rownames(z),z,check.names=FALSE) 
     if (reorder) { 
      ## if numeric order var, add 1 to account for level column 
      ov <- if (is.numeric(order.var)) order.var+1 else order.var 
      zz$level <- reorder(zz$level, zz[,order.var+1], FUN=identity) 
     } 
     ## Q-Q values, for each column separately 
     qq <- c(apply(z,2,function(y) { 
        qnorm(ppoints(nrow(z)))[order(order(y))] 
       })) 
     rownames(zz) <- NULL 
     pv <- attr(z, "postVar") 
     cols <- 1:(dim(pv)[1]) 
     se <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ]))) 
     ## n.b.: depends on explicit column-major ordering of se/melt 
     zzz <- cbind(melt(zz,id.vars="level",value.name="estimate"), 
        qq=qq,std.error=se) 
     ## reorder columns: 
     subset(zzz,select=c(variable, level, estimate, qq, std.error)) 
    } 
    dd <- ldply(x,tmpf,.id="grp") 
    ci.val <- -qnorm((1-ci.level)/2) 
    transform(dd, 
       p=2*pnorm(-abs(estimate/std.error)), ## 2-tailed p-val 
       lb=estimate-ci.val*std.error, 
       ub=estimate+ci.val*std.error) 
}