首页 文章

lme4计算协方差的置信区间

提问于
浏览 384
5

请参阅Ben Bolker 16/05/2016的答案,了解相应的解决方案 . OP下面 .


我正在使用lme4安装几个多级模型 . 我想报告随机效应的方差和协方差,并自动化这个过程 .

我知道我可以得到 as.data.frame(VarCorr(mymodel)) 的差异,我知道我可以用 confint(mymodel) 得到置信区间 . 很明显,我可以合并/绑定两个表,并通过简单地将 confint() 的输出平方在适当的行和列上,将置信区间放在方差周围,但是如果不是手动的话,我无法找到一个令人信服的方法来计算协方差 . .

confint 的结果是:

conf <- NULL
a <- c(6.2,-0.4,2.2,1.5,-0.4,-0.5,2.8,-0.9,1.3,3.9)
b <- c(6.8,-0.2,2.5,2.5,0.1,0.2,4.8,-0.7,2.3,5)
conf <- data.frame(a,b,row.names = c("sd_(Intercept)|ID","cor_Time.(Intercept)|ID","sd_Time|ID","sd_(Intercept)|Group","cor_Time.(Intercept)|Group","cor_I(Time^2).(Intercept)|Group","sd_Time|Group","cor_I(Time^2).Time|Group","sd_I(Time^2)|Group","sigma"))
colnames(conf) <- c("2.5%","97.5%")
conf

如何自动化各种乘法以获得协方差,如

cov.time.intercept <- conf[1,2]*conf[1,1]*conf[1,3]

我试过分割标准偏差和相关性,创建“ID”,“时间”,“我(时间^ 2)”和“(拦截)”变量然后匹配两列,但我没有得到任何结果 . 问题是,每次模型更改时,您可能会有不同数量的方差和协方差,以及不同的三角矩阵 .

感谢您的任何帮助,

ķ .

