R/Anova.R
In jonathon-love/car3: Companion to Applied Regression
#-------------------------------------------------------------------------------
# Revision history:
# 2009-01-05: bug fix in Anova.II.lm(). J. Fox
# 2009-01-16: Cox models with clusters now handled. J. Fox
# 2009-09-16: reworked glm and lm methods to handle aliased parameters. J. Fox
# 2009-09-30: renamed "Anova" to "Analysis of Deviance" in output for some methods. J. Fox
# 2009-12-22: modified Anova.mlm() to handle a user-supplied within-subject model matrix. J. Fox
# 2009-12-28: named the components of P in Anova.III.mlm(). John
# 2010-01-01: Anova.II.mlm() now hands off (again) to Anova.III.mlm() when there
# is only an intercept in the between-subjects model
# 2010-02-17: Fixed bug that caused some models with aliased coefficients to fail. J. Fox
# 2010-06-14: added wcrossprod and allow use of observation weights in Anova.mlm()
# 2010-06-28: Fixed Anova() tables for coxph and survreg models
# (failed because of changes in survival package.
# 2011-01-21: Added functions for mixed models. J. Fox
# 2011-01-25: Fixed Anova.polr() and Anova.multinom() to work with models with only one term. J. Fox
# 2011-05-19: local fixef() to avoid nlme/lme4 issues. J. Fox
# 2011-05-11: changed order of columns in ANOVA tables for mixed models. J. Fox
# 2011-11-27: added Anova.svyglm(). J. Fox
# 2011-12-31: fixed bug in Anova.II(and III).F.glm() when na.exclude used. J. Fox
# 2012-02-28: added test.statistic argument to Anova.mer(). J.Fox
# 2012-03-02: fixed test abbreviation of test.statistic argument to Anova.default()
# called by other Anova() methods. J. Fox
# 2013-06-17: modified summary.Anova.mlm(), introduced print.summary.Anova.mlm(),
# adapting code contributed by Gabriel Baud-Bovy. J. Fox
# 2013-06-20: added Anova.merMod() method. J. Fox
# 2013-06-22: tweaks to local fixef(). J. Fox
# 2013-06-22: test argument uniformly uses "Chisq" rather than "chisq". J. Fox
# 2013-08-19: replaced calls to print.anova(). J. Fox
# 2014-08-17: added calls to requireNamespace() and :: where needed (doesn't work for pbkrtest). J. Fox
# 2014-08-18: fixed bugs in Anova.survreg() for types II, III LR tests and Wald tests. J. Fox
# 2014-09-23: added Anova.rlm(). J. Fox
# 2014-10-10: removed MASS:: from calls to polr(). John
# 2014-12-18: check that residual df and SS are nonzero in Anova.lm(). John
# 2015-01-27: vcovAdj() and methods now imported from pbkrtest. John
# 2015-02-18: force evaluation of vcov. when it's a function. John
# 2015-04-30: don't allow error.estimate="dispersion" for F-tests in binomial
# and Poission GLMs. John
# 2015-08-29: fixed Anova() for coxph models with clusters. John
# 2015-09-04: added support for coxme models. John
# 2015-09-11: modified Anova.default() to work with vglm objects from VGAM. John
# 2015-09-15: fixed Anova.default() so that F-tests work again. John
# 2015-11-13: modify Anova.coxph() to take account of method/ties argument. John
# 2016-06-03: added SSP and SSPE args to print.summary.Anova.mlm(). John
# 2016-06-25: added code to optionally print univariate ANOVAs for a mlm. John
#-------------------------------------------------------------------------------
# Type II and III tests for linear, generalized linear, and other models (J. Fox)
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))]
}
Anova <- function(mod, ...){
UseMethod("Anova", mod)
}
# linear models
Anova.lm <- function(mod, error, type=c("II","III", 2, 3),
white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4"),
vcov.=NULL, singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
if (df.residual(mod) == 0) stop("residual df = 0")
if (deviance(mod) < sqrt(.Machine$double.eps)) stop("residual sum of squares is 0 (within rounding error)")
type <- as.character(type)
white.adjust <- as.character(white.adjust)
type <- match.arg(type)
white.adjust <- match.arg(white.adjust)
if (missing(singular.ok)){
singular.ok <- type == "2" || type == "II"
}
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
if (white.adjust != "FALSE"){
if (white.adjust == "TRUE") white.adjust <- "hc3"
return(Anova.default(mod, type=type, vcov.=hccm(mod, type=white.adjust), test.statistic="F",
singular.ok=singular.ok, ...))
}
else if (!is.null(vcov.)) return(Anova.default(mod, type=type, vcov.=vcov., test.statistic="F",
singular.ok=singular.ok, ...))
switch(type,
II=Anova.II.lm(mod, error, singular.ok=singular.ok, ...),
III=Anova.III.lm(mod, error, singular.ok=singular.ok, ...),
"2"=Anova.II.lm(mod, error, singular.ok=singular.ok, ...),
"3"=Anova.III.lm(mod, error, singular.ok=singular.ok,...))
}
Anova.aov <- function(mod, ...){
class(mod) <- "lm"
Anova.lm(mod, ...)
}
Anova.II.lm <- function(mod, error, singular.ok=TRUE, ...){
if (!missing(error)){
sumry <- summary(error, corr=FALSE)
s2 <- sumry$sigma^2
error.df <- error$df.residual
error.SS <- s2*error.df
}
SS.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(mod)))
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(SS=NA, df=0))
lh <- linearHypothesis(mod, hyp.matrix.term,
singular.ok=singular.ok, ...)
abs(c(SS=lh$"Sum of Sq"[2], df=lh$Df[2]))
}
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
I.p <- diag(length(coefficients(mod)))
assign <- mod$assign
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <-names[-1]
n.terms <- length(names)
p <- df <- f <- SS <- rep(0, n.terms + 1)
sumry <- summary(mod, corr = FALSE)
SS[n.terms + 1] <- if (missing(error)) sumry$sigma^2*mod$df.residual
else error.SS
df[n.terms + 1] <- if (missing(error)) mod$df.residual else error.df
p[n.terms + 1] <- f[n.terms + 1] <- NA
for (i in 1:n.terms){
ss <- SS.term(names[i])
SS[i] <- ss["SS"]
df[i] <- ss["df"]
f[i] <- df[n.terms+1]*SS[i]/(df[i]*SS[n.terms + 1])
p[i] <- pf(f[i], df[i], df[n.terms + 1], lower.tail = FALSE)
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(names,"Residuals")
names(result) <- c("Sum Sq", "Df", "F value", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Anova Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
# type III
Anova.III.lm <- function(mod, error, singular.ok=FALSE, ...){
if (!missing(error)){
error.df <- df.residual(error)
error.SS <- deviance(error)
}
else {
error.df <- df.residual(mod)
error.SS <- deviance(mod)
}
intercept <- has.intercept(mod)
I.p <- diag(length(coefficients(mod)))
Source <- term.names(mod)
n.terms <- length(Source)
p <- df <- f <- SS <- rep(0, n.terms + 1)
assign <- mod$assign
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1: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){
SS[term] <- NA
df[term] <- 0
f[term] <- NA
p[term] <- NA
}
else {
test <- if (missing(error)) linearHypothesis(mod, hyp.matrix,
singular.ok=singular.ok, ...)
else linearHypothesis(mod, hyp.matrix, error.SS=error.SS, error.df=error.df,
singular.ok=singular.ok, ...)
SS[term] <- test$"Sum of Sq"[2]
df[term] <- test$"Df"[2]
f[term] <- test$"F"[2]
p[term] <- test$"Pr(>F)"[2]
}
}
Source[n.terms + 1] <- "Residuals"
SS[n.terms + 1] <- error.SS
df[n.terms + 1] <- error.df
p[n.terms + 1] <- f[n.terms + 1] <- NA
result <- data.frame(SS, df, f, p)
row.names(result) <- Source
names(result) <- c("Sum Sq", "Df", "F value", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Anova Table (Type III tests)\n", paste("Response:", responseName(mod)))
result
}
# generalized linear models
Anova.glm <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald", "F"),
error, error.estimate=c("pearson", "dispersion", "deviance"), singular.ok, ...){
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
if (missing(singular.ok)){
singular.ok <- type == "2" || type == "II"
}
test.statistic <- match.arg(test.statistic)
error.estimate <- match.arg(error.estimate)
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="II", singular.ok=singular.ok),
F=Anova.II.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
III=switch(test.statistic,
LR=Anova.III.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="III", singular.ok=singular.ok),
F=Anova.III.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
"2"=switch(test.statistic,
LR=Anova.II.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="II", singular.ok=singular.ok),
F=Anova.II.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
"3"=switch(test.statistic,
LR=Anova.III.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="III", singular.ok=singular.ok),
F=Anova.III.F.glm(mod, error, error.estimate, singular.ok=singular.ok)))
}
# type III
# LR test
Anova.III.LR.glm <- function(mod, singular.ok=FALSE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
Source <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(Source)
p <- df <- LR <- rep(0, n.terms)
dispersion <- summary(mod, corr = FALSE)$dispersion
deviance <- deviance(mod)/dispersion
for (term in 1:n.terms){
mod.1 <- drop1(mod, scope=eval(parse(text=paste("~",Source[term]))))
LR[term] <- (mod.1$Deviance[2]/dispersion)-deviance
df[term] <- mod.1$Df[2]
p[term] <- pchisq(LR[term], df[term], lower.tail = FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- Source
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova","data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n", paste("Response:", responseName(mod)))
result
}
# F test
Anova.III.F.glm <- function(mod, error, error.estimate, singular.ok=FALSE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
fam <- family(mod)$family
if ((fam == "binomial" || fam == "poisson") && error.estimate == "dispersion"){
warning("dispersion parameter estimated from the Pearson residuals, not taken as 1")
error.estimate <- "pearson"
}
if (missing(error)) error <- mod
df.res <- df.residual(error)
error.SS <- switch(error.estimate,
pearson=sum(residuals(error, "pearson")^2, na.rm=TRUE),
dispersion=df.res*summary(error, corr = FALSE)$dispersion,
deviance=deviance(error))
Source <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(Source)
p <- df <- f <- SS <-rep(0, n.terms+1)
f[n.terms+1] <- p[n.terms+1] <- NA
df[n.terms+1] <- df.res
SS[n.terms+1] <- error.SS
dispersion <- error.SS/df.res
deviance <- deviance(mod)
for (term in 1:n.terms){
mod.1 <- drop1(mod, scope=eval(parse(text=paste("~",Source[term]))))
df[term] <- mod.1$Df[2]
SS[term] <- mod.1$Deviance[2] - deviance
f[term] <- (SS[term]/df[term])/dispersion
p[term] <- pf(f[term], df[term], df.res, lower.tail = FALSE)
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(Source, "Residuals")
names(result) <- c("SS", "Df", "F", "Pr(>F)")
class(result) <- c("anova","data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)),
paste("Error estimate based on",
switch(error.estimate,
pearson="Pearson residuals", dispersion="estimated dispersion",
deviance="deviance"), "\n"))
result
}
# type II
# LR test
Anova.II.LR.glm <- function(mod, singular.ok=TRUE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
# (some code adapted from drop1.glm)
which.nms <- function(name) which(asgn == which(names == name))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- mod$y
if (is.null(y)) y <- model.response(model.frame(mod), "numeric")
wt <- mod$prior.weights
if (is.null(wt)) wt <- rep(1, length(y))
asgn <- attr(X, 'assign')
df <- p <- LR <- rep(0, n.terms)
dispersion <- summary(mod, corr = FALSE)$dispersion
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- glm.fit(X[, -exclude.1, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
glm.fit(X[, -exclude.2, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
}
dev.2 <- deviance(mod.2)
df[term] <- df.residual(mod.1) - df.residual(mod.2)
if (df[term] == 0) LR[term] <- p[term] <- NA
else {
LR[term] <- (dev.1 - dev.2)/dispersion
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <-
c("Analysis of Deviance Table (Type II tests)\n", paste("Response:", responseName(mod)))
result
}
# F test
Anova.II.F.glm <- function(mod, error, error.estimate, singular.ok=TRUE, ...){
# (some code adapted from drop1.glm)
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
fam <- family(mod)$family
if ((fam == "binomial" || fam == "poisson") && error.estimate == "dispersion"){
warning("dispersion parameter estimated from the Pearson residuals, not taken as 1")
error.estimate <- "pearson"
}
which.nms <- function(name) which(asgn == which(names == name))
if (missing(error)) error <- mod
df.res <- df.residual(error)
error.SS <- switch(error.estimate,
pearson = sum(residuals(error, "pearson")^2, na.rm=TRUE),
dispersion = df.res*summary(error, corr = FALSE)$dispersion,
deviance = deviance(error))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- mod$y
if (is.null(y)) y <- model.response(model.frame(mod), "numeric")
wt <- mod$prior.weights
if (is.null(wt)) wt <- rep(1, length(y))
asgn <- attr(X, 'assign')
p <- df <- f <- SS <- rep(0, n.terms+1)
f[n.terms+1] <- p[n.terms+1] <- NA
df[n.terms+1] <- df.res
SS[n.terms+1] <- error.SS
dispersion <- error.SS/df.res
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- glm.fit(X[, -exclude.1, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
glm.fit(X[, -exclude.2, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
}
dev.2 <- deviance(mod.2)
df[term] <- df.residual(mod.1) - df.residual(mod.2)
if (df[term] == 0) SS[term] <- f[term] <- p[term] <- NA
else {
SS[term] <- dev.1 - dev.2
f[term] <- SS[term]/(dispersion*df[term])
p[term] <- pf(f[term], df[term], df.res, lower.tail=FALSE)
}
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(names, "Residuals")
names(result) <- c("SS", "Df", "F", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)),
paste("Error estimate based on",
switch(error.estimate,
pearson="Pearson residuals",
dispersion="estimated dispersion",
deviance="deviance"), "\n"))
result
}
# multinomial logit models (via multinom in the nnet package)
Anova.multinom <-
function (mod, type = c("II", "III", 2, 3), ...)
{
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
switch(type,
II = Anova.II.multinom(mod, ...),
III = Anova.III.multinom(mod, ...),
"2" = Anova.II.multinom(mod, ...),
"3" = Anova.III.multinom(mod, ...))
}
Anova.II.multinom <- function (mod, ...)
{
which.nms <- function(name) which(asgn == which(names ==
name))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- mod$weights
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms) {
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels),
which.nms)))
mod.1 <-if (n.terms > 1) multinom(y ~ X[, -c(1, exclude.1)], weights=wt, trace=FALSE)
else multinom(y ~ 1, weights=wt, race=FALSE)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0)
mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
multinom(y ~ X[, -c(1, exclude.2)], weights=wt, trace=FALSE)
}
dev.2 <- deviance(mod.2)
LR[term] <- dev.1 - dev.2
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.multinom <- function (mod, ...)
{
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- mod$weights
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
deviance <- deviance(mod)
for (term in 1:n.terms) {
mod.1 <- if (n.terms > 1) multinom(y ~ X[, term != asgn][, -1], weights=wt, trace=FALSE)
else multinom(y ~ 1, weights=wt, trace=FALSE)
LR[term] <- deviance(mod.1) - deviance
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
# proportional-odds logit models (via polr in the MASS package)
Anova.polr <- function (mod, type = c("II", "III", 2, 3), ...)
{
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
switch(type,
II = Anova.II.polr(mod, ...),
III = Anova.III.polr(mod, ...),
"2" = Anova.II.polr(mod, ...),
"3" = Anova.III.polr(mod, ...))
}
Anova.II.polr <- function (mod, ...)
{
if (!requireNamespace("MASS")) stop("MASS package is missing")
which.nms <- function(name) which(asgn == which(names ==
name))
fac <- attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- model.weights(model.frame(mod))
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms) {
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels),
which.nms)))
mod.1 <- if (n.terms > 1) polr(y ~ X[, -c(1, exclude.1)], weights=wt)
else polr(y ~ 1, weights=wt)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0)
mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
polr(y ~ X[, -c(1, exclude.2)], weights=wt)
}
dev.2 <- deviance(mod.2)
LR[term] <- dev.1 - dev.2
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.polr <- function (mod, ...)
{
if (!requireNamespace("MASS")) stop("MASS package is missing")
names <- term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- model.weights(model.frame(mod))
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
deviance <- deviance(mod)
for (term in 1:n.terms) {
mod.1 <- if (n.terms > 1) polr(y ~ X[, term != asgn][, -1], weights=wt)
else polr(y ~ 1, weights=wt)
LR[term] <- deviance(mod.1) - deviance
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
# multivariate linear models
# the following 3 functions copied from the stats package (not exported from stats)
Pillai <- function (eig, q, df.res) {
test <- sum(eig/(1 + eig))
p <- length(eig)
s <- min(p, q)
n <- 0.5 * (df.res - p - 1)
m <- 0.5 * (abs(p - q) - 1)
tmp1 <- 2 * m + s + 1
tmp2 <- 2 * n + s + 1
c(test, (tmp2/tmp1 * test)/(s - test), s * tmp1, s * tmp2)
}
Wilks <- function (eig, q, df.res) {
test <- prod(1/(1 + eig))
p <- length(eig)
tmp1 <- df.res - 0.5 * (p - q + 1)
tmp2 <- (p * q - 2)/4
tmp3 <- p^2 + q^2 - 5
tmp3 <- if (tmp3 > 0)
sqrt(((p * q)^2 - 4)/tmp3)
else 1
c(test, ((test^(-1/tmp3) - 1) * (tmp1 * tmp3 - 2 * tmp2))/p/q,
p * q, tmp1 * tmp3 - 2 * tmp2)
}
HL <- function (eig, q, df.res) {
test <- sum(eig)
p <- length(eig)
m <- 0.5 * (abs(p - q) - 1)
n <- 0.5 * (df.res - p - 1)
s <- min(p, q)
tmp1 <- 2 * m + s + 1
tmp2 <- 2 * (s * n + 1)
c(test, (tmp2 * test)/s/s/tmp1, s * tmp1, tmp2)
}
Roy <- function (eig, q, df.res) {
p <- length(eig)
test <- max(eig)
tmp1 <- max(p, q)
tmp2 <- df.res - tmp1 + q
c(test, (tmp2 * test)/tmp1, tmp1, tmp2)
}
has.intercept.mlm <- function (model, ...)
any(row.names(coefficients(model)) == "(Intercept)")
Anova.mlm <- function(mod, type=c("II","III", 2, 3), SSPE, error.df, idata,
idesign, icontrasts=c("contr.sum", "contr.poly"), imatrix,
test.statistic=c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"),...){
wts <- if (!is.null(mod$weights)) mod$weights else rep(1, nrow(model.matrix(mod)))
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(SSPE)) SSPE <- wcrossprod(residuals(mod), w=wts)
if (missing(idata)) {
idata <- NULL
idesign <- NULL
}
if (missing(imatrix)) imatrix <- NULL
error.df <- if (missing(error.df)) df.residual(mod)
else error.df
switch(type,
II=Anova.II.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
III=Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
"2"=Anova.II.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
"3"=Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...))
}
Anova.III.mlm <- function(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...){
intercept <- has.intercept(mod)
V <- solve(crossprod(model.matrix(mod)))
p <- nrow(coefficients(mod))
I.p <- diag(p)
terms <- term.names(mod)
n.terms <- length(terms)
assign <- mod$assign
if (is.null(idata) && is.null(imatrix)){
if ((n.terms == 0) && intercept) {
Test <- linearHypothesis(mod, 1, SSPE=SSPE, ...)
result <- list(SSP=Test$SSPH, SSPE=SSPE, df=1, error.df=error.df,
terms="(Intercept)", repeated=FALSE, type="III", test=test)
class(result) <- "Anova.mlm"
return(result)
}
SSP <- as.list(rep(0, n.terms))
df <- rep(0, n.terms)
names(df) <- names(SSP) <- terms
for (term in 1:n.terms){
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs,,drop=FALSE]
Test <- linearHypothesis(mod, hyp.matrix, SSPE=SSPE, ...)
SSP[[term]] <- Test$SSPH
df[term]<- length(subs)
}
result <- list(SSP=SSP, SSPE=SSPE, df=df, error.df=error.df, terms=terms,
repeated=FALSE, type="III", test=test)
}
else {
if (!is.null(imatrix)){
X.design <- do.call(cbind, imatrix)
ncols <- sapply(imatrix, ncol)
end <- cumsum(ncols)
start <- c(1, (end + 1))[-(length(end) + 1)]
cols <- mapply(seq, from=start, to=end)
iterms <- names(end)
names(cols) <- iterms
check.imatrix(X.design, iterms)
}
else {
if (is.null(idesign)) stop("idesign (intra-subject design) missing.")
for (i in 1:length(idata)){
if (is.null(attr(idata[,i], "contrasts"))){
contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2]
else icontrasts[1]
}
}
X.design <- model.matrix(idesign, data=idata)
i.intercept <- has.intercept(X.design)
iterms <- term.names(idesign)
if (i.intercept) iterms <- c("(Intercept)", iterms)
check.imatrix(X.design)
}
df <- rep(0, n.terms*length(iterms))
hnames <- rep("", length(df))
P <- SSPEH <- SSP <- as.list(df)
i <- 0
for (iterm in iterms){
for (term in 1:n.terms){
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs,,drop=FALSE]
i <- i + 1
Test <- linearHypothesis(mod, hyp.matrix, SSPE=SSPE,
idata=idata, idesign=idesign, icontrasts=icontrasts, iterms=iterm,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)
SSP[[i]] <- Test$SSPH
SSPEH[[i]] <- Test$SSPE
P[[i]] <- Test$P
df[i] <- length(subs)
hnames[i] <- if (iterm == "(Intercept)") terms[term]
else if (terms[term] == "(Intercept)") iterm
else paste(terms[term], ":", iterm, sep="")
}
}
names(df) <- names(SSP) <- names(SSPEH) <- names(P) <- hnames
result <- list(SSP=SSP, SSPE=SSPEH, P=P, df=df, error.df=error.df,
terms=hnames, repeated=TRUE, type="III", test=test,
idata=idata, idesign=idesign, icontrasts=icontrasts, imatrix=imatrix,
singular=Test$singular)
}
class(result) <- "Anova.mlm"
result
}
Anova.II.mlm <- function(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...){
wts <- if (!is.null(mod$weights)) mod$weights else rep(1, nrow(model.matrix(mod)))
V <- solve(wcrossprod(model.matrix(mod), w=wts))
SSP.term <- function(term, iterm){
which.term <- which(term == terms)
subs.term <- which(assign == which.term)
relatives <- relatives(term, terms, 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.2 <- I.p[c(subs.relatives, subs.term),,drop=FALSE]
if (missing(iterm)){
SSP1 <- if (length(subs.relatives) == 0) 0
else linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V, singular.ok=TRUE, ...)$SSPH
SSP2 <- linearHypothesis(mod, hyp.matrix.2, SSPE=SSPE, V=V, singular.ok=TRUE, ...)$SSPH
return(SSP2 - SSP1)
}
else {
SSP1 <- if (length(subs.relatives) == 0) 0
else linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts, P=imatrix[[iterm]], singular.ok=TRUE, ...)$SSPH
lh2 <- linearHypothesis(mod, hyp.matrix.2, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts, P=imatrix[[iterm]], singular.ok=TRUE, ...)
return(list(SSP = lh2$SSPH - SSP1, SSPE=lh2$SSPE, P=lh2$P, singular=lh2$singular))
}
}
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- nrow(coefficients(mod))
I.p <- diag(p)
assign <- mod$assign
terms <- term.names(mod)
if (intercept) terms <- terms[-1]
n.terms <- length(terms)
if (n.terms == 0){
message("Note: model has only an intercept; equivalent type-III tests substituted.")
return(Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...))
}
if (is.null(idata) && is.null(imatrix)){
SSP <- as.list(rep(0, n.terms))
df <- rep(0, n.terms)
names(df) <- names(SSP) <- terms
for (i in 1:n.terms){
SSP[[i]] <- SSP.term(terms[i])
df[i]<- df.terms(mod, terms[i])
}
result <- list(SSP=SSP, SSPE=SSPE, df=df, error.df=error.df, terms=terms,
repeated=FALSE, type="II", test=test)
}
else {
if (!is.null(imatrix)){
X.design <- do.call(cbind, imatrix)
ncols <- sapply(imatrix, ncol)
end <- cumsum(ncols)
start <- c(1, (end + 1))[-(length(end) + 1)]
cols <- mapply(seq, from=start, to=end)
iterms <- names(end)
names(cols) <- iterms
check.imatrix(X.design, iterms)
}
else {
if (is.null(idesign)) stop("idesign (intra-subject design) missing.")
for (i in 1:length(idata)){
if (is.null(attr(idata[,i], "contrasts"))){
contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2]
else icontrasts[1]
}
}
X.design <- model.matrix(idesign, data=idata)
iintercept <- has.intercept(X.design)
iterms <- term.names(idesign)
if (iintercept) iterms <- c("(Intercept)", iterms)
check.imatrix(X.design)
}
df <- rep(0, (n.terms + intercept)*length(iterms))
hnames <- rep("", length(df))
P <- SSPEH <- SSP <- as.list(df)
i <- 0
for (iterm in iterms){
if (intercept){
i <- i + 1
hyp.matrix.1 <- I.p[-1,,drop=FALSE]
SSP1 <- linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)$SSPH
lh2 <- linearHypothesis(mod, I.p, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)
SSP[[i]] <- lh2$SSPH - SSP1
SSPEH[[i]] <- lh2$SSPE
P[[i]] <- lh2$P
df[i] <- 1
hnames[i] <- iterm
}
for (term in 1:n.terms){
subs <- which(assign == term)
i <- i + 1
Test <- SSP.term(terms[term], iterm)
SSP[[i]] <- Test$SSP
SSPEH[[i]] <- Test$SSPE
P[[i]] <- Test$P
df[i]<- length(subs)
hnames[i] <- if (iterm == "(Intercept)") terms[term]
else paste(terms[term], ":", iterm, sep="")
}
}
names(df) <- names(P) <- names(SSP) <- names(SSPEH) <- hnames
result <- list(SSP=SSP, SSPE=SSPEH, P=P, df=df, error.df=error.df,
terms=hnames, repeated=TRUE, type="II", test=test,
idata=idata, idesign=idesign, icontrasts=icontrasts, imatrix=imatrix,
singular=Test$singular)
}
class(result) <- "Anova.mlm"
result
}
print.Anova.mlm <- function(x, ...){
if ((!is.null(x$singular)) && x$singular) stop("singular error SSP matrix; multivariate tests unavailable\ntry summary(object, multivariate=FALSE)")
test <- x$test
repeated <- x$repeated
ntests <- length(x$terms)
tests <- matrix(NA, ntests, 4)
if (!repeated) SSPE.qr <- qr(x$SSPE)
for (term in 1:ntests){
# some of the code here adapted from stats:::summary.manova
eigs <- Re(eigen(qr.coef(if (repeated) qr(x$SSPE[[term]]) else SSPE.qr,
x$SSP[[term]]), symmetric = FALSE)$values)
tests[term, 1:4] <- switch(test,
Pillai = Pillai(eigs, x$df[term], x$error.df),
Wilks = Wilks(eigs, x$df[term], x$error.df),
"Hotelling-Lawley" = HL(eigs, x$df[term], x$error.df),
Roy = Roy(eigs, x$df[term], x$error.df))
}
ok <- tests[, 2] >= 0 & tests[, 3] > 0 & tests[, 4] > 0
ok <- !is.na(ok) & ok
tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3], tests[ok, 4],
lower.tail = FALSE))
rownames(tests) <- x$terms
colnames(tests) <- c("Df", "test stat", "approx F", "num Df", "den Df", "Pr(>F)")
tests <- structure(as.data.frame(tests),
heading = paste("\nType ", x$type, if (repeated) " Repeated Measures",
" MANOVA Tests: ", test, " test statistic", sep=""),
class = c("anova", "data.frame"))
print(tests)
invisible(x)
}
# path <- "D:/R-package-sources/car/R"
# files <- list.files(path, pattern=".*\\.R")
# files <- paste(path, files, sep="/")
# for (file in files) source(file)
# summary.Anova.mlm and print.summary.Anova.mlm methods
# with contributions from Gabriel Baud-Bovy
summary.Anova.mlm <- function (object, test.statistic, univariate=object$repeated, multivariate=TRUE, p.adjust.method, ...) {
GG <- function(SSPE, P) { # Greenhouse-Geisser correction
p <- nrow(SSPE)
if (p < 2)
return(NA)
lambda <- eigen(SSPE %*% solve(t(P) %*% P))$values
lambda <- lambda[lambda > 0]
((sum(lambda)/p)^2)/(sum(lambda^2)/p)
}
HF <- function(gg, error.df, p) { # Huynh-Feldt correction
((error.df + 1) * p * gg - 2)/(p * (error.df - p * gg))
}
mauchly <- function(SSD, P, df) {
# most of this function borrowed from stats:::mauchly.test.SSD
if (nrow(SSD) < 2)
return(c(NA, NA))
Tr <- function(X) sum(diag(X))
p <- nrow(P)
I <- diag(p)
Psi <- t(P) %*% I %*% P
B <- SSD
pp <- nrow(SSD)
U <- solve(Psi, B)
n <- df
logW <- log(det(U)) - pp * log(Tr(U/pp))
rho <- 1 - (2 * pp^2 + pp + 2)/(6 * pp * n)
w2 <- (pp + 2) * (pp - 1) * (pp - 2) * (2 * pp^3 + 6 *
pp^2 + 3 * p + 2)/(288 * (n * pp * rho)^2)
z <- -n * rho * logW
f <- pp * (pp + 1)/2 - 1
Pr1 <- pchisq(z, f, lower.tail = FALSE)
Pr2 <- pchisq(z, f + 4, lower.tail = FALSE)
pval <- Pr1 + w2 * (Pr2 - Pr1)
c(statistic = c(W = exp(logW)), p.value = pval)
}
if (missing(test.statistic))
test.statistic <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy")
test.statistic <- match.arg(test.statistic, c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), several.ok = TRUE)
nterms <- length(object$terms)
summary.object <- list(type=object$type, repeated=object$repeated,
multivariate.tests=NULL, univariate.tests=NULL,
pval.adjustments=NULL, sphericity.tests=NULL)
if (multivariate){
summary.object$multivariate.tests <- vector(nterms, mode="list")
names(summary.object$multivariate.tests) <- object$terms
summary.object$SSPE <- object$SSPE
for (term in 1:nterms) {
hyp <- list(SSPH = object$SSP[[term]],
SSPE = if (object$repeated) object$SSPE[[term]] else object$SSPE,
P = if (object$repeated) object$P[[term]] else NULL,
test = test.statistic, df = object$df[term],
df.residual = object$error.df, title = object$terms[term])
class(hyp) <- "linearHypothesis.mlm"
summary.object$multivariate.tests[[term]] <- hyp
}
}
if (object$repeated && univariate) {
singular <- object$singular
error.df <- object$error.df
table <- matrix(0, nterms, 6)
table2 <- matrix(0, nterms, 4)
table3 <- matrix(0, nterms, 2)
rownames(table3) <- rownames(table2) <- rownames(table) <- object$terms
colnames(table) <- c("SS", "num Df", "Error SS", "den Df", "F", "Pr(>F)")
colnames(table2) <- c("GG eps", "Pr(>F[GG])", "HF eps","Pr(>F[HF])")
colnames(table3) <- c("Test statistic", "p-value")
if (singular)
warning("Singular error SSP matrix:\nnon-sphericity test and corrections not available")
for (term in 1:nterms) {
SSP <- object$SSP[[term]]
SSPE <- object$SSPE[[term]]
P <- object$P[[term]]
p <- ncol(P)
PtPinv <- solve(t(P) %*% P)
gg <- if (!singular) GG(SSPE, P) else NA
table[term, "SS"] <- sum(diag(SSP %*% PtPinv))
table[term, "Error SS"] <- sum(diag(SSPE %*% PtPinv))
table[term, "num Df"] <- object$df[term] * p
table[term, "den Df"] <- error.df * p
table[term, "F"] <- (table[term, "SS"]/table[term, "num Df"])/
(table[term, "Error SS"]/table[term, "den Df"])
table[term, "Pr(>F)"] <- pf(table[term, "F"], table[term, "num Df"], table[term, "den Df"],
lower.tail = FALSE)
table2[term, "GG eps"] <- gg
table2[term, "HF eps"] <- if (!singular) HF(gg, error.df, p) else NA
table3[term, ] <- if (!singular) mauchly(SSPE, P, object$error.df) else NA
}
table3 <- na.omit(table3)
if (nrow(table3) > 0) {
table2[, "Pr(>F[GG])"] <- pf(table[, "F"], table2[, "GG eps"] *
table[, "num Df"], table2[, "GG eps"] * table[, "den Df"],
lower.tail = FALSE)
table2[, "Pr(>F[HF])"] <- pf(table[, "F"], pmin(1, table2[, "HF eps"]) *
table[, "num Df"], pmin(1, table2[, "HF eps"]) * table[, "den Df"],
lower.tail = FALSE)
table2 <- na.omit(table2)
if (any(table2[, "HF eps"] > 1)) warning("HF eps > 1 treated as 1")
}
class(table3) <- class(table) <- "anova"
summary.object$univariate.tests <- table
summary.object$pval.adjustments <- table2
summary.object$sphericity.tests <- table3
}
if (!object$repeated && univariate) {
SS <- sapply(object$SSP, diag)
SSE <- diag(object$SSPE)
df <- object$df
dfe <- object$error.df
F <- (SS/df)/(SSE/dfe)
SS <- cbind(SS, residuals=SSE)
SS <- rbind(df=c(df, residuals=dfe), SS)
p <- pf(F, df, dfe, lower.tail=FALSE)
result <- list(SS=t(SS), F=t(F), p=t(p), type=object$type)
if (!missing(p.adjust.method)){
if (isTRUE(p.adjust.method)) p.adjust.method <- "holm"
p.adj <- apply(p, 2, p.adjust, method=p.adjust.method)
result$p.adjust <- t(p.adj)
result$p.adjust.method <- p.adjust.method
}
class(result) = "univaov"
summary.object$univaov <- result
}
class(summary.object) <- "summary.Anova.mlm"
summary.object
}
print.summary.Anova.mlm <- function(x, digits = getOption("digits"), SSP=TRUE, SSPE=SSP, ... ) {
if (!is.null(x$multivariate.tests)) {
cat(paste("\nType ", x$type, if (x$repeated)
" Repeated Measures", " MANOVA Tests:\n", sep = ""))
if ((!x$repeated) && SSPE) {
cat("\nSum of squares and products for error:\n")
print(x$SSPE, digits = digits)
}
for (term in 1:length(x$multivariate.tests)) {
cat(paste("\n------------------------------------------\n",
"\nTerm:", names(x$multivariate.tests)[term], "\n"))
print(x$multivariate.tests[[term]], digits = digits, SSP=SSP, SSPE=FALSE, ...)
}
}
if (!is.null(x$univariate.tests)) {
cat("\nUnivariate Type", x$type, "Repeated-Measures ANOVA Assuming Sphericity\n\n")
print(x$univariate.tests)
if (nrow(x$sphericity.tests) > 0) {
cat("\n\nMauchly Tests for Sphericity\n\n")
print(x$sphericity.tests)
cat("\n\nGreenhouse-Geisser and Huynh-Feldt Corrections\n",
"for Departure from Sphericity\n\n")
table <- x$pval.adjustments[, 1:2, drop = FALSE]
class(table) <- "anova"
print(table)
cat("\n")
table <- x$pval.adjustments[, 3:4, drop = FALSE]
class(table)
print(table)
}
}
if (!is.null(x$univaov)){
print(x$univaov, ...)
}
invisible(x)
}
print.univaov <- function(x, digits = max(getOption("digits") - 2L, 3L),
style=c("wide", "long"),
by=c("response", "term"),
...){
style <- match.arg(style)
if (style == "wide") {
cat("\n Type", x$type, "Sums of Squares\n")
print(x$SS, digits=digits)
cat("\n F-tests\n")
F <- x$F
print(round(F, 2))
cat("\n p-values\n")
p <- format.pval(x$p)
p <- matrix(p, nrow=nrow(F))
rownames(p) <- rownames(F)
colnames(p) <- colnames(F)
print(p, quote=FALSE)
if (!is.null(x$p.adjust)){
cat("\n p-values adjusted (by term) for simultaneous inference by", x$p.adjust.method, "method\n")
p.adjust <- format.pval(x$p.adjust)
p.adjust <- matrix(p.adjust, nrow=nrow(F))
rownames(p.adjust) <- rownames(F)
colnames(p.adjust) <- colnames(F)
print(p.adjust, quote=FALSE)
}
}
else {
x.df <- as.data.frame(x, by=by)
x.df$F <- round(x.df$F, 2)
x.df$p <- format.pval(x.df$p)
if (!is.null(x$p.adjust)) x.df$"adjusted p" <- format.pval(x.df$"adjusted p")
cat("\n Type", x$type, "Sums of Squares and F tests\n")
print(x.df, quote=FALSE, digits=digits)
}
invisible(x)
}
as.data.frame.univaov <- function(x, row.names, optional, by=c("response", "term"), ...) {
melt <- function(data, varnames = names(dimnames(data)), value.name = "value") {
dn <- dimnames(data)
labels <- expand.grid( dn[[1]], dn[[2]])
colnames(labels) <- varnames
value_df <- setNames(data.frame(as.vector(data)), value.name)
cbind(labels, value_df)
}
nv <- ncol(x$F)
nt <- nrow(x$F)
by <- match.arg(by)
if (by=="response") {
vn <- c("term", "response")
df <- matrix(x$SS[1:nt, "df", drop=FALSE], nrow=nt, ncol=nv)
SS <- melt(x$SS[1:nt, -1, drop=FALSE], varnames=vn, value.name="SS")
F <- melt(x$F, varnames=vn, value.name="F")
p <- melt(x$p, varnames=vn, value.name="p")
if (!is.null(x$p.adjust)) p.adjust <- melt(x$p.adjust, varnames=vn, value.name="adjusted p")
}
else {
vn <- rev(c("term", "response"))
df <- t(matrix(x$SS[1:nt, "df", drop=FALSE], nrow=nt, ncol=nv))
SS <- melt(t(x$SS[1:nt, -1, drop=FALSE]), varnames=vn, value.name="SS")
F <- melt(t(x$F), varnames=vn, value.name="F")
p <- melt(t(x$p), varnames=vn, value.name="p")
if (!is.null(x$p.adjust)) p.adjust <- melt(x$p.adjust, varnames=vn, value.name="adjusted p")
}
result <- cbind(SS[,c(2,1,3)], df=c(df), F=F[,"F"], p=p[,"p"])
if (!is.null(x$p.adjust)) result <- cbind(result, "adjusted p"=p.adjust[, "adjusted p"])
result
}
Anova.manova <- function(mod, ...){
class(mod) <- c("mlm", "lm")
Anova(mod, ...)
}
Manova <- function(mod, ...){
UseMethod("Manova")
}
Manova.mlm <- function(mod, ...){
Anova(mod, ...)
}
# Cox regression models
df.residual.coxph <- function(object, ...){
object$n - sum(!is.na(coef(object)))
}
alias.coxph <- function(model){
if(any(which <- is.na(coef(model)))) return(list(Complete=which))
else list()
}
logLik.coxph <- function(object, ...) object$loglik[2]
Anova.coxph <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (length((mod$rscore) > 0) && (test.statistic == "LR")){
warning("LR tests unavailable with robust variances\nWald tests substituted")
test.statistic <- "Wald"
}
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.coxph(mod),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", vcov.=vcov(mod))),
III=switch(test.statistic,
LR=Anova.III.LR.coxph(mod),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", vcov.=vcov(mod))),
"2"=switch(test.statistic,
LR=Anova.II.LR.coxph(mod),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", vcov.=vcov(mod))),
"3"=switch(test.statistic,
LR=Anova.III.LR.coxph(mod),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", vcov.=vcov(mod))))
}
Anova.II.LR.coxph <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
method <- mod$method
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- survival::coxph(mod$y ~ X[, -exclude.1, drop = FALSE], method=method)
loglik.1 <- logLik(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
survival::coxph(mod$y ~ X[, -exclude.2, drop = FALSE], method=method)
}
loglik.2 <- logLik(mod.2)
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
Anova.III.LR.coxph <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
method <- mod$method
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
df <- df.terms(mod)
LR <- p <- rep(0, n.terms)
loglik1 <- logLik(mod)
for (term in 1:n.terms){
mod.0 <- survival::coxph(mod$y ~ X[, -which.nms(names[term])], method=method)
LR[term] <- -2*(logLik(mod.0) - loglik1)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df","Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result,"heading") <- "Analysis of Deviance Table (Type III tests)"
result
}
# parametric survival regression models
alias.survreg <- function(model){
if(any(which <- diag(vcov(model)) < 1e-10)) return(list(Complete=which))
else list()
}
logLik.survreg <- function(object, ...) object$loglik[2]
Anova.survreg <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (length((mod$rscore) > 0) && (test.statistic == "LR")){
warning("LR tests unavailable with robust variances\nWald tests substituted")
test.statistic <- "Wald"
}
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.survreg(mod),
Wald=Anova.II.Wald.survreg(mod)),
III=switch(test.statistic,
LR=Anova.III.LR.survreg(mod),
Wald=Anova.III.Wald.survreg(mod)),
"2"=switch(test.statistic,
LR=Anova.II.LR.survreg(mod),
Wald=Anova.II.Wald.survreg(mod)),
"3"=switch(test.statistic,
LR=Anova.III.LR.survreg(mod),
Wald=Anova.III.Wald.survreg(mod)))
}
Anova.II.LR.survreg <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
dist <- mod$dist
scale <- mod$call$scale
weights <- model.frame(mod)$"(weights)"
arg.list <- list(dist=dist)
if (!is.null(scale)) arg.list$scale <- scale
if (!is.null(weights)) arg.list$weights <- weights
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
asgn <- asgn[asgn != 0]
if (has.intercept(mod)){
int <- which(names == "(Intercept)")
X <- X[, -int]
names <- names[-int]
}
n.terms <- length(names)
if (n.terms < 2) return(anova(mod))
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
y <- model.frame(mod)[,1]
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
arg.list$formula <- y ~ X[, -exclude.1, drop = FALSE]
mod.1 <- do.call(survival::survreg, arg.list)
# mod.1 <- survival::survreg(y ~ X[, -exclude.1, drop = FALSE])
loglik.1 <- logLik(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
arg.list$formula <- y ~ X[, -exclude.2, drop = FALSE]
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
do.call(survival::survreg, arg.list)
# survival::survreg(y ~ X[, -exclude.2, drop = FALSE])
}
loglik.2 <- logLik(mod.2)
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
Anova.III.LR.survreg <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
dist <- mod$dist
scale <- mod$call$scale
weights <- model.frame(mod)$"(weights)"
arg.list <- list(dist=dist)
if (!is.null(scale)) arg.list$scale <- scale
if (!is.null(weights)) arg.list$weights <- weights
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
asgn <- asgn[asgn != 0]
if (has.intercept(mod)){
int <- which(names == "(Intercept)")
X <- X[, -int]
names <- names[-int]
}
n.terms <- length(names)
if (n.terms < 2) return(anova(mod))
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
y <- model.frame(mod)[,1]
loglik1 <- logLik(mod)
for (term in 1:n.terms){
arg.list$formula <- y ~ X[, -which.nms(names[term])]
mod.0 <- do.call(survival::survreg, arg.list)
# mod.0 <- survival::survreg(y ~ X[, -which.nms(names[term])])
LR[term] <- -2*(logLik(mod.0) - loglik1)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df","Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result,"heading") <- "Analysis of Deviance Table (Type III tests)"
result
}
Anova.II.Wald.survreg <- function(mod){
V <- vcov(mod)
p <- which(rownames(V) == "Log(scale)")
if (length(p) > 0) V <- V[-p, -p]
Anova.II.default(mod, V, test="Chisq")
}
Anova.III.Wald.survreg <- function(mod){
V <- vcov(mod)
p <- which(rownames(V) == "Log(scale)")
if (length(p) > 0) V <- V[-p, -p]
Anova.III.default(mod, V, test="Chisq")
}
# Default Anova() method: requires methods for vcov() (if vcov. argument not specified) and coef().
Anova.default <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
switch(type,
II=Anova.II.default(mod, vcov., test.statistic, singular.ok=singular.ok),
III=Anova.III.default(mod, vcov., test.statistic, singular.ok=singular.ok),
"2"=Anova.II.default(mod, vcov., test.statistic, singular.ok=singular.ok),
"3"=Anova.III.default(mod, vcov., test.statistic, singular.ok=singular.ok))
}
assignVector <- function(model){
m <- model.matrix(model)
assign <- attr(m, "assign")
if (!is.null(assign)) return (assign)
m <- model.matrix(formula(model), data=model.frame(model))
assign <- attr(m, "assign")
if (!has.intercept(model)) assign <- assign[assign != 0]
assign
}
Anova.II.default <- function(mod, vcov., test, singular.ok=TRUE, ...){
hyp.term <- function(term){
which.term <- which(term==names)
subs.term <- if (is.list(assign)) assign[[which.term]] else which(assign == which.term)
relatives <- relatives(term, names, fac)
subs.relatives <- NULL
for (relative in relatives){
sr <- if (is.list(assign)) assign[[relative]] else which(assign == relative)
subs.relatives <- c(subs.relatives, sr)
}
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.default(mod, hyp.matrix.term,
vcov.=vcov., test=test, singular.ok=singular.ok, ...)
if (test=="Chisq") c(statistic=hyp$Chisq[2], df=hyp$Df[2])
else c(statistic=hyp$F[2], df=hyp$Df[2])
}
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
assign <- assignVector(mod) # attr(model.matrix(mod), "assign")
if (!is.list(assign)) assign[!not.aliased] <- NA
else if (intercept) assign <- assign[-1]
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
df <- c(rep(0, n.terms), df.residual(mod))
if (inherits(mod, "coxph")){
assign <- assign[assign != 0]
clusters <- grep("^cluster\\(", names)
if (length(clusters) > 0) {
names <- names[-clusters]
df <- df[-clusters]
n.terms <- n.terms - length(clusters)
}
}
# if (inherits(mod, "plm")) assign <- assign[assign != 0]
p <- teststat <- rep(0, n.terms + 1)
teststat[n.terms + 1] <- p[n.terms + 1] <- NA
for (i in 1:n.terms){
hyp <- hyp.term(names[i])
teststat[i] <- abs(hyp["statistic"])
df[i] <- abs(hyp["df"])
p[i] <- if (test == "Chisq")
pchisq(teststat[i], df[i], lower.tail=FALSE)
else pf(teststat[i], df[i], df[n.terms + 1], lower.tail=FALSE)
}
result <- if (test == "Chisq"){
if (length(df) == n.terms + 1) df <- df[1:n.terms]
data.frame(df, teststat[!is.na(teststat)], p[!is.na(teststat)])
}
else data.frame(df, teststat, p)
if (nrow(result) == length(names) + 1) names <- c(names,"Residuals")
row.names(result) <- names
names(result) <- c ("Df", test, if (test == "Chisq") "Pr(>Chisq)"
else "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.default <- function(mod, vcov., test, singular.ok=FALSE, ...){
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- assignVector(mod) # attr(model.matrix(mod), "assign")
df <- c(rep(0, n.terms), df.residual(mod))
if (inherits(mod, "coxph")){
if (intercept) names <- names[-1]
assign <- assign[assign != 0]
clusters <- grep("^cluster\\(", names)
if (length(clusters) > 0) {
names <- names[-clusters]
df <- df[-clusters]
n.terms <- n.terms - length(clusters)
}
}
# if (inherits(mod, "plm")) assign <- assign[assign != 0]
if (intercept) df[1] <- sum(grepl("^\\(Intercept\\)", names(coef(mod))))
teststat <- rep(0, n.terms + 1)
p <- rep(0, n.terms + 1)
teststat[n.terms + 1] <- p[n.terms + 1] <- NA
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1:n.terms){
subs <- if (is.list(assign)) assign[[term]] else 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.default(mod, hyp.matrix,
vcov.=vcov., test=test, singular.ok=singular.ok, ...)
teststat[term] <- if (test=="Chisq") hyp$Chisq[2] else hyp$F[2]
df[term] <- abs(hyp$Df[2])
p[term] <- if (test == "Chisq")
pchisq(teststat[term], df[term], lower.tail=FALSE)
else pf(teststat[term], df[term], df[n.terms + 1], lower.tail=FALSE)
}
}
result <- if (test == "Chisq"){
if (length(df) == n.terms + 1) df <- df[1:n.terms]
data.frame(df, teststat[!is.na(teststat)], p[!is.na(teststat)])
}
else data.frame(df, teststat, p)
if (nrow(result) == length(names) + 1) names <- c(names,"Residuals")
row.names(result) <- names
names(result) <- c ("Df", test, if (test == "Chisq") "Pr(>Chisq)"
else "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
## functions for mixed models
# the following function, not exported, to make car consistent with CRAN and development versions of lme4 and with nlme
fixef <- function (object){
if (isS4(object)) {
if (!inherits(object, "merMod")) object@fixef
else lme4::fixef(object)
}
else object$coefficients$fixed
}
Anova.merMod <- function(mod, type=c("II","III", 2, 3),
test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
Anova.mer(mod=mod, type=type, test.statistic=test.statistic, vcov.=vcov.,
singular.ok=singular.ok, ...)
}
Anova.mer <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
switch(type,
II=Anova.II.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
III=Anova.III.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
"2"=Anova.II.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
"3"=Anova.III.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok))
}
Anova.II.mer <- function(mod, vcov., singular.ok=TRUE, test=c("Chisq", "F"), ...){
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(mod, hyp.matrix.term,
vcov.=vcov., singular.ok=singular.ok, test=test, ...)
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]))
}
test <- match.arg(test)
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
if (!missing(vcov.)){
vcov. <- if (test == "F"){
# if (!require("pbkrtest")) stop("pbkrtest package required for F-tests on linear mixed model")
as.matrix(pbkrtest::vcovAdj(mod, details=0))
}
else vcov(mod)
}
assign <- attr(model.matrix(mod), "assign")
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
p <- teststat <- df <- res.df <- rep(0, n.terms)
for (i in 1: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] <- if (test == "Chisq") pchisq(teststat[i], df[i], lower.tail=FALSE)
else pf(teststat[i], df[i], res.df[i], lower.tail=FALSE)
}
if (test=="Chisq"){
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(mod)))
}
else {
result <- data.frame(teststat, df, res.df, p)
row.names(result) <- names
names(result) <- c ("F", "Df", "Df.res", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)\n",
paste("Response:", responseName(mod)))
}
result
}
Anova.III.mer <- function(mod, vcov., singular.ok=FALSE, test=c("Chisq", "F"), ...){
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- attr(model.matrix(mod), "assign")
p <- teststat <- df <- res.df <- rep(0, n.terms)
if (intercept) df[1] <- 1
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
if (!missing(vcov.)){
vcov. <- if (test == "F"){
# if (!require("pbkrtest")) stop("pbkrtest package required for F-tests on linear mixed model")
as.matrix(pbkrtest::vcovAdj(mod, details=0))
}
else vcov(mod)
}
for (term in 1: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(mod, hyp.matrix, test=test,
vcov.=vcov., singular.ok=singular.ok, ...)
if (test == "Chisq"){
teststat[term] <- hyp$Chisq[2]
df[term] <- abs(hyp$Df[2])
p[term] <- pchisq(teststat[term], df[term], lower.tail=FALSE)
}
else{
teststat[term] <- hyp$F[2]
df[term] <- abs(hyp$Df[2])
res.df[term]=hyp$Res.Df[2]
p[term] <- pf(teststat[term], df[term], res.df[term], lower.tail=FALSE)
}
}
}
if (test == "Chisq"){
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(mod)))
}
else {
result <- data.frame(teststat, df, res.df, p)
row.names(result) <- names
names(result) <- c ("F", "Df", "Df.res", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)\n",
paste("Response:", responseName(mod)))
}
result
}
Anova.lme <- function(mod, type=c("II","III", 2, 3),
vcov.=vcov(mod), singular.ok, ...){
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"
switch(type,
II=Anova.II.lme(mod, vcov., singular.ok=singular.ok),
III=Anova.III.lme(mod, vcov., singular.ok=singular.ok),
"2"=Anova.II.lme(mod, vcov., singular.ok=singular.ok),
"3"=Anova.III.lme(mod, vcov., singular.ok=singular.ok))
}
Anova.II.lme <- function(mod, vcov., singular.ok=TRUE, ...){
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(mod, hyp.matrix.term,
vcov.=vcov., singular.ok=singular.ok, ...)
c(statistic=hyp$Chisq[2], df=hyp$Df[2])
}
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
assign <- attr(model.matrix(mod), "assign")
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
p <- teststat <- df <- rep(0, n.terms)
for (i in 1:n.terms){
hyp <- hyp.term(names[i])
teststat[i] <- abs(hyp["statistic"])
df[i] <- abs(hyp["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 tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.lme <- function(mod, vcov., singular.ok=FALSE, ...){
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- attr(model.matrix(mod), "assign")
df <- rep(0, n.terms)
if (intercept) df[1] <- 1
p <- teststat <-rep(0, n.terms)
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1: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(mod, hyp.matrix,
vcov.=vcov., singular.ok=singular.ok, ...)
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 tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.svyglm <- function(mod, ...) Anova.default(mod, ...)
Anova.rlm <- function(mod, ...) Anova.default(mod, test.statistic="F", ...)
Anova.coxme <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Wald", "LR"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.coxme(mod, ...),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", ...)),
III=switch(test.statistic,
LR=stop("type-III LR tests not available for coxme models"),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", ...)),
"2"=switch(test.statistic,
LR=Anova.II.LR.coxme(mod, ...),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", ...)),
"3"=switch(test.statistic,
LR=stop("type-III LR tests not available for coxme models"),
Wald=Anova.default(mod, type="III", test.statistic="Chisq")))
}
Anova.II.LR.coxme <- function(mod, ...){
if (!requireNamespace("coxme")) stop("coxme package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
random <- mod$formulaList$random
random <- sapply(random, as.character)[2, ]
random <- paste(paste0("(", random, ")"), collapse=" + ")
fixed <- as.character(mod$formulaList$fixed)[3]
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
formula <- paste0(". ~ . - ", paste(c(names[term], rels), collapse=" - "), " + ", random)
mod.1 <- update(mod, as.formula(formula))
loglik.1 <- logLik(mod.1, type="integrated")
mod.2 <- if (length(rels) == 0) mod
else {
formula <- paste0(". ~ . - ", paste(rels, collapse=" - "), " + ", random)
update(mod, as.formula(formula))
}
loglik.2 <- logLik(mod.2, type="integrated")
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
jonathon-love/car3 documentation built on May 19, 2019, 7:28 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#-------------------------------------------------------------------------------
# Revision history:
# 2009-01-05: bug fix in Anova.II.lm(). J. Fox
# 2009-01-16: Cox models with clusters now handled. J. Fox
# 2009-09-16: reworked glm and lm methods to handle aliased parameters. J. Fox
# 2009-09-30: renamed "Anova" to "Analysis of Deviance" in output for some methods. J. Fox
# 2009-12-22: modified Anova.mlm() to handle a user-supplied within-subject model matrix. J. Fox
# 2009-12-28: named the components of P in Anova.III.mlm(). John
# 2010-01-01: Anova.II.mlm() now hands off (again) to Anova.III.mlm() when there
# is only an intercept in the between-subjects model
# 2010-02-17: Fixed bug that caused some models with aliased coefficients to fail. J. Fox
# 2010-06-14: added wcrossprod and allow use of observation weights in Anova.mlm()
# 2010-06-28: Fixed Anova() tables for coxph and survreg models
# (failed because of changes in survival package.
# 2011-01-21: Added functions for mixed models. J. Fox
# 2011-01-25: Fixed Anova.polr() and Anova.multinom() to work with models with only one term. J. Fox
# 2011-05-19: local fixef() to avoid nlme/lme4 issues. J. Fox
# 2011-05-11: changed order of columns in ANOVA tables for mixed models. J. Fox
# 2011-11-27: added Anova.svyglm(). J. Fox
# 2011-12-31: fixed bug in Anova.II(and III).F.glm() when na.exclude used. J. Fox
# 2012-02-28: added test.statistic argument to Anova.mer(). J.Fox
# 2012-03-02: fixed test abbreviation of test.statistic argument to Anova.default()
# called by other Anova() methods. J. Fox
# 2013-06-17: modified summary.Anova.mlm(), introduced print.summary.Anova.mlm(),
# adapting code contributed by Gabriel Baud-Bovy. J. Fox
# 2013-06-20: added Anova.merMod() method. J. Fox
# 2013-06-22: tweaks to local fixef(). J. Fox
# 2013-06-22: test argument uniformly uses "Chisq" rather than "chisq". J. Fox
# 2013-08-19: replaced calls to print.anova(). J. Fox
# 2014-08-17: added calls to requireNamespace() and :: where needed (doesn't work for pbkrtest). J. Fox
# 2014-08-18: fixed bugs in Anova.survreg() for types II, III LR tests and Wald tests. J. Fox
# 2014-09-23: added Anova.rlm(). J. Fox
# 2014-10-10: removed MASS:: from calls to polr(). John
# 2014-12-18: check that residual df and SS are nonzero in Anova.lm(). John
# 2015-01-27: vcovAdj() and methods now imported from pbkrtest. John
# 2015-02-18: force evaluation of vcov. when it's a function. John
# 2015-04-30: don't allow error.estimate="dispersion" for F-tests in binomial
# and Poission GLMs. John
# 2015-08-29: fixed Anova() for coxph models with clusters. John
# 2015-09-04: added support for coxme models. John
# 2015-09-11: modified Anova.default() to work with vglm objects from VGAM. John
# 2015-09-15: fixed Anova.default() so that F-tests work again. John
# 2015-11-13: modify Anova.coxph() to take account of method/ties argument. John
# 2016-06-03: added SSP and SSPE args to print.summary.Anova.mlm(). John
# 2016-06-25: added code to optionally print univariate ANOVAs for a mlm. John
#-------------------------------------------------------------------------------
# Type II and III tests for linear, generalized linear, and other models (J. Fox)
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))]
}
Anova <- function(mod, ...){
UseMethod("Anova", mod)
}
# linear models
Anova.lm <- function(mod, error, type=c("II","III", 2, 3),
white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4"),
vcov.=NULL, singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
if (df.residual(mod) == 0) stop("residual df = 0")
if (deviance(mod) < sqrt(.Machine$double.eps)) stop("residual sum of squares is 0 (within rounding error)")
type <- as.character(type)
white.adjust <- as.character(white.adjust)
type <- match.arg(type)
white.adjust <- match.arg(white.adjust)
if (missing(singular.ok)){
singular.ok <- type == "2" || type == "II"
}
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
if (white.adjust != "FALSE"){
if (white.adjust == "TRUE") white.adjust <- "hc3"
return(Anova.default(mod, type=type, vcov.=hccm(mod, type=white.adjust), test.statistic="F",
singular.ok=singular.ok, ...))
}
else if (!is.null(vcov.)) return(Anova.default(mod, type=type, vcov.=vcov., test.statistic="F",
singular.ok=singular.ok, ...))
switch(type,
II=Anova.II.lm(mod, error, singular.ok=singular.ok, ...),
III=Anova.III.lm(mod, error, singular.ok=singular.ok, ...),
"2"=Anova.II.lm(mod, error, singular.ok=singular.ok, ...),
"3"=Anova.III.lm(mod, error, singular.ok=singular.ok,...))
}
Anova.aov <- function(mod, ...){
class(mod) <- "lm"
Anova.lm(mod, ...)
}
Anova.II.lm <- function(mod, error, singular.ok=TRUE, ...){
if (!missing(error)){
sumry <- summary(error, corr=FALSE)
s2 <- sumry$sigma^2
error.df <- error$df.residual
error.SS <- s2*error.df
}
SS.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(mod)))
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(SS=NA, df=0))
lh <- linearHypothesis(mod, hyp.matrix.term,
singular.ok=singular.ok, ...)
abs(c(SS=lh$"Sum of Sq"[2], df=lh$Df[2]))
}
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
I.p <- diag(length(coefficients(mod)))
assign <- mod$assign
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <-names[-1]
n.terms <- length(names)
p <- df <- f <- SS <- rep(0, n.terms + 1)
sumry <- summary(mod, corr = FALSE)
SS[n.terms + 1] <- if (missing(error)) sumry$sigma^2*mod$df.residual
else error.SS
df[n.terms + 1] <- if (missing(error)) mod$df.residual else error.df
p[n.terms + 1] <- f[n.terms + 1] <- NA
for (i in 1:n.terms){
ss <- SS.term(names[i])
SS[i] <- ss["SS"]
df[i] <- ss["df"]
f[i] <- df[n.terms+1]*SS[i]/(df[i]*SS[n.terms + 1])
p[i] <- pf(f[i], df[i], df[n.terms + 1], lower.tail = FALSE)
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(names,"Residuals")
names(result) <- c("Sum Sq", "Df", "F value", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Anova Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
# type III
Anova.III.lm <- function(mod, error, singular.ok=FALSE, ...){
if (!missing(error)){
error.df <- df.residual(error)
error.SS <- deviance(error)
}
else {
error.df <- df.residual(mod)
error.SS <- deviance(mod)
}
intercept <- has.intercept(mod)
I.p <- diag(length(coefficients(mod)))
Source <- term.names(mod)
n.terms <- length(Source)
p <- df <- f <- SS <- rep(0, n.terms + 1)
assign <- mod$assign
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1: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){
SS[term] <- NA
df[term] <- 0
f[term] <- NA
p[term] <- NA
}
else {
test <- if (missing(error)) linearHypothesis(mod, hyp.matrix,
singular.ok=singular.ok, ...)
else linearHypothesis(mod, hyp.matrix, error.SS=error.SS, error.df=error.df,
singular.ok=singular.ok, ...)
SS[term] <- test$"Sum of Sq"[2]
df[term] <- test$"Df"[2]
f[term] <- test$"F"[2]
p[term] <- test$"Pr(>F)"[2]
}
}
Source[n.terms + 1] <- "Residuals"
SS[n.terms + 1] <- error.SS
df[n.terms + 1] <- error.df
p[n.terms + 1] <- f[n.terms + 1] <- NA
result <- data.frame(SS, df, f, p)
row.names(result) <- Source
names(result) <- c("Sum Sq", "Df", "F value", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Anova Table (Type III tests)\n", paste("Response:", responseName(mod)))
result
}
# generalized linear models
Anova.glm <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald", "F"),
error, error.estimate=c("pearson", "dispersion", "deviance"), singular.ok, ...){
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
if (missing(singular.ok)){
singular.ok <- type == "2" || type == "II"
}
test.statistic <- match.arg(test.statistic)
error.estimate <- match.arg(error.estimate)
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="II", singular.ok=singular.ok),
F=Anova.II.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
III=switch(test.statistic,
LR=Anova.III.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="III", singular.ok=singular.ok),
F=Anova.III.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
"2"=switch(test.statistic,
LR=Anova.II.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="II", singular.ok=singular.ok),
F=Anova.II.F.glm(mod, error, error.estimate, singular.ok=singular.ok)),
"3"=switch(test.statistic,
LR=Anova.III.LR.glm(mod, singular.ok=singular.ok),
Wald=Anova.default(mod, type="III", singular.ok=singular.ok),
F=Anova.III.F.glm(mod, error, error.estimate, singular.ok=singular.ok)))
}
# type III
# LR test
Anova.III.LR.glm <- function(mod, singular.ok=FALSE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
Source <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(Source)
p <- df <- LR <- rep(0, n.terms)
dispersion <- summary(mod, corr = FALSE)$dispersion
deviance <- deviance(mod)/dispersion
for (term in 1:n.terms){
mod.1 <- drop1(mod, scope=eval(parse(text=paste("~",Source[term]))))
LR[term] <- (mod.1$Deviance[2]/dispersion)-deviance
df[term] <- mod.1$Df[2]
p[term] <- pchisq(LR[term], df[term], lower.tail = FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- Source
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova","data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n", paste("Response:", responseName(mod)))
result
}
# F test
Anova.III.F.glm <- function(mod, error, error.estimate, singular.ok=FALSE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
fam <- family(mod)$family
if ((fam == "binomial" || fam == "poisson") && error.estimate == "dispersion"){
warning("dispersion parameter estimated from the Pearson residuals, not taken as 1")
error.estimate <- "pearson"
}
if (missing(error)) error <- mod
df.res <- df.residual(error)
error.SS <- switch(error.estimate,
pearson=sum(residuals(error, "pearson")^2, na.rm=TRUE),
dispersion=df.res*summary(error, corr = FALSE)$dispersion,
deviance=deviance(error))
Source <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(Source)
p <- df <- f <- SS <-rep(0, n.terms+1)
f[n.terms+1] <- p[n.terms+1] <- NA
df[n.terms+1] <- df.res
SS[n.terms+1] <- error.SS
dispersion <- error.SS/df.res
deviance <- deviance(mod)
for (term in 1:n.terms){
mod.1 <- drop1(mod, scope=eval(parse(text=paste("~",Source[term]))))
df[term] <- mod.1$Df[2]
SS[term] <- mod.1$Deviance[2] - deviance
f[term] <- (SS[term]/df[term])/dispersion
p[term] <- pf(f[term], df[term], df.res, lower.tail = FALSE)
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(Source, "Residuals")
names(result) <- c("SS", "Df", "F", "Pr(>F)")
class(result) <- c("anova","data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)),
paste("Error estimate based on",
switch(error.estimate,
pearson="Pearson residuals", dispersion="estimated dispersion",
deviance="deviance"), "\n"))
result
}
# type II
# LR test
Anova.II.LR.glm <- function(mod, singular.ok=TRUE, ...){
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
# (some code adapted from drop1.glm)
which.nms <- function(name) which(asgn == which(names == name))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- mod$y
if (is.null(y)) y <- model.response(model.frame(mod), "numeric")
wt <- mod$prior.weights
if (is.null(wt)) wt <- rep(1, length(y))
asgn <- attr(X, 'assign')
df <- p <- LR <- rep(0, n.terms)
dispersion <- summary(mod, corr = FALSE)$dispersion
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- glm.fit(X[, -exclude.1, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
glm.fit(X[, -exclude.2, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
}
dev.2 <- deviance(mod.2)
df[term] <- df.residual(mod.1) - df.residual(mod.2)
if (df[term] == 0) LR[term] <- p[term] <- NA
else {
LR[term] <- (dev.1 - dev.2)/dispersion
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <-
c("Analysis of Deviance Table (Type II tests)\n", paste("Response:", responseName(mod)))
result
}
# F test
Anova.II.F.glm <- function(mod, error, error.estimate, singular.ok=TRUE, ...){
# (some code adapted from drop1.glm)
if (!singular.ok && any(is.na(coef(mod))))
stop("there are aliased coefficients in the model")
fam <- family(mod)$family
if ((fam == "binomial" || fam == "poisson") && error.estimate == "dispersion"){
warning("dispersion parameter estimated from the Pearson residuals, not taken as 1")
error.estimate <- "pearson"
}
which.nms <- function(name) which(asgn == which(names == name))
if (missing(error)) error <- mod
df.res <- df.residual(error)
error.SS <- switch(error.estimate,
pearson = sum(residuals(error, "pearson")^2, na.rm=TRUE),
dispersion = df.res*summary(error, corr = FALSE)$dispersion,
deviance = deviance(error))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- mod$y
if (is.null(y)) y <- model.response(model.frame(mod), "numeric")
wt <- mod$prior.weights
if (is.null(wt)) wt <- rep(1, length(y))
asgn <- attr(X, 'assign')
p <- df <- f <- SS <- rep(0, n.terms+1)
f[n.terms+1] <- p[n.terms+1] <- NA
df[n.terms+1] <- df.res
SS[n.terms+1] <- error.SS
dispersion <- error.SS/df.res
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- glm.fit(X[, -exclude.1, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
glm.fit(X[, -exclude.2, drop = FALSE], y, wt, offset = mod$offset,
family = mod$family, control = mod$control)
}
dev.2 <- deviance(mod.2)
df[term] <- df.residual(mod.1) - df.residual(mod.2)
if (df[term] == 0) SS[term] <- f[term] <- p[term] <- NA
else {
SS[term] <- dev.1 - dev.2
f[term] <- SS[term]/(dispersion*df[term])
p[term] <- pf(f[term], df[term], df.res, lower.tail=FALSE)
}
}
result <- data.frame(SS, df, f, p)
row.names(result) <- c(names, "Residuals")
names(result) <- c("SS", "Df", "F", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)),
paste("Error estimate based on",
switch(error.estimate,
pearson="Pearson residuals",
dispersion="estimated dispersion",
deviance="deviance"), "\n"))
result
}
# multinomial logit models (via multinom in the nnet package)
Anova.multinom <-
function (mod, type = c("II", "III", 2, 3), ...)
{
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
switch(type,
II = Anova.II.multinom(mod, ...),
III = Anova.III.multinom(mod, ...),
"2" = Anova.II.multinom(mod, ...),
"3" = Anova.III.multinom(mod, ...))
}
Anova.II.multinom <- function (mod, ...)
{
which.nms <- function(name) which(asgn == which(names ==
name))
fac <- attr(terms(mod), "factors")
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- mod$weights
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms) {
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels),
which.nms)))
mod.1 <-if (n.terms > 1) multinom(y ~ X[, -c(1, exclude.1)], weights=wt, trace=FALSE)
else multinom(y ~ 1, weights=wt, race=FALSE)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0)
mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
multinom(y ~ X[, -c(1, exclude.2)], weights=wt, trace=FALSE)
}
dev.2 <- deviance(mod.2)
LR[term] <- dev.1 - dev.2
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.multinom <- function (mod, ...)
{
names <- if (has.intercept(mod)) term.names(mod)[-1]
else term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- mod$weights
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
deviance <- deviance(mod)
for (term in 1:n.terms) {
mod.1 <- if (n.terms > 1) multinom(y ~ X[, term != asgn][, -1], weights=wt, trace=FALSE)
else multinom(y ~ 1, weights=wt, trace=FALSE)
LR[term] <- deviance(mod.1) - deviance
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
# proportional-odds logit models (via polr in the MASS package)
Anova.polr <- function (mod, type = c("II", "III", 2, 3), ...)
{
type <- as.character(type)
type <- match.arg(type)
if (has.intercept(mod) && length(coef(mod)) == 1
&& (type == "2" || type == "II")) {
type <- "III"
warning("the model contains only an intercept: Type III test substituted")
}
switch(type,
II = Anova.II.polr(mod, ...),
III = Anova.III.polr(mod, ...),
"2" = Anova.II.polr(mod, ...),
"3" = Anova.III.polr(mod, ...))
}
Anova.II.polr <- function (mod, ...)
{
if (!requireNamespace("MASS")) stop("MASS package is missing")
which.nms <- function(name) which(asgn == which(names ==
name))
fac <- attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- model.weights(model.frame(mod))
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms) {
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels),
which.nms)))
mod.1 <- if (n.terms > 1) polr(y ~ X[, -c(1, exclude.1)], weights=wt)
else polr(y ~ 1, weights=wt)
dev.1 <- deviance(mod.1)
mod.2 <- if (length(rels) == 0)
mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
polr(y ~ X[, -c(1, exclude.2)], weights=wt)
}
dev.2 <- deviance(mod.2)
LR[term] <- dev.1 - dev.2
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.polr <- function (mod, ...)
{
if (!requireNamespace("MASS")) stop("MASS package is missing")
names <- term.names(mod)
n.terms <- length(names)
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
wt <- model.weights(model.frame(mod))
asgn <- attr(X, "assign")
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
deviance <- deviance(mod)
for (term in 1:n.terms) {
mod.1 <- if (n.terms > 1) polr(y ~ X[, term != asgn][, -1], weights=wt)
else polr(y ~ 1, weights=wt)
LR[term] <- deviance(mod.1) - deviance
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
# multivariate linear models
# the following 3 functions copied from the stats package (not exported from stats)
Pillai <- function (eig, q, df.res) {
test <- sum(eig/(1 + eig))
p <- length(eig)
s <- min(p, q)
n <- 0.5 * (df.res - p - 1)
m <- 0.5 * (abs(p - q) - 1)
tmp1 <- 2 * m + s + 1
tmp2 <- 2 * n + s + 1
c(test, (tmp2/tmp1 * test)/(s - test), s * tmp1, s * tmp2)
}
Wilks <- function (eig, q, df.res) {
test <- prod(1/(1 + eig))
p <- length(eig)
tmp1 <- df.res - 0.5 * (p - q + 1)
tmp2 <- (p * q - 2)/4
tmp3 <- p^2 + q^2 - 5
tmp3 <- if (tmp3 > 0)
sqrt(((p * q)^2 - 4)/tmp3)
else 1
c(test, ((test^(-1/tmp3) - 1) * (tmp1 * tmp3 - 2 * tmp2))/p/q,
p * q, tmp1 * tmp3 - 2 * tmp2)
}
HL <- function (eig, q, df.res) {
test <- sum(eig)
p <- length(eig)
m <- 0.5 * (abs(p - q) - 1)
n <- 0.5 * (df.res - p - 1)
s <- min(p, q)
tmp1 <- 2 * m + s + 1
tmp2 <- 2 * (s * n + 1)
c(test, (tmp2 * test)/s/s/tmp1, s * tmp1, tmp2)
}
Roy <- function (eig, q, df.res) {
p <- length(eig)
test <- max(eig)
tmp1 <- max(p, q)
tmp2 <- df.res - tmp1 + q
c(test, (tmp2 * test)/tmp1, tmp1, tmp2)
}
has.intercept.mlm <- function (model, ...)
any(row.names(coefficients(model)) == "(Intercept)")
Anova.mlm <- function(mod, type=c("II","III", 2, 3), SSPE, error.df, idata,
idesign, icontrasts=c("contr.sum", "contr.poly"), imatrix,
test.statistic=c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"),...){
wts <- if (!is.null(mod$weights)) mod$weights else rep(1, nrow(model.matrix(mod)))
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(SSPE)) SSPE <- wcrossprod(residuals(mod), w=wts)
if (missing(idata)) {
idata <- NULL
idesign <- NULL
}
if (missing(imatrix)) imatrix <- NULL
error.df <- if (missing(error.df)) df.residual(mod)
else error.df
switch(type,
II=Anova.II.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
III=Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
"2"=Anova.II.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...),
"3"=Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test.statistic, ...))
}
Anova.III.mlm <- function(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...){
intercept <- has.intercept(mod)
V <- solve(crossprod(model.matrix(mod)))
p <- nrow(coefficients(mod))
I.p <- diag(p)
terms <- term.names(mod)
n.terms <- length(terms)
assign <- mod$assign
if (is.null(idata) && is.null(imatrix)){
if ((n.terms == 0) && intercept) {
Test <- linearHypothesis(mod, 1, SSPE=SSPE, ...)
result <- list(SSP=Test$SSPH, SSPE=SSPE, df=1, error.df=error.df,
terms="(Intercept)", repeated=FALSE, type="III", test=test)
class(result) <- "Anova.mlm"
return(result)
}
SSP <- as.list(rep(0, n.terms))
df <- rep(0, n.terms)
names(df) <- names(SSP) <- terms
for (term in 1:n.terms){
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs,,drop=FALSE]
Test <- linearHypothesis(mod, hyp.matrix, SSPE=SSPE, ...)
SSP[[term]] <- Test$SSPH
df[term]<- length(subs)
}
result <- list(SSP=SSP, SSPE=SSPE, df=df, error.df=error.df, terms=terms,
repeated=FALSE, type="III", test=test)
}
else {
if (!is.null(imatrix)){
X.design <- do.call(cbind, imatrix)
ncols <- sapply(imatrix, ncol)
end <- cumsum(ncols)
start <- c(1, (end + 1))[-(length(end) + 1)]
cols <- mapply(seq, from=start, to=end)
iterms <- names(end)
names(cols) <- iterms
check.imatrix(X.design, iterms)
}
else {
if (is.null(idesign)) stop("idesign (intra-subject design) missing.")
for (i in 1:length(idata)){
if (is.null(attr(idata[,i], "contrasts"))){
contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2]
else icontrasts[1]
}
}
X.design <- model.matrix(idesign, data=idata)
i.intercept <- has.intercept(X.design)
iterms <- term.names(idesign)
if (i.intercept) iterms <- c("(Intercept)", iterms)
check.imatrix(X.design)
}
df <- rep(0, n.terms*length(iterms))
hnames <- rep("", length(df))
P <- SSPEH <- SSP <- as.list(df)
i <- 0
for (iterm in iterms){
for (term in 1:n.terms){
subs <- which(assign == term - intercept)
hyp.matrix <- I.p[subs,,drop=FALSE]
i <- i + 1
Test <- linearHypothesis(mod, hyp.matrix, SSPE=SSPE,
idata=idata, idesign=idesign, icontrasts=icontrasts, iterms=iterm,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)
SSP[[i]] <- Test$SSPH
SSPEH[[i]] <- Test$SSPE
P[[i]] <- Test$P
df[i] <- length(subs)
hnames[i] <- if (iterm == "(Intercept)") terms[term]
else if (terms[term] == "(Intercept)") iterm
else paste(terms[term], ":", iterm, sep="")
}
}
names(df) <- names(SSP) <- names(SSPEH) <- names(P) <- hnames
result <- list(SSP=SSP, SSPE=SSPEH, P=P, df=df, error.df=error.df,
terms=hnames, repeated=TRUE, type="III", test=test,
idata=idata, idesign=idesign, icontrasts=icontrasts, imatrix=imatrix,
singular=Test$singular)
}
class(result) <- "Anova.mlm"
result
}
Anova.II.mlm <- function(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...){
wts <- if (!is.null(mod$weights)) mod$weights else rep(1, nrow(model.matrix(mod)))
V <- solve(wcrossprod(model.matrix(mod), w=wts))
SSP.term <- function(term, iterm){
which.term <- which(term == terms)
subs.term <- which(assign == which.term)
relatives <- relatives(term, terms, 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.2 <- I.p[c(subs.relatives, subs.term),,drop=FALSE]
if (missing(iterm)){
SSP1 <- if (length(subs.relatives) == 0) 0
else linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V, singular.ok=TRUE, ...)$SSPH
SSP2 <- linearHypothesis(mod, hyp.matrix.2, SSPE=SSPE, V=V, singular.ok=TRUE, ...)$SSPH
return(SSP2 - SSP1)
}
else {
SSP1 <- if (length(subs.relatives) == 0) 0
else linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts, P=imatrix[[iterm]], singular.ok=TRUE, ...)$SSPH
lh2 <- linearHypothesis(mod, hyp.matrix.2, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts, P=imatrix[[iterm]], singular.ok=TRUE, ...)
return(list(SSP = lh2$SSPH - SSP1, SSPE=lh2$SSPE, P=lh2$P, singular=lh2$singular))
}
}
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- nrow(coefficients(mod))
I.p <- diag(p)
assign <- mod$assign
terms <- term.names(mod)
if (intercept) terms <- terms[-1]
n.terms <- length(terms)
if (n.terms == 0){
message("Note: model has only an intercept; equivalent type-III tests substituted.")
return(Anova.III.mlm(mod, SSPE, error.df, idata, idesign, icontrasts, imatrix, test, ...))
}
if (is.null(idata) && is.null(imatrix)){
SSP <- as.list(rep(0, n.terms))
df <- rep(0, n.terms)
names(df) <- names(SSP) <- terms
for (i in 1:n.terms){
SSP[[i]] <- SSP.term(terms[i])
df[i]<- df.terms(mod, terms[i])
}
result <- list(SSP=SSP, SSPE=SSPE, df=df, error.df=error.df, terms=terms,
repeated=FALSE, type="II", test=test)
}
else {
if (!is.null(imatrix)){
X.design <- do.call(cbind, imatrix)
ncols <- sapply(imatrix, ncol)
end <- cumsum(ncols)
start <- c(1, (end + 1))[-(length(end) + 1)]
cols <- mapply(seq, from=start, to=end)
iterms <- names(end)
names(cols) <- iterms
check.imatrix(X.design, iterms)
}
else {
if (is.null(idesign)) stop("idesign (intra-subject design) missing.")
for (i in 1:length(idata)){
if (is.null(attr(idata[,i], "contrasts"))){
contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2]
else icontrasts[1]
}
}
X.design <- model.matrix(idesign, data=idata)
iintercept <- has.intercept(X.design)
iterms <- term.names(idesign)
if (iintercept) iterms <- c("(Intercept)", iterms)
check.imatrix(X.design)
}
df <- rep(0, (n.terms + intercept)*length(iterms))
hnames <- rep("", length(df))
P <- SSPEH <- SSP <- as.list(df)
i <- 0
for (iterm in iterms){
if (intercept){
i <- i + 1
hyp.matrix.1 <- I.p[-1,,drop=FALSE]
SSP1 <- linearHypothesis(mod, hyp.matrix.1, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)$SSPH
lh2 <- linearHypothesis(mod, I.p, SSPE=SSPE, V=V,
idata=idata, idesign=idesign, iterms=iterm, icontrasts=icontrasts,
check.imatrix=FALSE, P=imatrix[[iterm]], singular.ok=TRUE, ...)
SSP[[i]] <- lh2$SSPH - SSP1
SSPEH[[i]] <- lh2$SSPE
P[[i]] <- lh2$P
df[i] <- 1
hnames[i] <- iterm
}
for (term in 1:n.terms){
subs <- which(assign == term)
i <- i + 1
Test <- SSP.term(terms[term], iterm)
SSP[[i]] <- Test$SSP
SSPEH[[i]] <- Test$SSPE
P[[i]] <- Test$P
df[i]<- length(subs)
hnames[i] <- if (iterm == "(Intercept)") terms[term]
else paste(terms[term], ":", iterm, sep="")
}
}
names(df) <- names(P) <- names(SSP) <- names(SSPEH) <- hnames
result <- list(SSP=SSP, SSPE=SSPEH, P=P, df=df, error.df=error.df,
terms=hnames, repeated=TRUE, type="II", test=test,
idata=idata, idesign=idesign, icontrasts=icontrasts, imatrix=imatrix,
singular=Test$singular)
}
class(result) <- "Anova.mlm"
result
}
print.Anova.mlm <- function(x, ...){
if ((!is.null(x$singular)) && x$singular) stop("singular error SSP matrix; multivariate tests unavailable\ntry summary(object, multivariate=FALSE)")
test <- x$test
repeated <- x$repeated
ntests <- length(x$terms)
tests <- matrix(NA, ntests, 4)
if (!repeated) SSPE.qr <- qr(x$SSPE)
for (term in 1:ntests){
# some of the code here adapted from stats:::summary.manova
eigs <- Re(eigen(qr.coef(if (repeated) qr(x$SSPE[[term]]) else SSPE.qr,
x$SSP[[term]]), symmetric = FALSE)$values)
tests[term, 1:4] <- switch(test,
Pillai = Pillai(eigs, x$df[term], x$error.df),
Wilks = Wilks(eigs, x$df[term], x$error.df),
"Hotelling-Lawley" = HL(eigs, x$df[term], x$error.df),
Roy = Roy(eigs, x$df[term], x$error.df))
}
ok <- tests[, 2] >= 0 & tests[, 3] > 0 & tests[, 4] > 0
ok <- !is.na(ok) & ok
tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3], tests[ok, 4],
lower.tail = FALSE))
rownames(tests) <- x$terms
colnames(tests) <- c("Df", "test stat", "approx F", "num Df", "den Df", "Pr(>F)")
tests <- structure(as.data.frame(tests),
heading = paste("\nType ", x$type, if (repeated) " Repeated Measures",
" MANOVA Tests: ", test, " test statistic", sep=""),
class = c("anova", "data.frame"))
print(tests)
invisible(x)
}
# path <- "D:/R-package-sources/car/R"
# files <- list.files(path, pattern=".*\\.R")
# files <- paste(path, files, sep="/")
# for (file in files) source(file)
# summary.Anova.mlm and print.summary.Anova.mlm methods
# with contributions from Gabriel Baud-Bovy
summary.Anova.mlm <- function (object, test.statistic, univariate=object$repeated, multivariate=TRUE, p.adjust.method, ...) {
GG <- function(SSPE, P) { # Greenhouse-Geisser correction
p <- nrow(SSPE)
if (p < 2)
return(NA)
lambda <- eigen(SSPE %*% solve(t(P) %*% P))$values
lambda <- lambda[lambda > 0]
((sum(lambda)/p)^2)/(sum(lambda^2)/p)
}
HF <- function(gg, error.df, p) { # Huynh-Feldt correction
((error.df + 1) * p * gg - 2)/(p * (error.df - p * gg))
}
mauchly <- function(SSD, P, df) {
# most of this function borrowed from stats:::mauchly.test.SSD
if (nrow(SSD) < 2)
return(c(NA, NA))
Tr <- function(X) sum(diag(X))
p <- nrow(P)
I <- diag(p)
Psi <- t(P) %*% I %*% P
B <- SSD
pp <- nrow(SSD)
U <- solve(Psi, B)
n <- df
logW <- log(det(U)) - pp * log(Tr(U/pp))
rho <- 1 - (2 * pp^2 + pp + 2)/(6 * pp * n)
w2 <- (pp + 2) * (pp - 1) * (pp - 2) * (2 * pp^3 + 6 *
pp^2 + 3 * p + 2)/(288 * (n * pp * rho)^2)
z <- -n * rho * logW
f <- pp * (pp + 1)/2 - 1
Pr1 <- pchisq(z, f, lower.tail = FALSE)
Pr2 <- pchisq(z, f + 4, lower.tail = FALSE)
pval <- Pr1 + w2 * (Pr2 - Pr1)
c(statistic = c(W = exp(logW)), p.value = pval)
}
if (missing(test.statistic))
test.statistic <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy")
test.statistic <- match.arg(test.statistic, c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), several.ok = TRUE)
nterms <- length(object$terms)
summary.object <- list(type=object$type, repeated=object$repeated,
multivariate.tests=NULL, univariate.tests=NULL,
pval.adjustments=NULL, sphericity.tests=NULL)
if (multivariate){
summary.object$multivariate.tests <- vector(nterms, mode="list")
names(summary.object$multivariate.tests) <- object$terms
summary.object$SSPE <- object$SSPE
for (term in 1:nterms) {
hyp <- list(SSPH = object$SSP[[term]],
SSPE = if (object$repeated) object$SSPE[[term]] else object$SSPE,
P = if (object$repeated) object$P[[term]] else NULL,
test = test.statistic, df = object$df[term],
df.residual = object$error.df, title = object$terms[term])
class(hyp) <- "linearHypothesis.mlm"
summary.object$multivariate.tests[[term]] <- hyp
}
}
if (object$repeated && univariate) {
singular <- object$singular
error.df <- object$error.df
table <- matrix(0, nterms, 6)
table2 <- matrix(0, nterms, 4)
table3 <- matrix(0, nterms, 2)
rownames(table3) <- rownames(table2) <- rownames(table) <- object$terms
colnames(table) <- c("SS", "num Df", "Error SS", "den Df", "F", "Pr(>F)")
colnames(table2) <- c("GG eps", "Pr(>F[GG])", "HF eps","Pr(>F[HF])")
colnames(table3) <- c("Test statistic", "p-value")
if (singular)
warning("Singular error SSP matrix:\nnon-sphericity test and corrections not available")
for (term in 1:nterms) {
SSP <- object$SSP[[term]]
SSPE <- object$SSPE[[term]]
P <- object$P[[term]]
p <- ncol(P)
PtPinv <- solve(t(P) %*% P)
gg <- if (!singular) GG(SSPE, P) else NA
table[term, "SS"] <- sum(diag(SSP %*% PtPinv))
table[term, "Error SS"] <- sum(diag(SSPE %*% PtPinv))
table[term, "num Df"] <- object$df[term] * p
table[term, "den Df"] <- error.df * p
table[term, "F"] <- (table[term, "SS"]/table[term, "num Df"])/
(table[term, "Error SS"]/table[term, "den Df"])
table[term, "Pr(>F)"] <- pf(table[term, "F"], table[term, "num Df"], table[term, "den Df"],
lower.tail = FALSE)
table2[term, "GG eps"] <- gg
table2[term, "HF eps"] <- if (!singular) HF(gg, error.df, p) else NA
table3[term, ] <- if (!singular) mauchly(SSPE, P, object$error.df) else NA
}
table3 <- na.omit(table3)
if (nrow(table3) > 0) {
table2[, "Pr(>F[GG])"] <- pf(table[, "F"], table2[, "GG eps"] *
table[, "num Df"], table2[, "GG eps"] * table[, "den Df"],
lower.tail = FALSE)
table2[, "Pr(>F[HF])"] <- pf(table[, "F"], pmin(1, table2[, "HF eps"]) *
table[, "num Df"], pmin(1, table2[, "HF eps"]) * table[, "den Df"],
lower.tail = FALSE)
table2 <- na.omit(table2)
if (any(table2[, "HF eps"] > 1)) warning("HF eps > 1 treated as 1")
}
class(table3) <- class(table) <- "anova"
summary.object$univariate.tests <- table
summary.object$pval.adjustments <- table2
summary.object$sphericity.tests <- table3
}
if (!object$repeated && univariate) {
SS <- sapply(object$SSP, diag)
SSE <- diag(object$SSPE)
df <- object$df
dfe <- object$error.df
F <- (SS/df)/(SSE/dfe)
SS <- cbind(SS, residuals=SSE)
SS <- rbind(df=c(df, residuals=dfe), SS)
p <- pf(F, df, dfe, lower.tail=FALSE)
result <- list(SS=t(SS), F=t(F), p=t(p), type=object$type)
if (!missing(p.adjust.method)){
if (isTRUE(p.adjust.method)) p.adjust.method <- "holm"
p.adj <- apply(p, 2, p.adjust, method=p.adjust.method)
result$p.adjust <- t(p.adj)
result$p.adjust.method <- p.adjust.method
}
class(result) = "univaov"
summary.object$univaov <- result
}
class(summary.object) <- "summary.Anova.mlm"
summary.object
}
print.summary.Anova.mlm <- function(x, digits = getOption("digits"), SSP=TRUE, SSPE=SSP, ... ) {
if (!is.null(x$multivariate.tests)) {
cat(paste("\nType ", x$type, if (x$repeated)
" Repeated Measures", " MANOVA Tests:\n", sep = ""))
if ((!x$repeated) && SSPE) {
cat("\nSum of squares and products for error:\n")
print(x$SSPE, digits = digits)
}
for (term in 1:length(x$multivariate.tests)) {
cat(paste("\n------------------------------------------\n",
"\nTerm:", names(x$multivariate.tests)[term], "\n"))
print(x$multivariate.tests[[term]], digits = digits, SSP=SSP, SSPE=FALSE, ...)
}
}
if (!is.null(x$univariate.tests)) {
cat("\nUnivariate Type", x$type, "Repeated-Measures ANOVA Assuming Sphericity\n\n")
print(x$univariate.tests)
if (nrow(x$sphericity.tests) > 0) {
cat("\n\nMauchly Tests for Sphericity\n\n")
print(x$sphericity.tests)
cat("\n\nGreenhouse-Geisser and Huynh-Feldt Corrections\n",
"for Departure from Sphericity\n\n")
table <- x$pval.adjustments[, 1:2, drop = FALSE]
class(table) <- "anova"
print(table)
cat("\n")
table <- x$pval.adjustments[, 3:4, drop = FALSE]
class(table)
print(table)
}
}
if (!is.null(x$univaov)){
print(x$univaov, ...)
}
invisible(x)
}
print.univaov <- function(x, digits = max(getOption("digits") - 2L, 3L),
style=c("wide", "long"),
by=c("response", "term"),
...){
style <- match.arg(style)
if (style == "wide") {
cat("\n Type", x$type, "Sums of Squares\n")
print(x$SS, digits=digits)
cat("\n F-tests\n")
F <- x$F
print(round(F, 2))
cat("\n p-values\n")
p <- format.pval(x$p)
p <- matrix(p, nrow=nrow(F))
rownames(p) <- rownames(F)
colnames(p) <- colnames(F)
print(p, quote=FALSE)
if (!is.null(x$p.adjust)){
cat("\n p-values adjusted (by term) for simultaneous inference by", x$p.adjust.method, "method\n")
p.adjust <- format.pval(x$p.adjust)
p.adjust <- matrix(p.adjust, nrow=nrow(F))
rownames(p.adjust) <- rownames(F)
colnames(p.adjust) <- colnames(F)
print(p.adjust, quote=FALSE)
}
}
else {
x.df <- as.data.frame(x, by=by)
x.df$F <- round(x.df$F, 2)
x.df$p <- format.pval(x.df$p)
if (!is.null(x$p.adjust)) x.df$"adjusted p" <- format.pval(x.df$"adjusted p")
cat("\n Type", x$type, "Sums of Squares and F tests\n")
print(x.df, quote=FALSE, digits=digits)
}
invisible(x)
}
as.data.frame.univaov <- function(x, row.names, optional, by=c("response", "term"), ...) {
melt <- function(data, varnames = names(dimnames(data)), value.name = "value") {
dn <- dimnames(data)
labels <- expand.grid( dn[[1]], dn[[2]])
colnames(labels) <- varnames
value_df <- setNames(data.frame(as.vector(data)), value.name)
cbind(labels, value_df)
}
nv <- ncol(x$F)
nt <- nrow(x$F)
by <- match.arg(by)
if (by=="response") {
vn <- c("term", "response")
df <- matrix(x$SS[1:nt, "df", drop=FALSE], nrow=nt, ncol=nv)
SS <- melt(x$SS[1:nt, -1, drop=FALSE], varnames=vn, value.name="SS")
F <- melt(x$F, varnames=vn, value.name="F")
p <- melt(x$p, varnames=vn, value.name="p")
if (!is.null(x$p.adjust)) p.adjust <- melt(x$p.adjust, varnames=vn, value.name="adjusted p")
}
else {
vn <- rev(c("term", "response"))
df <- t(matrix(x$SS[1:nt, "df", drop=FALSE], nrow=nt, ncol=nv))
SS <- melt(t(x$SS[1:nt, -1, drop=FALSE]), varnames=vn, value.name="SS")
F <- melt(t(x$F), varnames=vn, value.name="F")
p <- melt(t(x$p), varnames=vn, value.name="p")
if (!is.null(x$p.adjust)) p.adjust <- melt(x$p.adjust, varnames=vn, value.name="adjusted p")
}
result <- cbind(SS[,c(2,1,3)], df=c(df), F=F[,"F"], p=p[,"p"])
if (!is.null(x$p.adjust)) result <- cbind(result, "adjusted p"=p.adjust[, "adjusted p"])
result
}
Anova.manova <- function(mod, ...){
class(mod) <- c("mlm", "lm")
Anova(mod, ...)
}
Manova <- function(mod, ...){
UseMethod("Manova")
}
Manova.mlm <- function(mod, ...){
Anova(mod, ...)
}
# Cox regression models
df.residual.coxph <- function(object, ...){
object$n - sum(!is.na(coef(object)))
}
alias.coxph <- function(model){
if(any(which <- is.na(coef(model)))) return(list(Complete=which))
else list()
}
logLik.coxph <- function(object, ...) object$loglik[2]
Anova.coxph <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (length((mod$rscore) > 0) && (test.statistic == "LR")){
warning("LR tests unavailable with robust variances\nWald tests substituted")
test.statistic <- "Wald"
}
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.coxph(mod),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", vcov.=vcov(mod))),
III=switch(test.statistic,
LR=Anova.III.LR.coxph(mod),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", vcov.=vcov(mod))),
"2"=switch(test.statistic,
LR=Anova.II.LR.coxph(mod),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", vcov.=vcov(mod))),
"3"=switch(test.statistic,
LR=Anova.III.LR.coxph(mod),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", vcov.=vcov(mod))))
}
Anova.II.LR.coxph <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
method <- mod$method
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
mod.1 <- survival::coxph(mod$y ~ X[, -exclude.1, drop = FALSE], method=method)
loglik.1 <- logLik(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
survival::coxph(mod$y ~ X[, -exclude.2, drop = FALSE], method=method)
}
loglik.2 <- logLik(mod.2)
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
Anova.III.LR.coxph <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
method <- mod$method
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
df <- df.terms(mod)
LR <- p <- rep(0, n.terms)
loglik1 <- logLik(mod)
for (term in 1:n.terms){
mod.0 <- survival::coxph(mod$y ~ X[, -which.nms(names[term])], method=method)
LR[term] <- -2*(logLik(mod.0) - loglik1)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df","Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result,"heading") <- "Analysis of Deviance Table (Type III tests)"
result
}
# parametric survival regression models
alias.survreg <- function(model){
if(any(which <- diag(vcov(model)) < 1e-10)) return(list(Complete=which))
else list()
}
logLik.survreg <- function(object, ...) object$loglik[2]
Anova.survreg <- function(mod, type=c("II","III", 2, 3), test.statistic=c("LR", "Wald"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (length((mod$rscore) > 0) && (test.statistic == "LR")){
warning("LR tests unavailable with robust variances\nWald tests substituted")
test.statistic <- "Wald"
}
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.survreg(mod),
Wald=Anova.II.Wald.survreg(mod)),
III=switch(test.statistic,
LR=Anova.III.LR.survreg(mod),
Wald=Anova.III.Wald.survreg(mod)),
"2"=switch(test.statistic,
LR=Anova.II.LR.survreg(mod),
Wald=Anova.II.Wald.survreg(mod)),
"3"=switch(test.statistic,
LR=Anova.III.LR.survreg(mod),
Wald=Anova.III.Wald.survreg(mod)))
}
Anova.II.LR.survreg <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
dist <- mod$dist
scale <- mod$call$scale
weights <- model.frame(mod)$"(weights)"
arg.list <- list(dist=dist)
if (!is.null(scale)) arg.list$scale <- scale
if (!is.null(weights)) arg.list$weights <- weights
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
asgn <- asgn[asgn != 0]
if (has.intercept(mod)){
int <- which(names == "(Intercept)")
X <- X[, -int]
names <- names[-int]
}
n.terms <- length(names)
if (n.terms < 2) return(anova(mod))
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
y <- model.frame(mod)[,1]
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
exclude.1 <- as.vector(unlist(sapply(c(names[term], rels), which.nms)))
arg.list$formula <- y ~ X[, -exclude.1, drop = FALSE]
mod.1 <- do.call(survival::survreg, arg.list)
# mod.1 <- survival::survreg(y ~ X[, -exclude.1, drop = FALSE])
loglik.1 <- logLik(mod.1)
mod.2 <- if (length(rels) == 0) mod
else {
arg.list$formula <- y ~ X[, -exclude.2, drop = FALSE]
exclude.2 <- as.vector(unlist(sapply(rels, which.nms)))
do.call(survival::survreg, arg.list)
# survival::survreg(y ~ X[, -exclude.2, drop = FALSE])
}
loglik.2 <- logLik(mod.2)
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
Anova.III.LR.survreg <- function(mod, ...){
if (!requireNamespace("survival")) stop("survival package is missing")
dist <- mod$dist
scale <- mod$call$scale
weights <- model.frame(mod)$"(weights)"
arg.list <- list(dist=dist)
if (!is.null(scale)) arg.list$scale <- scale
if (!is.null(weights)) arg.list$weights <- weights
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
asgn <- asgn[asgn != 0]
if (has.intercept(mod)){
int <- which(names == "(Intercept)")
X <- X[, -int]
names <- names[-int]
}
n.terms <- length(names)
if (n.terms < 2) return(anova(mod))
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
y <- model.frame(mod)[,1]
loglik1 <- logLik(mod)
for (term in 1:n.terms){
arg.list$formula <- y ~ X[, -which.nms(names[term])]
mod.0 <- do.call(survival::survreg, arg.list)
# mod.0 <- survival::survreg(y ~ X[, -which.nms(names[term])])
LR[term] <- -2*(logLik(mod.0) - loglik1)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df","Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result,"heading") <- "Analysis of Deviance Table (Type III tests)"
result
}
Anova.II.Wald.survreg <- function(mod){
V <- vcov(mod)
p <- which(rownames(V) == "Log(scale)")
if (length(p) > 0) V <- V[-p, -p]
Anova.II.default(mod, V, test="Chisq")
}
Anova.III.Wald.survreg <- function(mod){
V <- vcov(mod)
p <- which(rownames(V) == "Log(scale)")
if (length(p) > 0) V <- V[-p, -p]
Anova.III.default(mod, V, test="Chisq")
}
# Default Anova() method: requires methods for vcov() (if vcov. argument not specified) and coef().
Anova.default <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
switch(type,
II=Anova.II.default(mod, vcov., test.statistic, singular.ok=singular.ok),
III=Anova.III.default(mod, vcov., test.statistic, singular.ok=singular.ok),
"2"=Anova.II.default(mod, vcov., test.statistic, singular.ok=singular.ok),
"3"=Anova.III.default(mod, vcov., test.statistic, singular.ok=singular.ok))
}
assignVector <- function(model){
m <- model.matrix(model)
assign <- attr(m, "assign")
if (!is.null(assign)) return (assign)
m <- model.matrix(formula(model), data=model.frame(model))
assign <- attr(m, "assign")
if (!has.intercept(model)) assign <- assign[assign != 0]
assign
}
Anova.II.default <- function(mod, vcov., test, singular.ok=TRUE, ...){
hyp.term <- function(term){
which.term <- which(term==names)
subs.term <- if (is.list(assign)) assign[[which.term]] else which(assign == which.term)
relatives <- relatives(term, names, fac)
subs.relatives <- NULL
for (relative in relatives){
sr <- if (is.list(assign)) assign[[relative]] else which(assign == relative)
subs.relatives <- c(subs.relatives, sr)
}
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.default(mod, hyp.matrix.term,
vcov.=vcov., test=test, singular.ok=singular.ok, ...)
if (test=="Chisq") c(statistic=hyp$Chisq[2], df=hyp$Df[2])
else c(statistic=hyp$F[2], df=hyp$Df[2])
}
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
assign <- assignVector(mod) # attr(model.matrix(mod), "assign")
if (!is.list(assign)) assign[!not.aliased] <- NA
else if (intercept) assign <- assign[-1]
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
df <- c(rep(0, n.terms), df.residual(mod))
if (inherits(mod, "coxph")){
assign <- assign[assign != 0]
clusters <- grep("^cluster\\(", names)
if (length(clusters) > 0) {
names <- names[-clusters]
df <- df[-clusters]
n.terms <- n.terms - length(clusters)
}
}
# if (inherits(mod, "plm")) assign <- assign[assign != 0]
p <- teststat <- rep(0, n.terms + 1)
teststat[n.terms + 1] <- p[n.terms + 1] <- NA
for (i in 1:n.terms){
hyp <- hyp.term(names[i])
teststat[i] <- abs(hyp["statistic"])
df[i] <- abs(hyp["df"])
p[i] <- if (test == "Chisq")
pchisq(teststat[i], df[i], lower.tail=FALSE)
else pf(teststat[i], df[i], df[n.terms + 1], lower.tail=FALSE)
}
result <- if (test == "Chisq"){
if (length(df) == n.terms + 1) df <- df[1:n.terms]
data.frame(df, teststat[!is.na(teststat)], p[!is.na(teststat)])
}
else data.frame(df, teststat, p)
if (nrow(result) == length(names) + 1) names <- c(names,"Residuals")
row.names(result) <- names
names(result) <- c ("Df", test, if (test == "Chisq") "Pr(>Chisq)"
else "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.default <- function(mod, vcov., test, singular.ok=FALSE, ...){
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- assignVector(mod) # attr(model.matrix(mod), "assign")
df <- c(rep(0, n.terms), df.residual(mod))
if (inherits(mod, "coxph")){
if (intercept) names <- names[-1]
assign <- assign[assign != 0]
clusters <- grep("^cluster\\(", names)
if (length(clusters) > 0) {
names <- names[-clusters]
df <- df[-clusters]
n.terms <- n.terms - length(clusters)
}
}
# if (inherits(mod, "plm")) assign <- assign[assign != 0]
if (intercept) df[1] <- sum(grepl("^\\(Intercept\\)", names(coef(mod))))
teststat <- rep(0, n.terms + 1)
p <- rep(0, n.terms + 1)
teststat[n.terms + 1] <- p[n.terms + 1] <- NA
not.aliased <- !is.na(coef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1:n.terms){
subs <- if (is.list(assign)) assign[[term]] else 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.default(mod, hyp.matrix,
vcov.=vcov., test=test, singular.ok=singular.ok, ...)
teststat[term] <- if (test=="Chisq") hyp$Chisq[2] else hyp$F[2]
df[term] <- abs(hyp$Df[2])
p[term] <- if (test == "Chisq")
pchisq(teststat[term], df[term], lower.tail=FALSE)
else pf(teststat[term], df[term], df[n.terms + 1], lower.tail=FALSE)
}
}
result <- if (test == "Chisq"){
if (length(df) == n.terms + 1) df <- df[1:n.terms]
data.frame(df, teststat[!is.na(teststat)], p[!is.na(teststat)])
}
else data.frame(df, teststat, p)
if (nrow(result) == length(names) + 1) names <- c(names,"Residuals")
row.names(result) <- names
names(result) <- c ("Df", test, if (test == "Chisq") "Pr(>Chisq)"
else "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III tests)\n",
paste("Response:", responseName(mod)))
result
}
## functions for mixed models
# the following function, not exported, to make car consistent with CRAN and development versions of lme4 and with nlme
fixef <- function (object){
if (isS4(object)) {
if (!inherits(object, "merMod")) object@fixef
else lme4::fixef(object)
}
else object$coefficients$fixed
}
Anova.merMod <- function(mod, type=c("II","III", 2, 3),
test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
Anova.mer(mod=mod, type=type, test.statistic=test.statistic, vcov.=vcov.,
singular.ok=singular.ok, ...)
}
Anova.mer <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Chisq", "F"),
vcov.=vcov(mod), singular.ok, ...){
if (is.function(vcov.)) vcov. <- vcov.(mod)
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
if (missing(singular.ok))
singular.ok <- type == "2" || type == "II"
switch(type,
II=Anova.II.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
III=Anova.III.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
"2"=Anova.II.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok),
"3"=Anova.III.mer(mod, test=test.statistic, vcov., singular.ok=singular.ok))
}
Anova.II.mer <- function(mod, vcov., singular.ok=TRUE, test=c("Chisq", "F"), ...){
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(mod, hyp.matrix.term,
vcov.=vcov., singular.ok=singular.ok, test=test, ...)
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]))
}
test <- match.arg(test)
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
if (!missing(vcov.)){
vcov. <- if (test == "F"){
# if (!require("pbkrtest")) stop("pbkrtest package required for F-tests on linear mixed model")
as.matrix(pbkrtest::vcovAdj(mod, details=0))
}
else vcov(mod)
}
assign <- attr(model.matrix(mod), "assign")
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
p <- teststat <- df <- res.df <- rep(0, n.terms)
for (i in 1: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] <- if (test == "Chisq") pchisq(teststat[i], df[i], lower.tail=FALSE)
else pf(teststat[i], df[i], res.df[i], lower.tail=FALSE)
}
if (test=="Chisq"){
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(mod)))
}
else {
result <- data.frame(teststat, df, res.df, p)
row.names(result) <- names
names(result) <- c ("F", "Df", "Df.res", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)\n",
paste("Response:", responseName(mod)))
}
result
}
Anova.III.mer <- function(mod, vcov., singular.ok=FALSE, test=c("Chisq", "F"), ...){
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- attr(model.matrix(mod), "assign")
p <- teststat <- df <- res.df <- rep(0, n.terms)
if (intercept) df[1] <- 1
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
if (!missing(vcov.)){
vcov. <- if (test == "F"){
# if (!require("pbkrtest")) stop("pbkrtest package required for F-tests on linear mixed model")
as.matrix(pbkrtest::vcovAdj(mod, details=0))
}
else vcov(mod)
}
for (term in 1: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(mod, hyp.matrix, test=test,
vcov.=vcov., singular.ok=singular.ok, ...)
if (test == "Chisq"){
teststat[term] <- hyp$Chisq[2]
df[term] <- abs(hyp$Df[2])
p[term] <- pchisq(teststat[term], df[term], lower.tail=FALSE)
}
else{
teststat[term] <- hyp$F[2]
df[term] <- abs(hyp$Df[2])
res.df[term]=hyp$Res.Df[2]
p[term] <- pf(teststat[term], df[term], res.df[term], lower.tail=FALSE)
}
}
}
if (test == "Chisq"){
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(mod)))
}
else {
result <- data.frame(teststat, df, res.df, p)
row.names(result) <- names
names(result) <- c ("F", "Df", "Df.res", "Pr(>F)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- c("Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)\n",
paste("Response:", responseName(mod)))
}
result
}
Anova.lme <- function(mod, type=c("II","III", 2, 3),
vcov.=vcov(mod), singular.ok, ...){
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"
switch(type,
II=Anova.II.lme(mod, vcov., singular.ok=singular.ok),
III=Anova.III.lme(mod, vcov., singular.ok=singular.ok),
"2"=Anova.II.lme(mod, vcov., singular.ok=singular.ok),
"3"=Anova.III.lme(mod, vcov., singular.ok=singular.ok))
}
Anova.II.lme <- function(mod, vcov., singular.ok=TRUE, ...){
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(mod, hyp.matrix.term,
vcov.=vcov., singular.ok=singular.ok, ...)
c(statistic=hyp$Chisq[2], df=hyp$Df[2])
}
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
fac <- attr(terms(mod), "factors")
intercept <- has.intercept(mod)
p <- length(fixef(mod))
I.p <- diag(p)
assign <- attr(model.matrix(mod), "assign")
assign[!not.aliased] <- NA
names <- term.names(mod)
if (intercept) names <- names[-1]
n.terms <- length(names)
p <- teststat <- df <- rep(0, n.terms)
for (i in 1:n.terms){
hyp <- hyp.term(names[i])
teststat[i] <- abs(hyp["statistic"])
df[i] <- abs(hyp["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 tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.III.lme <- function(mod, vcov., singular.ok=FALSE, ...){
intercept <- has.intercept(mod)
p <- length(coefficients(mod))
I.p <- diag(p)
names <- term.names(mod)
n.terms <- length(names)
assign <- attr(model.matrix(mod), "assign")
df <- rep(0, n.terms)
if (intercept) df[1] <- 1
p <- teststat <-rep(0, n.terms)
not.aliased <- !is.na(fixef(mod))
if (!singular.ok && !all(not.aliased))
stop("there are aliased coefficients in the model")
for (term in 1: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(mod, hyp.matrix,
vcov.=vcov., singular.ok=singular.ok, ...)
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 tests)\n",
paste("Response:", responseName(mod)))
result
}
Anova.svyglm <- function(mod, ...) Anova.default(mod, ...)
Anova.rlm <- function(mod, ...) Anova.default(mod, test.statistic="F", ...)
Anova.coxme <- function(mod, type=c("II","III", 2, 3), test.statistic=c("Wald", "LR"), ...){
type <- as.character(type)
type <- match.arg(type)
test.statistic <- match.arg(test.statistic)
switch(type,
II=switch(test.statistic,
LR=Anova.II.LR.coxme(mod, ...),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", ...)),
III=switch(test.statistic,
LR=stop("type-III LR tests not available for coxme models"),
Wald=Anova.default(mod, type="III", test.statistic="Chisq", ...)),
"2"=switch(test.statistic,
LR=Anova.II.LR.coxme(mod, ...),
Wald=Anova.default(mod, type="II", test.statistic="Chisq", ...)),
"3"=switch(test.statistic,
LR=stop("type-III LR tests not available for coxme models"),
Wald=Anova.default(mod, type="III", test.statistic="Chisq")))
}
Anova.II.LR.coxme <- function(mod, ...){
if (!requireNamespace("coxme")) stop("coxme package is missing")
which.nms <- function(name) which(asgn == which(names == name))
fac <-attr(terms(mod), "factors")
names <- term.names(mod)
n.terms <- length(names)
if (n.terms < 2) return(anova(mod, test="Chisq"))
X <- model.matrix(mod)
asgn <- attr(X, 'assign')
p <- LR <- rep(0, n.terms)
df <- df.terms(mod)
random <- mod$formulaList$random
random <- sapply(random, as.character)[2, ]
random <- paste(paste0("(", random, ")"), collapse=" + ")
fixed <- as.character(mod$formulaList$fixed)[3]
for (term in 1:n.terms){
rels <- names[relatives(names[term], names, fac)]
formula <- paste0(". ~ . - ", paste(c(names[term], rels), collapse=" - "), " + ", random)
mod.1 <- update(mod, as.formula(formula))
loglik.1 <- logLik(mod.1, type="integrated")
mod.2 <- if (length(rels) == 0) mod
else {
formula <- paste0(". ~ . - ", paste(rels, collapse=" - "), " + ", random)
update(mod, as.formula(formula))
}
loglik.2 <- logLik(mod.2, type="integrated")
LR[term] <- -2*(loglik.1 - loglik.2)
p[term] <- pchisq(LR[term], df[term], lower.tail=FALSE)
}
result <- data.frame(LR, df, p)
row.names(result) <- names
names(result) <- c("LR Chisq", "Df", "Pr(>Chisq)")
class(result) <- c("anova", "data.frame")
attr(result, "heading") <- "Analysis of Deviance Table (Type II tests)"
result
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.