Гусеничный сюжет только «значительных» случайных эффектов из модели смешанных эффектов

Question

У меня были большие впечатления, прося о помощи здесь раньше, и я надеюсь снова получить помощь.

Я оцениваю довольно большую модель смешанных эффектов, в которой один из случайных эффектов имеет более 150 различных уровней. Это сделало бы стандартный гусеничный сюжет довольно нечитаемым.

Я хотел бы, если вообще возможно, получить участок гусеницы только уровней случайного эффекта, которые являются для отсутствия лучшего термина «значительными». То есть: мне нужен гусеничный график, в котором либо случайный перехват , либо случайный уклон для изменяющегося коэффициента имеет «доверительный интервал» (я знаю, что это не совсем то, что это), который не включает нуль.

Рассмотрите эту стандартную модель из данных sleepstudy, которая является стандартной с lme4.

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

Я бы выбрал этот участок гусеницы.

Гусеница участок я использую приходит от this code. Заметьте, я склонен использовать менее консервативные оценки для интервалов (то есть 1,645 * se, а не 1,96 * se).

В принципе, я хочу гусеница сюжет, который будет содержать только уровни для 308, 309, 310, 330, 331, 335, 337, 349, 350, 352, и 370, так как эти уровни были либо перехватывает или склоны интервалы которых не равны нулю. Я спрашиваю, потому что мой гусеничный сюжет более 150 разных уровней нечитабельен, и я думаю, что это может быть полезным решением.

Воспроизводимый код следует. Я искренне ценю любую помощь.

# 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

Во-первых, спасибо за ввод «значительный» в кавычки ... каждый читающий это должен помнить, что значение не имеет статистический значение в этом контексте (это может быть лучше использовать Z- статистика (значение/std.error) критерий, такой как | Z |> 1,5 или | Z |> 1,75 вместо этого, просто чтобы подчеркнуть, что это не выведенный порог ...)

я в конечном итоге получить немного увлекся ... Я решил, что было бы лучше немного рефакторировать/модулировать вещи, поэтому я написал метод augment (предназначен для работы с broom), который создает полезные кадры данных из объектов ranef.mer ... как только это будет сделано, требуемые манипуляции довольно просты.

Я положил код augment.ranef.mer в конце своего ответа - он немного длинный (вам нужно будет указать его, прежде чем вы сможете запустить код здесь).

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

Применить метод augment к объекту RE:

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

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

Теперь ggplot код довольно простой. Я использую geom_errorbarh(height=0), а не geom_pointrange()+coord_flip(), потому что ggplot2 не может использовать coord_flip с 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") 

Теперь найти уровни, которые вы хотите сохранить:

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

Update гусеничный участок только с уровнями с «значительными» откосах или перехватывает:

g1 %+% aa3 

Если вы только хотели подчеркнуть «значительных» уровней, а не полностью исключить «несущественные» уровни

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