Last active
May 13, 2022 12:37
-
-
Save k-barton/4b54e687dee6859eb10dfe4e27d4d5cd to your computer and use it in GitHub Desktop.
Calculates model-averaged random effects std. dev.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment