# written by C. Scherber on 19th Jan 2015 # based on Code by John Fox (Anova function in "car" package) # the Anova function now works also for multinom models with large (user-specified) weights # usage: mymodel=multinom(...,MaxNWts=xxx); Anova.multinom2(mymodel,MaxNWts=xxx) responseName=function (model, ...) { UseMethod("responseName") } responseName.default=function (model, ...) deparse(attr(terms(model), "variables")[[2]]) 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))] } term.names=function (model, ...) { UseMethod("term.names") } term.names.default=function (model, ...) { term.names <- labels(terms(model)) c("(Intercept)", term.names) } df.terms=function (model, ...) { UseMethod("df.terms") } is.aliased=function (model) { !is.null(alias(model)$Complete) } alias=function (object, ...) UseMethod("alias") alias.formula=function (object, data, ...) { lm.obj <- if (missing(data)) aov(object) else aov(object, data) alias(lm.obj, ...) } alias.multinom=function (object, complete = TRUE, partial = FALSE, partial.pattern = FALSE, ...) { CompPatt <- function(x, ...) { x[abs(x) < 1e-06] <- 0 MASS::fractions(x) } PartPatt <- function(x) { z <- zapsmall(x) != 0 if (any(z)) { xx <- abs(signif(x[z], 2)) ll <- length(unique(xx)) if (ll > 10L) xx <- cut(xx, 9L) else if (ll == 1L) x[] <- 1 x[z] <- paste0(ifelse(x[z] > 0, " ", "-"), xx) } x[!z] <- "" collabs <- colnames(x) collabs <- if (length(collabs)) abbreviate(sub("\\.", "", collabs), 3L) else 1L:ncol(x) colnames(x) <- collabs class(x) <- "mtable" x } Model <- object$terms attributes(Model) <- NULL value <- list(Model = Model) R <- qr(object)$qr R <- R[1L:min(dim(R)), , drop = FALSE] R[lower.tri(R)] <- 0 d <- dim(R) rank <- object$rank p <- d[2L] if (complete) { value$Complete <- if (is.null(p) || rank == p) NULL else { p1 <- 1L:rank X <- R[p1, p1] Y <- R[p1, -p1, drop = FALSE] beta12 <- as.matrix(qr.coef(qr(X), Y)) CompPatt(t(beta12)) } } if (partial) { tmp <- suppressWarnings(summary.lm(object)$cov.unscaled) ses <- sqrt(diag(tmp)) beta11 <- tmp/outer(ses, ses) beta11[row(beta11) >= col(beta11)] <- 0 beta11[abs(beta11) < 1e-06] <- 0 if (all(beta11 == 0)) beta11 <- NULL else if (partial.pattern) beta11 <- PartPatt(beta11) value$Partial <- beta11 } class(value) <- "listof" value } df.terms.default=function (model, term, ...) { #if (is.aliased(model)) # stop("Model has aliased term(s); df ambiguous.") if (!missing(term) && 1 == length(term)) { assign <- attr(model.matrix(model), "assign") which.term <- which(term == labels(terms(model))) if (0 == length(which.term)) stop(paste(term, "is not in the model.")) sum(assign == which.term) } else { terms <- if (missing(term)) labels(terms(model)) else term result <- numeric(0) for (term in terms) result <- c(result, Recall(model, term)) names(result) <- terms result } } Anova.multinom2=function (mod, MaxNWts,...) { which.nms <- function(name) which(asgn == which(names == name)) fac <- attr(mod$terms, "factors") names <- term.names(mod)[-1] 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,MaxNWts=MaxNWts, trace = FALSE) else multinom(y ~ 1, weights = wt, race = FALSE,MaxNWts=MaxNWts) 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, MaxNWts=MaxNWts, 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 }