Gąsienicowa fabuła tylko "znaczących" efektów losowych z mieszanego modelu efektów

Question

Miałem świetne doświadczenia, prosząc o pomoc tutaj i mam nadzieję, że znowu pomogę.

Oceniam dość duży mieszany model efektów, w którym jeden z efektów losowych ma ponad 150 różnych poziomów. To spowodowałoby, że standardowy spisek gąsienicowy byłby zupełnie nieczytelny.

Chciałbym, o ile to możliwe, uzyskać wykres gąsienicy o samych poziomach losowego efektu, które z braku lepszego określenia są "znaczące". To znaczy: Chcę wykresu gąsienicy, w którym losowe przechwycenie o losowym nachyleniu dla zmiennego współczynnika ma "przedział ufności" (wiem, że to nie jest to, co to jest), który nie zawiera zera.

Należy wziąć pod uwagę ten standardowy model z danych sleepstudy, który jest standardem z lme4.

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

Zdobędę tę historię gąsienicy.

Działka gąsienica używam pochodzi z this code. Zauważ, że używam mniej konserwatywnych granic dla interwałów (tj. 1.645 * se, a nie 1.96 * se).

Zasadniczo chcę działkę gąsienica, która właśnie zawierać poziomy dla 308, 309, 310, 330, 331, 335, 337, 349, 350, 352 i 370, ponieważ te poziomy były albo przechwytuje lub stoki w których interwałach nie było zera. Pytam, ponieważ moja gąsienica o ponad 150 różnych poziomach jest nieczytelna i myślę, że to może być warte tego rozwiązanie.

Powtarzalny kod jest następujący. Naprawdę doceniam każdą pomoc.

# 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

Po pierwsze, dzięki za umieszczenie „znaczące” w cudzysłowie ... każdy czyta to powinien pamiętać, że znaczenie ma nie mieć żadnego statystycznych znaczenie w tym kontekście (może lepiej byłoby użyć Z- statystyka (wartość/std.error) kryterium takie jak | Z |> 1.5 lub | Z |> 1.75 zamiast tego, aby podkreślić, że jest to nie próg inferencyjny ...)

Skończyło się na tym, że otrzymałem trochę poniechany ... Postanowiłem, że lepiej będzie refaktoryzować/modularyzować rzeczy, więc napisałem metodę augment (zaprojektowaną do pracy zPakiet), który tworzy użyteczne ramki danych z obiektów ranef.mer ... po wykonaniu tej operacji pożądane manipulacje są dość łatwe.

Na końcu mojej odpowiedzi umieszczam kod augment.ranef.mer - jest on trochę długi (musisz go pobrać, zanim będziesz mógł uruchomić kod tutaj).

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

Zastosować metodę augment do obiektu RE:

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

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

Teraz kod ggplot jest dość podstawowe. Używam geom_errorbarh(height=0) zamiast geom_pointrange()+coord_flip(), ponieważ ggplot2 nie może używać coord_flip z 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") 

teraz znaleźć poziomy, które chcesz zachować:

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

Aktualizacja gąsienica działkę tylko poziomach z „znaczące” zboczach lub przechwytuje:

g1 %+% aa3 

Gdybyś tylko chciał, aby podświetlić "istotne" poziomy zamiast usuwania "nieistotnych" poziomów całkowicie

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