R/r.squaredGLMM.R In MuMIn: Multi-Model Inference

Defines functions ranform

````r.squaredGLMM` <-
function(x)
UseMethod("r.squaredGLMM")

`r.squaredGLMM.default` <-
function(x) .NotYetImplemented()

`r.squaredGLMM.lme` <-
function(x) {
VarFx <- var(fitted(x, level = 0L))
mmRE <- model.matrix(x\$modelStruct\$reStruct,
data = x\$data[rownames(x\$fitted), , drop = FALSE])
n <- nrow(mmRE)

sigma2 <- x\$sigma^2
reStruct <- x\$modelStruct\$reStruct
if((m <- length(reStruct)) > 1L) {
nams <- names(reStruct)
for(i in seq.int(m))
attr(reStruct[[i]], "Dimnames")[[2L]] <-
paste(nams[[i]], attr(reStruct[[i]], "Dimnames")[[2L]], sep = ".")
}

varRe <- sum(sapply(reStruct, function(z) {
sig <- nlme::pdMatrix(z) * sigma2
mm1 <-  mmRE[, rownames(sig), drop = FALSE]
#sum(diag(mm1 %*% sig %*% t(mm1))) / n
sum(matmultdiag(mm1 %*% sig, ty = mm1)) / n
}))

varTot <- sum(VarFx, varRe)
res <- c(VarFx, varTot) / (varTot + sigma2)
names(res) <- c("R2m", "R2c")
res
}

## extracts random effect formula. e.g:
ranform <- function (form) {
### XXX: would give an error: values must be length 1 ...
###      for very long RE formulas
ans <- update.formula(reformulate(vapply(lapply(.findbars(form),
"[[", 2L), deparse, "",  width.cutoff = 500L)), ~ . + 1)
environment(ans) <- environment(form)
ans
}

`r.squaredGLMM.merMod` <-
function(x) {
fam <- family(x)
useObsLevVar <- (fam\$family == "poisson" && fam\$link == "log") || fam\$family == "binomial"
## for poisson(log) and binomial(*), update 'x' to include individual-level
## variance (1 | 1:nobs(x)):
if (useObsLevVar && !any(sapply(x@flist, nlevels) == nobs(x))) {
cl <- get_call(x)
frm <- formula(x)
nRowData <- eval(call("eval", as.expression(call("NROW", cl\$formula[[2L]])),
envir = cl\$data), envir = environment(frm),
enclos = parent.frame())
fl <- length(frm)
frx <- . ~ . + 1
frx[[3L]][[3L]] <- call("(", call("|", 1, call("gl", nRowData, 1)))
cl\$formula <- update.formula(frm, frx)
x <- tryCatch(eval(cl, envir = environment(frm), enclos = parent.frame()),
error = function(e) {
cry(conditionCall(e), conditionMessage(e), warn = TRUE)
cry(cl, "fitting model with the observation-level random effect term failed. Add the term manually")
})
message("The result is correct only if all data used by the model ",
"has not changed since model was fitted.", domain = "R-MuMIn")
}

mmAll <- model.matrix(ranform(formula(x)), data = model.frame(x))
##Note: Argument 'contrasts' can only be specified for fixed effects
##contrasts.arg = eval(cl\$contrasts, envir = environment(formula(x))))

vc <- lme4::VarCorr(x)

n <- nrow(mmAll)
fx <- lme4::fixef(x) # fixed effect estimates
fxpred <- as.vector(model.matrix(x) %*% fx)

if (useObsLevVar) {
vname <- names(x@flist)[sapply(x@flist, nlevels) == n][1L]

if(! vname %in% names(vc)) vname <- make.names(vname)
stopifnot(vname %in% names(vc))
varResid <- vc[[vname]][1L]
beta0 <- mean(fxpred)
vc <- vc[names(vc) != vname]
} else {
varResid <- attr(vc, "sc")^2
beta0 <- NULL
}

if(!all(unlist(sapply(vc, rownames), use.names = FALSE) %in% colnames(mmAll)))
stop("random term names do not match those in model matrix. \n",
"Have 'options(contrasts)' changed since the model was fitted?")

varRe <- if(length(vc) == 0L) 0L else
sum(sapply(vc, function(sig) {
mm1 <-  mmAll[, rownames(sig), drop = FALSE]
# sum(matmult(mm1 %*% sig, t(mm1), diag.only = TRUE)) / n
sum(matmultdiag(mm1 %*% sig, ty = mm1)) / n
#sum(diag(mm1 %*% sig %*% t(mm1))) / n
}))

#varRe0 <- if(length(vc) == 0L) 0L else
#          sum(sapply(vc, function(sig) sig[[1]]))

.rsqGLMM(fam = family(x), varFx = var(fxpred), varRe = varRe,
varResid = varResid, beta0 = beta0)
}

`r.squaredGLMM.glmmML` <-
function(x) {
if(is.null(x\$x))
stop("glmmML must be fitted with 'x = TRUE'")

fam <- family(x)
useObsLevVar <- (fam\$family == "poisson" && fam\$link == "log") || fam\$family == "binomial"
if(useObsLevVar) {
cry(, "cannot calculate 'unit variance' in glmmML")
}
fxpred <- as.vector(x\$x %*% coef(x))
.rsqGLMM(family(x), varFx = var(fxpred), varRe = x\$sigma^2, varResid = NULL,
beta0 = mean(fxpred))
}

`r.squaredGLMM.lm` <-
function(x) {
fam <- family(x)
.rsqGLMM(fam,
varFx = var(as.vector(model.matrix(x) %*% coef(x))),
#varFx = var(fitted(x)),
varRe = 0,
varResid = sum(if(is.null(x\$weights)) resid(x)^2 else
resid(x)^2 * x\$weights) / df.residual(x),
beta0 = if(fam\$family == "poisson" && fam\$link == "log")
log(mean(model.response(model.frame(x)))) else
NULL
)
}

`.rsqGLMM` <-
function (fam, varFx, varRe, varResid, beta0) {
varDistr <- switch(paste(fam\$family, fam\$link, sep = "."),
gaussian.identity = 0,
binomial.logit = 3.28986813369645,
binomial.probit = 1,
poisson.log = {
expBeta0 <- exp(beta0)
if (expBeta0 < 6) cry(sys.call(-1L), "exp(beta0) of %0.1f is too close to zero, estimate may be unreliable \n",
expBeta0, warn = TRUE)
log1p(1 / expBeta0)
},
poisson.sqrt = 0.25,
cry(sys.call(-1L), "do not know how to calculate variance for this family/link combination")
)

#print(c(Sf = varFx, Sl = varRe, Se = varResid, Sd = varDistr))
#  total.var <- [Sf + Sl] + [Se + Sd]
varTot <- sum(varFx, varRe)
res <- c(varFx, varTot) / (varTot + varDistr + varResid)
names(res) <- c("R2m", "R2c")
res
}
```

Try the MuMIn package in your browser

Any scripts or data that you put into this service are public.

MuMIn documentation built on Jan. 31, 2018, 1:03 a.m.