Una trama bruciante di soli effetti "significativi" casuali da un modello a effetti misti

Question

Ho avuto esperienze fantastiche per chiedere aiuto qui prima e spero di ricevere ancora aiuto.

Sto stimando un modello di effetti misti piuttosto ampio in cui uno degli effetti casuali ha oltre 150 livelli diversi. Ciò renderebbe la trama di un bruco standard piuttosto illeggibile.

Vorrei, se possibile, ottenere una trama di bruco di solo i livelli dell'effetto casuale che sono, per mancanza di termini migliori, "significativi". Cioè: voglio una trama di bruco in cui sia l'intercetta casuale o la pendenza casuale per un coefficiente variabile abbia un "intervallo di confidenza" (so che non è proprio quello che è) che non include zero.

Considerare questo modello standard dai dati sleepstudy standard con lme4.

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

vorrei ottenere questa trama bruco.

La trama bruco che uso viene da this code. Si noti che io tendo ad usare limiti meno conservativi per gli intervalli (cioè 1.645 * se e non 1.96 * se).

Fondamentalmente, io voglio una trama bruco che sarebbe solo includere i livelli per 308, 309, 310, 330, 331, 335, 337, 349, 350, 352, e 370, perché quei livelli avevano entrambi intercetta o pendenze i cui intervalli non includevano zero. Lo chiedo perché la mia trama di bruco di oltre 150 livelli diversi è un po 'illeggibile e penso che potrebbe essere una soluzione utile.

Segue il codice riproducibile. Apprezzo sinceramente qualsiasi aiuto.

# 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

In primo luogo, grazie per mettere "significativa" tra virgolette ... tutti la lettura di questo dovrebbe ricordare che il significato non ha statistica significato in questo contesto (che potrebbe essere meglio utilizzare un Z- statistica (valore/std.error) criterio come | z |> 1.5 o | Z |> 1,75 invece, proprio per sottolineare che si tratta non una soglia inferenziale ...)

ho finito per ottenere un po ' portato via ... Ho deciso che sarebbe stato meglio rifattare/modulare un po 'le cose, quindi ho scritto un metodo augment (progettato per funzionare conPacchetto) che costruisce frame di dati utili dagli oggetti ranef.mer ... una volta eseguita questa operazione, le manipolazioni desiderate sono piuttosto semplici.

Inserisco il codice augment.ranef.mer alla fine della mia risposta: è un po 'lungo (è necessario procurarselo prima di poter eseguire il codice qui).

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

applicare il metodo augment all'oggetto RE:

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

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

Ora il codice ggplot è piuttosto semplice. Sto usando geom_errorbarh(height=0) anziché geom_pointrange()+coord_flip() perché ggplot2 non è possibile utilizzare coord_flip con facet_wrap(...,scales="free") ...

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

Ora trovare i livelli che si desidera conservare:

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

Aggiornamento trama bruco con solo i livelli con pendenze "significative" o intercettazioni:

g1 %+% aa3 

Se si voleva solo mettere in evidenza livelli "significativi" piuttosto che rimuovere completamente livelli "non significativi"

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