Calculates model-averaged random effects std. dev.

sd-ranef-avg.R

sd.ranef.avg <-
function(object) {
    
    if(!inherits(object, "averaging"))
        stop("not an \"averaging\" object")

    mo <- get.models(object, TRUE)
    re <-  lapply(mo, function(x) as.data.frame(ranef(x, condVar = TRUE)))
    
    cols <- c("grpvar", "term", "grp")
    
    if(all(sapply(re, function(x) "component" %in% colnames(x))))
        cols <- c("component", cols)
    
    # check for RE consistency:
    if(!all(sapply(re[-1], function(x)
        all(re[[1]][, cols] == re[[2]][, cols])
        )))
        stop("inconsistent random effects between component models")
    
    if(any(re[[1L]][, "term"] != "(Intercept)"))
        warning("result is correct only for intercept random effects")

    wts <- Weights(object)
    ngroups <- nrow(re[[1]])
    avg.re <- array(NA_real_, dim = c(ngroups, 2L), dimnames = ) 
    for(i in 1:ngroups) {
        p <- do.call("rbind", lapply(re, "[", i, c("condval", "condsd")))
        avg.re[i, ] <- MuMIn::par.avg(p[, "condval"], p[, "condsd"], weight = wts)[1:2]
    }
    res <- re[[1]]
    res[, c("condval", "condsd")] <- avg.re
    #sqrt(var(res$condval) + res$condsd[1]^2)
    
    cols2 <- cols[-length(cols)]
    rval <- sapply(split(
        res[, c("condval", "condsd")],
        res[, cols2], drop = TRUE),
        function(x) sqrt(var(x$condval) + x$condsd[1]^2)
        )
    attr(rval, "ranef") <- res
    rval
}