मिश्रित प्रभाव मॉडल से केवल "महत्वपूर्ण" यादृच्छिक प्रभावों का एक कैटरपिलर साजिश

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 * से और 1.96 * से नहीं)।

मूल रूप से, मैं एक कैटरपिलर साजिश है कि के लिए 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) जैसे कसौटी | जेड |> 1.5 या | जेड |> 1.75 के बजाय, बस पर जोर देना है कि इस नहीं एक आनुमानिक सीमा है ...)

मैं एक छोटे से हो रही समाप्त हो गया दूर ले जाया गया ... मैंने फैसला किया कि चीजों को थोड़ा सा रेफैक्टर/मॉड्यूलर करना बेहतर होगा, इसलिए मैंने augment विधि लिखा (के साथ काम करने के लिए डिज़ाइन किया गयापैकेज) जो ranef.mer ऑब्जेक्ट्स से उपयोगी डेटा फ्रेम बनाता है ... एक बार ऐसा करने के बाद, आप जो मैनिप्ल्यूशन चाहते हैं वह बहुत आसान है।

मैंने अपने उत्तर के अंत में augment.ranef.mer कोड डाला - यह थोड़ा लंबा है (आपको यहां कोड चलाने से पहले इसे स्रोत करना होगा)।

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

आरई वस्तु को augment विधि लागू करें:

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() के बजाय कर रहा हूं क्योंकि ggplot2coord_flipfacet_wrap(...,scales="free") के साथ उपयोग नहीं कर सकता ...के साथ "महत्वपूर्ण" ढलानों या अवरोध केवल स्तर के साथ

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

अद्यतन कैटरपिलर साजिश:

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

अब स्तरों आप रखना चाहते हैं खोजने के

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