Nothing
## Type II and III tests for linear, generalized linear, and other models (J. Fox)
## most of what's below is copied from car::Anova.R
## main changes are (1) absence of F-test (K-R, Satterthwaite df) capability;
## (2) use of [[component]] to pick out relevant fixed effect parameters/v-cov matrix
## copied unchanged (?); unexported utilities from car
responseName.default <- function (model, ...) deparse(attr(terms(model), "variables")[[2]])
term.names.default <- function (model, component="cond", ...) {
term.names <- labels(terms(model, component=component))
if (has.intercept(model)) c("(Intercept)", term.names)
else term.names
}
has.intercept <- function (model, ...) {
UseMethod("has.intercept")
}
ConjComp <- function(X, Z = diag( nrow(X)), ip = diag(nrow(X))) {
## This function by Georges Monette
## finds the conjugate complement of the proj of X in span(Z) wrt
## inner product ip
## - assumes Z is of full column rank
## - projects X conjugately wrt ip into span Z
xq <- qr(t(Z) %*% ip %*% X)
if (xq$rank == 0) return(Z)
Z %*% qr.Q(xq, complete = TRUE) [ ,-(1:xq$rank)]
}
relatives <- function(term, names, factors){
is.relative <- function(term1, term2) {
all(!(factors[,term1]&(!factors[,term2])))
}
if(length(names) == 1) return(NULL)
which.term <- which(term==names)
(1:length(names))[-which.term][sapply(names[-which.term],
function(term2) is.relative(term, term2))]
}
## modified
has.intercept.glmmTMB <- function (model, component="cond", ...) {
nms <- names(fixef(model)[[component]])
any(grepl("\\(Intercept\\)",nms))
}
## n.b. rawNamespace spec must now be on a single line!
##' @rdname downstream_methods
##' @rawNamespace if(getRversion() >= "3.6.0") { S3method(car::Anova, glmmTMB) } else { export(Anova.glmmTMB) }
##' @param vcov. variance-covariance matrix (usually extracted automatically)
##' @param test.statistic unused: only valid choice is "Chisq" (i.e., Wald chi-squared test)
##' @param singular.ok OK to do ANOVA with singular models (unused) ?
##' @param type type of test, \code{"II"}, \code{"III"}, \code{2}, or \code{3}. Roman numerals are equivalent to the corresponding Arabic numerals. See \code{\link[car]{Anova}} for details.
##' @param include.rankdef.cols include all columns of a rank-deficient model matrix?
Anova.glmmTMB <- function (mod, type = c("II", "III", 2, 3),
test.statistic = c("Chisq","F"),
component="cond",
vcov. = vcov(mod)[[component]],
singular.ok,
include.rankdef.cols = FALSE,
...) {
ff <- fixef(mod)[[component]]
if (trivialFixef(names(ff),component)) {
stop(sprintf("trivial fixed effect for component %s: can't compute Anova table", sQuote(component)))
}
test.statistic <- match.arg(test.statistic)
if (test.statistic=="F") {
stop("F tests currently unavailable")
}
if (is.function(vcov.))
vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
afun <- switch(type,
`2` = , II = Anova.II.glmmTMB,
`3` = , III = Anova.III.glmmTMB)
afun(mod, vcov., test=test.statistic, singular.ok = singular.ok,
component = component, include.rankdef.cols = include.rankdef.cols)
}
## defined as a function, not a method, so we can hand the object
## off to car::linearHypothesis.default (not exported)
linearHypothesis_glmmTMB <- function (model, hypothesis.matrix,
rhs = NULL, test = c("Chisq", "F"),
vcov. = NULL, singular.ok = FALSE, verbose = FALSE,
coef. = NULL, component="cond", ...)
{
## what's the least ugly way to do this?
## match.call?
test <- match.arg(test)
## call linearHypothesis.default (not exported)
if (!requireNamespace("car")) {
stop("please install (if necessary) and load the car package")
}
if (utils::packageVersion("car")<"3.0.6") {
stop("please install a more recent version of the car package (>= 3.0.6)")
}
car::linearHypothesis(model=model,
hypothesis.matrix=hypothesis.matrix,
rhs=rhs,
test=test,
vcov. = vcov.,
singular.ok = singular.ok,
verbose = verbose,
coef. = fixef(model)[[component]],
...)
}
Anova.II.glmmTMB <- function(mod, vcov., singular.ok=TRUE, test="Chisq",
component="cond", include.rankdef.cols = FALSE, ...){
## would feel cleaner to have this external, but it uses
## lots of variable from the function environment ...
hyp.term <- function(term) {
which.term <- which(term==names)
subs.term <- which(assign==which.term)
relatives <- relatives(term, names, fac)
subs.relatives <- NULL
for (relative in relatives)
subs.relatives <- c(subs.relatives, which(assign==relative))
hyp.matrix.1 <- I.p[subs.relatives,,drop=FALSE]
hyp.matrix.1 <- hyp.matrix.1[, not.aliased, drop=FALSE]
hyp.matrix.2 <- I.p[c(subs.relatives,subs.term),,drop=FALSE]
hyp.matrix.2 <- hyp.matrix.2[, not.aliased, drop=FALSE]
hyp.matrix.term <- if (nrow(hyp.matrix.1) == 0) {
hyp.matrix.2
} else {
t(ConjComp(t(hyp.matrix.1),
t(hyp.matrix.2), vcov.))
}
hyp.matrix.term <- hyp.matrix.term[!apply(hyp.matrix.term, 1,
function(x) all(x == 0)), , drop=FALSE]
if (nrow(hyp.matrix.term) == 0)
return(c(statistic=NA, df=0))
hyp <- linearHypothesis_glmmTMB(mod, hyp.matrix.term,
vcov.=vcov.,
singular.ok=singular.ok,
test=test,
component=component, ...)
if (test == "Chisq") return(c(statistic=hyp$Chisq[2], df=hyp$Df[2]))
else return(c(statistic=hyp$F[2], df=hyp$Df[2], res.df=hyp$Res.Df[2]))
} ## hyp.term()
not.aliased <- !is.na(fixef(mod)[[component]])
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod, component=component), "factors")
intercept <- has.intercept(mod)
p <- length(fixef(mod)[[component]])
I.p <- diag(p)
## FIXME:: missing or !missing ???
if (missing(vcov.)){
vcov. <- vcov(mod, complete=FALSE)[[component]]
}
vcov. <- vcov.[not.aliased, not.aliased]
assign <- attr(model.matrix(mod, component=component, include_rankdef = include.rankdef.cols), "assign")
assign[!not.aliased] <- NA
names <- term.names.default(mod, component=component)
if (intercept) names <- names[-1]
n.terms <- length(names)
p <- teststat <- df <- res.df <- rep(0, n.terms)
for (i in seq_len(n.terms)) {
hyp <- hyp.term(names[i])
teststat[i] <- abs(hyp["statistic"])
df[i] <- abs(hyp["df"])
res.df[i] <- hyp["res.df"]
p[i] <- pchisq(teststat[i], df[i], lower.tail=FALSE)
}
result <- data.frame(teststat, df, p)
row.names(result) <- names
names(result) <- c ("Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II Wald chisquare tests)\n",
paste("Response:", responseName.default(mod)))
return(result)
}
Anova.III.glmmTMB <- function(mod, vcov., singular.ok=FALSE, test="Chisq",
component="cond", include.rankdef.cols = include.rankdef.cols, ...){
intercept <- has.intercept(mod)
p <- length(fixef(mod)[[component]])
I.p <- diag(p)
names <- term.names.default(mod, component=component)
n.terms <- length(names)
assign <- attr(model.matrix(mod, component=component, include.rankdef.cols = include.rankdef.cols), "assign")
p <- teststat <- df <- res.df <- rep(0, n.terms)
if (intercept) df[1] <- 1
not.aliased <- !is.na(fixef(mod)[[component]])
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
if (missing(vcov.)){
vcov. <- vcov(mod, complete=FALSE)[[component]]
}
vcov. <- vcov.[not.aliased, not.aliased]
assign <- attr(model.matrix(mod, component=component, include_rankdef = include.rankdef.cols), "assign")
assign[!not.aliased] <- NA
for (term in seq_len(n.terms)){
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs,,drop=FALSE]
hyp.matrix <- hyp.matrix[, not.aliased, drop=FALSE]
hyp.matrix <- hyp.matrix[!apply(hyp.matrix, 1, function(x) all(x == 0)), , drop=FALSE]
if (nrow(hyp.matrix) == 0){
teststat[term] <- NA
df[term] <- 0
p[term] <- NA
}
else {
hyp <- linearHypothesis_glmmTMB(mod, hyp.matrix, test=test,
vcov.=vcov., singular.ok=singular.ok,
component=component, ...)
teststat[term] <- hyp$Chisq[2]
df[term] <- abs(hyp$Df[2])
p[term] <- pchisq(teststat[term], df[term], lower.tail=FALSE)
}
result <- data.frame(teststat, df, p)
row.names(result) <- names
names(result) <- c ("Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III Wald chisquare tests)\n",
paste("Response:", responseName.default(mod)))
}
result
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.