3 回答

  • 2

    请注意, lme4 摘要中随机效应的标准差是 NOT 方差的标准误差!这只是方差的平方根!

    如果你需要对随机效应方差的置信区间,那么你需要 profile() 的可能性 . 见 ?lme4::profile .

  • 1

    你的计算似乎给出了合理的答案,但它(对我来说;我随时准备纠正/开悟......) . 假设 cov = corr*var1*var2 . 假设 ci(.) 是数量的(下限或上限)置信限 . ci(cov) = ci(corr)*ci(var1)*ci(var2) (有趣的是你获得了合理的答案;我认为当数量大致不相关时,这很有可能发生......)如果你有每个成分的方差和它们之间的协方差(我并不是指随机效应方差和协方差本身,而是它们的采样方差/协方差)你可以使用delta方法近似传播它们,但这些很难获得(参见here) .

    据我所知,做到这一点的“正确”方法是在方差 - 协方差量表上进行似然概率计算而不是标准偏差 - 相关量表 . 这在以前是不可能的,但它现在(使用Github上的开发版本) .

    安装最新版本:

    library(remotes) ## for install_github (or library(devtools))
    install_github("lme4/lme4")
    

    预赛:

    chap12 <- foreign::read.dta(file = "ch12.dta")
    library(lme4)
    snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher),
                     data = chap12)
    
    as.data.frame(VarCorr(snijders))
    ##        grp        var1 var2        vcov      sdcor
    ## 1  teacher (Intercept) <NA>  0.15617962  0.3951957
    ## 2  teacher         occ <NA>  0.01205317  0.1097869
    ## 3  teacher (Intercept)  occ -0.03883458 -0.8950676
    ## 4 Residual        <NA> <NA>  0.04979762  0.2231538
    

    在比较结果时我们必须要小心,因为 profile.merMod ,我们看起来像这样会产生巨大的差异 .

    s2 <- refitML(snijders)
    as.data.frame(VarCorr(s2))
    ##        grp        var1 var2        vcov      sdcor
    ## 1  teacher (Intercept) <NA>  0.15426049  0.3927601
    ## 2  teacher         occ <NA>  0.01202631  0.1096645
    ## 3  teacher (Intercept)  occ -0.03884427 -0.9018483
    ## 4 Residual        <NA> <NA>  0.04955549  0.2226106
    
    p.sd <- profile(s2,which="theta_",
                  signames=FALSE)
    p.vcov <- profile(s2,which="theta_",prof.scale="varcov",
                  signames=FALSE)
    

    我们收到一些关于非单调轮廓的警告......

    confint(p.vcov)
    ##                                    2.5 %     97.5 %
    ## var_(Intercept)|teacher      0.08888931  0.26131067
    ## cov_occ.(Intercept)|teacher -0.07553263 -0.01589043
    ## var_occ|teacher              0.00000000  0.02783863
    ## sigma                        0.03463184  0.07258777
    

    如果我们检查相关(sd / variance)元素的平方怎么办?

    confint(p.sd)[c(1,3,4),]^2
    ##                              2.5 %     97.5 %
    ## sd_(Intercept)|teacher 0.089089363 0.26130970
    ## sd_occ|teacher         0.002467408 0.02779329
    ## sigma                  0.034631759 0.07263869
    

    这些匹配很好,除了 occ 方差的下限;它们也符合您上面的结果 . 但是,协方差结果(我声称这是很困难的)给了我(-0.0755,-0.0159),对你来说(-0.0588,-0.0148),大约有20%的差异 . 这可能不是什么大问题,取决于你想要做什么 .

    我们也试试蛮力:

    sumfun <- function(x) {
        vv <- as.data.frame(VarCorr(x),order="lower.tri")[,"vcov"]
        ## cheating a bit here, using internal lme4 naming functions ...
        return(setNames(vv,
           c(lme4:::tnames(x,old=FALSE,prefix=c("var","cov")),
             "sigmasq")))
    }
    
    cc <- confint(s2,method="boot",nsim=1000,FUN=sumfun,seed=101,
            .progress="txt", PBargs=list(style=3))
    ## .progress/PBargs just cosmetic ...
    
    ##                                    2.5 %      97.5 %
    ## var_(Intercept)|teacher      0.079429623  0.24053633
    ## cov_occ.(Intercept)|teacher -0.067063911 -0.01479572
    ## var_occ|teacher              0.002733402  0.02378310
    ## sigmasq                      0.031952508  0.06736664
    

    这里的“黄金标准”似乎在我的 Profiles 结果和结果之间:协方差的下限是-0.067,而-0.0755( Profiles )或-0.0588 .

  • 1

    解决了,谢谢你的贡献 . 我会更新最初的帖子 . 可以使用Snijders&Bosker的数据集here测试结果 .

    导入

    library(foreign)
    chap12 <- read.dta(file = "<your path>/ch12.dta")
    

    即兴模型:

    snijders <- lmer(prox_pup ~ 1 + prox_sel + (1 + occ|teacher), data = chap12)
    

    来源功能:

    ExtractVarCovCI <- function(Model) {
    
    v <- NULL
    v <- as.data.frame(VarCorr(Model),order = "lower.tri") #Extract variances and covariances
    
    conf <- confint(Model, parm  ="theta_", oldNames = F) #extract CIs
    
    v.conf <- cbind(v,conf) #bind confidence intervals
    
    covs <- as.data.frame(v.conf[!is.na(v[,3]),]) #separate variance from covariance components
    vars <- as.data.frame(v.conf[is.na(v[,3]),]) #separate variance from covariance components
    vars.sq <- vars[,6:7]^2 #calculate square of variance components
    colnames(vars.sq) <- sub("[%]", "% sq.", colnames(vars.sq))
    
    vars2 <- cbind(vars,vars.sq) #bind squares of variance components
    covs$`2.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later
    covs$`97.5 % sq.` <- c(rep(NA,nrow(covs))) #create empty columns for later
    
    lcovs <- length(row.names(covs)) #now we re-organise the table so that each covariance is below the variance of its variables
    k <- NULL
    for (i in seq(1:lcovs)) {
      k <- rbind(k,vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,2],],vars2[vars2$grp %in% covs[i,1] & vars2$var1 %in% covs[i,3],],covs[i,])
    }
    
    k2 <- rbind(k,vars2["sigma",]) #bind the level-1 residuals at the end
    
    k2.covrow <- grep("^cor",rownames(k2)) # isolate covariance row position
    k2[k2.covrow,8] <- k2[k2.covrow,6]*k2[k2.covrow-1,6]*k2[k2.covrow-2,6] #calculate covariance 2.5%
    k2[k2.covrow,9] <- k2[k2.covrow,7]*k2[k2.covrow-1,7]*k2[k2.covrow-2,7] #calculate covariance 97.5%
    
    p <- NULL
    p <- k2[,c(4,8:9)] #retain only the estimates and the confidence intervals
    rownames(p) <- sub("^sd","var",rownames(p)) #now it's clear that we have proper variances and covariances
    rownames(p) <- sub("^cor","cov",rownames(p)) #now it's clear that we have proper variances and covariances
    colnames(p) <- c("Estimate", "2.5%", "97.5%")
    
    return(p)
    }
    

    运行功能:

    ExtractVarCovCI(snijders)
    

    我的输出是:

    Estimate         2.5%       97.5%
    var_(Intercept)|teacher      0.15617962  0.089020350  0.26130969
    var_occ|teacher              0.01205317  0.002467408  0.02779329
    cov_occ.(Intercept)|teacher -0.03883458 -0.014820577 -0.05887660
    sigma                        0.04979762  0.034631759  0.07263837
    

    现在我们有一个方差 - 协方差表,它使用非标准化的随机效应及其上限和下限置信区间 . 我相信有更好的方法可以做到这一点,但这是一个开始......

    ķ .

相关问题