1#--------------------------------------------------------------------------------------- 2# Revision history: 3# 2009-01-16: replaced unlist(options("foo")) with getOption("foo") 4# 2009-09-16: optionally allow models with aliased coefficients. J. Fox 5# 2009-12-10: modification by A. Zeileis to allow wider range of coef. names. 6# 2009-12-22: small changes to linearHypothesis.mlm() to handle user-specified 7# within-subjects designs in Anova() 8# 2010-05-21: linearHypothesis.default() and .lm() changed so that differences 9# in df, etc. will be postive. 10# 2010-06-12: linearHypothesis.mlm() changed to allow observation weights 11# 2010-06-22: fixed bug in linearHypothesis.lm caused by 2010-05-21 revision 12# 2010-01-21: added methods for mixed models; added matchCoefs() and methods. J. Fox 13# 2011-05-03: fixed bug in displaying numbers starting with "-1" or "+1" in printed representation. J. Fox 14# 2011-06-09: added matchCoefs.mlm(). J. Fox 15# 2011-11-27: added linearHypothesis.svyglm(). John 16# 2011-12-27: fixed printing bug in linearHypothesis(). John 17# 2012-02-28: added F-test to linearHypothesis.mer(). John 18# 2012-03-07: singular.ok argument added to linearHypothesis.mlm(). J. Fox 19# 2012-08-20: Fixed p-value bug for chisq test in .mer method. John 20# 2012-09-17: updated linearHypothesis.mer for pkrtest 0.3-2. John 21# 2012-11-21: test for NULL rhs to avoid warning in R 2.16.0. John 22# 2013-01-28: hypotheses can now contain newlines and tabs 23# 2013-02-14: fixed bug in printing constants of the form 1.x*. John 24# 2013-06-20: added .merMod() method. John 25# 2013-06-22: tweaks for lme4. John 26# 2013-06-22: test argument uniformly uses "Chisq" rather than "chisq". J. Fox 27# 2013-08-19: removed calls to unexported functions in stats. J. Fox 28# 2014-08-17: added call to requireNamespace() and :: as needed (doesn't work for pbkrtest). J. Fox 29# 2014-08-18: fixed bug in linearHypothesis.survreg(). J. Fox 30# 2014-09-23: added linearHypothesis.rlm. J. Fox 31# 2014-12-18: check that residual df nonzero in Anova.lm() and Anova.default 32# and residual SS nonzero in Anova.lm(). John 33# 2015-01-27: KRmodcomp() and methods now imported from pbkrtest. John 34# 2015-02-03: Check for NULL df before 0 df in default method. John 35# 2016-06-29: added "value" and "vcov" attributes to returned object, print vcov when verbose. John 36# 2017-02-16: replaced KRmodcomp() with pbkrtest::KRmodcomp(). John 37# 2017-11-07: added complete=FALSE to vcov() calls. John 38# 2019-06-06: remove vcov.default(), which is no longer needed, suggestion of Pavel Krivitsky. John 39# 2020-05-27: tweak to linearHypothesis.survreg(). John 40# 2020-12-21: regularize handling of vcov. arg. Sandy and John 41# 2020-12-21: new matchCoefs.lmList() method, which covers nlsList objects. John 42# 2020-12-21: added linearHypothesis.lmList(). John 43#---------------------------------------------------------------------------------------------------- 44 45# vcov.default <- function(object, ...){ 46# stop(paste("there is no vcov() method for models of class", 47# paste(class(object), collapse=", "))) 48# } 49 50has.intercept.matrix <- function (model, ...) { 51 "(Intercept)" %in% colnames(model) 52} 53 54 55makeHypothesis <- function(cnames, hypothesis, rhs = NULL){ 56 parseTerms <- function(terms){ 57 component <- gsub("^[-\\ 0-9\\.]+", "", terms) 58 component <- gsub(" ", "", component, fixed=TRUE) 59 component 60 } 61 stripchars <- function(x) { 62 x <- gsub("\\n", " ", x) 63 x <- gsub("\\t", " ", x) 64 x <- gsub(" ", "", x, fixed = TRUE) 65 x <- gsub("*", "", x, fixed = TRUE) 66 x <- gsub("-", "+-", x, fixed = TRUE) 67 x <- strsplit(x, "+", fixed = TRUE)[[1]] 68 x <- x[x!=""] 69 x 70 } 71 char2num <- function(x) { 72 x[x == ""] <- "1" 73 x[x == "-"] <- "-1" 74 as.numeric(x) 75 } 76 constants <- function(x, y) { 77 with.coef <- unique(unlist(sapply(y, 78 function(z) which(z == parseTerms(x))))) 79 if (length(with.coef) > 0) x <- x[-with.coef] 80 x <- if (is.null(x)) 0 else sum(as.numeric(x)) 81 if (any(is.na(x))) 82 stop('The hypothesis "', hypothesis, 83 '" is not well formed: contains bad coefficient/variable names.') 84 x 85 } 86 coefvector <- function(x, y) { 87 rv <- gsub(" ", "", x, fixed=TRUE) == 88 parseTerms(y) 89 if (!any(rv)) return(0) 90 if (sum(rv) > 1) stop('The hypothesis "', hypothesis, 91 '" is not well formed.') 92 rv <- sum(char2num(unlist(strsplit(y[rv], x, fixed=TRUE)))) 93 if (is.na(rv)) 94 stop('The hypothesis "', hypothesis, 95 '" is not well formed: contains non-numeric coefficients.') 96 rv 97 } 98 99 if (!is.null(rhs)) rhs <- rep(rhs, length.out = length(hypothesis)) 100 if (length(hypothesis) > 1) 101 return(rbind(Recall(cnames, hypothesis[1], rhs[1]), 102 Recall(cnames, hypothesis[-1], rhs[-1]))) 103 104 cnames_symb <- sapply(c("@", "#", "~"), function(x) length(grep(x, cnames)) < 1) 105 106 if(any(cnames_symb)) { 107 cnames_symb <- head(c("@", "#", "~")[cnames_symb], 1) 108 cnames_symb <- paste(cnames_symb, seq_along(cnames), cnames_symb, sep = "") 109 hypothesis_symb <- hypothesis 110 for(i in order(nchar(cnames), decreasing = TRUE)) 111 hypothesis_symb <- gsub(cnames[i], cnames_symb[i], hypothesis_symb, fixed = TRUE) 112 } else { 113 stop('The hypothesis "', hypothesis, 114 '" is not well formed: contains non-standard coefficient names.') 115 } 116 117 lhs <- strsplit(hypothesis_symb, "=", fixed=TRUE)[[1]] 118 if (is.null(rhs)) { 119 if (length(lhs) < 2) rhs <- "0" 120 else if (length(lhs) == 2) { 121 rhs <- lhs[2] 122 lhs <- lhs[1] 123 } 124 else stop('The hypothesis "', hypothesis, 125 '" is not well formed: contains more than one = sign.') 126 } 127 else { 128 if (length(lhs) < 2) as.character(rhs) 129 else stop('The hypothesis "', hypothesis, 130 '" is not well formed: contains a = sign although rhs was specified.') 131 } 132 lhs <- stripchars(lhs) 133 rhs <- stripchars(rhs) 134 rval <- sapply(cnames_symb, coefvector, y = lhs) - sapply(cnames_symb, coefvector, y = rhs) 135 rval <- c(rval, constants(rhs, cnames_symb) - constants(lhs, cnames_symb)) 136 names(rval) <- c(cnames, "*rhs*") 137 rval 138} 139 140printHypothesis <- function(L, rhs, cnames){ 141 hyp <- rep("", nrow(L)) 142 for (i in 1:nrow(L)){ 143 sel <- L[i,] != 0 144 h <- L[i, sel] 145 h <- ifelse(h < 0, as.character(h), paste("+", h, sep="")) 146 nms <- cnames[sel] 147 h <- paste(h, nms) 148 h <- gsub("-", " - ", h) 149 h <- gsub("+", " + ", h, fixed=TRUE) 150 h <- paste(h, collapse="") 151 h <- gsub(" ", " ", h, fixed=TRUE) 152 h <- sub("^\\ \\+", "", h) 153 h <- sub("^\\ ", "", h) 154 h <- sub("^-\\ ", "-", h) 155 h <- paste(" ", h, sep="") 156 h <- paste(h, "=", rhs[i]) 157 h <- gsub(" 1([^[:alnum:]_.]+)[ *]*", "", 158 gsub("-1([^[:alnum:]_.]+)[ *]*", "-", 159 gsub("- +1 +", "-1 ", h))) 160 h <- sub("Intercept)", "(Intercept)", h) 161 h <- gsub("-", " - ", h) 162 h <- gsub("+", " + ", h, fixed=TRUE) 163 h <- gsub(" ", " ", h, fixed=TRUE) 164 h <- sub("^ *", "", h) 165 hyp[i] <- h 166 } 167 hyp 168} 169 170linearHypothesis <- function (model, ...) 171 UseMethod("linearHypothesis") 172 173lht <- function (model, ...) 174 UseMethod("linearHypothesis") 175 176linearHypothesis.lmList <- function(model, ..., vcov.=vcov, coef.=coef){ 177 vcov.List <- function(object, ...) { 178 vlist <- lapply(object, vcov) 179 ng <- length(vlist) 180 nv <- dim(vlist[[1]])[1] 181 v <- matrix(0, nrow=ng*nv, ncol=ng*nv) 182 for (j in 1:ng){ 183 cells <- ((j-1)*nv + 1):(j*nv) 184 v[cells, cells] <- vlist[[j]] 185 } 186 v 187 } 188 suppress.vcov.msg <- missing(vcov.) 189 if (!is.function(vcov.)) stop("vcov. must be a function") 190 if (!is.function(coef.)) stop("coef. must be a function") 191 linearHypothesis.default(model, vcov.=vcov.List(model), 192 coef.=unlist(lapply(model, coef.)), suppress.vcov.msg = suppress.vcov.msg, ...) 193 } 194 195linearHypothesis.nlsList <- function(model, ..., vcov.=vcov, coef.=coef){ 196 NextMethod() 197} 198 199linearHypothesis.default <- function(model, hypothesis.matrix, rhs=NULL, 200 test=c("Chisq", "F"), vcov.=NULL, singular.ok=FALSE, verbose=FALSE, 201 coef. = coef(model), suppress.vcov.msg=FALSE, ...){ 202 df <- df.residual(model) 203 if (is.null(df)) df <- Inf ## if no residual df available 204 if (df == 0) stop("residual df = 0") 205 V <- if (is.null(vcov.)) vcov(model, complete=FALSE) 206 else if (is.function(vcov.)) vcov.(model) else vcov. 207 b <- coef. 208 if (any(aliased <- is.na(b)) && !singular.ok) 209 stop("there are aliased coefficients in the model") 210 b <- b[!aliased] 211 if (is.null(b)) stop(paste("there is no coef() method for models of class", 212 paste(class(model), collapse=", "))) 213 if (is.character(hypothesis.matrix)) { 214 L <- makeHypothesis(names(b), hypothesis.matrix, rhs) 215 if (is.null(dim(L))) L <- t(L) 216 rhs <- L[, NCOL(L)] 217 L <- L[, -NCOL(L), drop = FALSE] 218 rownames(L) <- hypothesis.matrix 219 } 220 else { 221 L <- if (is.null(dim(hypothesis.matrix))) t(hypothesis.matrix) 222 else hypothesis.matrix 223 if (is.null(rhs)) rhs <- rep(0, nrow(L)) 224 } 225 q <- NROW(L) 226 value.hyp <- L %*% b - rhs 227 vcov.hyp <- L %*% V %*% t(L) 228 if (verbose){ 229 cat("\nHypothesis matrix:\n") 230 print(L) 231 cat("\nRight-hand-side vector:\n") 232 print(rhs) 233 cat("\nEstimated linear function (hypothesis.matrix %*% coef - rhs)\n") 234 print(drop(value.hyp)) 235 cat("\n") 236 if (length(vcov.hyp) == 1) cat("\nEstimated variance of linear function\n") 237 else cat("\nEstimated variance/covariance matrix for linear function\n") 238 print(drop(vcov.hyp)) 239 cat("\n") 240 } 241 SSH <- as.vector(t(value.hyp) %*% solve(vcov.hyp) %*% value.hyp) 242 test <- match.arg(test) 243 if (!(is.finite(df) && df > 0)) test <- "Chisq" 244 name <- try(formula(model), silent = TRUE) 245 if (inherits(name, "try-error")) name <- substitute(model) 246 title <- "Linear hypothesis test\n\nHypothesis:" 247 topnote <- paste("Model 1: restricted model","\n", "Model 2: ", 248 paste(deparse(name), collapse = "\n"), sep = "") 249 note <- if (is.null(vcov.) || suppress.vcov.msg) "" 250 else "\nNote: Coefficient covariance matrix supplied.\n" 251 rval <- matrix(rep(NA, 8), ncol = 4) 252 colnames(rval) <- c("Res.Df", "Df", test, paste("Pr(>", test, ")", sep = "")) 253 rownames(rval) <- 1:2 254 rval[,1] <- c(df+q, df) 255 if (test == "F") { 256 f <- SSH/q 257 p <- pf(f, q, df, lower.tail = FALSE) 258 rval[2, 2:4] <- c(q, f, p) 259 } 260 else { 261 p <- pchisq(SSH, q, lower.tail = FALSE) 262 rval[2, 2:4] <- c(q, SSH, p) 263 } 264 if (!(is.finite(df) && df > 0)) rval <- rval[,-1] 265 result <- structure(as.data.frame(rval), 266 heading = c(title, printHypothesis(L, rhs, names(b)), "", topnote, note), 267 class = c("anova", "data.frame")) 268 attr(result, "value") <- value.hyp 269 attr(result, "vcov") <- vcov.hyp 270 result 271} 272 273linearHypothesis.glm <- function(model, ...) 274 linearHypothesis.default(model, ...) 275 276linearHypothesis.lm <- function(model, hypothesis.matrix, rhs=NULL, 277 test=c("F", "Chisq"), vcov.=NULL, 278 white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4"), 279 singular.ok=FALSE, ...){ 280 if (df.residual(model) == 0) stop("residual df = 0") 281 if (deviance(model) < sqrt(.Machine$double.eps)) stop("residual sum of squares is 0 (within rounding error)") 282 if (!singular.ok && is.aliased(model)) 283 stop("there are aliased coefficients in the model.") 284 test <- match.arg(test) 285 white.adjust <- as.character(white.adjust) 286 white.adjust <- match.arg(white.adjust) 287 if (white.adjust != "FALSE"){ 288 if (white.adjust == "TRUE") white.adjust <- "hc3" 289 vcov. <- hccm(model, type=white.adjust) 290 } 291 rval <- linearHypothesis.default(model, hypothesis.matrix, rhs = rhs, 292 test = test, vcov. = vcov., singular.ok=singular.ok, ...) 293 if (is.null(vcov.)) { 294 rval2 <- matrix(rep(NA, 4), ncol = 2) 295 colnames(rval2) <- c("RSS", "Sum of Sq") 296 SSH <- rval[2,test] 297 if (test == "F") SSH <- SSH * abs(rval[2, "Df"]) 298 df <- rval[2, "Res.Df"] 299 error.SS <- deviance(model) 300 rval2[,1] <- c(error.SS + SSH * error.SS/df, error.SS) 301 rval2[2,2] <- abs(diff(rval2[,1])) 302 rval2 <- cbind(rval, rval2)[,c(1, 5, 2, 6, 3, 4)] 303 class(rval2) <- c("anova", "data.frame") 304 attr(rval2, "heading") <- attr(rval, "heading") 305 attr(rval2, "value") <- attr(rval, "value") 306 attr(rval2, "vcov") <- attr(rval, "vcov") 307 rval <- rval2 308 } 309 rval 310} 311 312 313check.imatrix <- function(X, terms){ 314# check block orthogonality of within-subjects model matrix 315 XX <- crossprod(X) 316 if (missing(terms)) terms <- attr(X, "assign") 317 for (term in unique(terms)){ 318 subs <- term == terms 319 XX[subs, subs] <- 0 320 } 321 if (any(abs(XX) > sqrt(.Machine$double.eps))) 322 stop("Terms in the intra-subject model matrix are not orthogonal.") 323} 324 325linearHypothesis.mlm <- function(model, hypothesis.matrix, rhs=NULL, SSPE, V, 326 test, idata, icontrasts=c("contr.sum", "contr.poly"), idesign, iterms, 327 check.imatrix=TRUE, P=NULL, title="", singular.ok=FALSE, verbose=FALSE, ...){ 328 if (missing(test)) test <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy") 329 test <- match.arg(test, c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), 330 several.ok=TRUE) 331 df.residual <- df.residual(model) 332 wts <- if (!is.null(model$weights)) model$weights else rep(1,nrow(model.matrix(model))) 333 # V = (X'WX)^{-1} 334 if (missing (V)) V <- solve(wcrossprod(model.matrix(model), w=wts)) 335 B <- coef(model) 336 if (is.character(hypothesis.matrix)) { 337 L <- makeHypothesis(rownames(B), hypothesis.matrix, rhs) 338 if (is.null(dim(L))) L <- t(L) 339 L <- L[, -NCOL(L), drop = FALSE] 340 rownames(L) <- hypothesis.matrix 341 } 342 else { 343 L <- if (is.null(dim(hypothesis.matrix))) t(hypothesis.matrix) 344 else hypothesis.matrix 345 } 346 # SSPE = E'WE 347 if (missing(SSPE)) SSPE <- wcrossprod(residuals(model),w=wts) 348 if (missing(idata)) idata <- NULL 349 if (missing(idesign)) idesign <- NULL 350 if (!is.null(idata)){ 351 for (i in 1:length(idata)){ 352 if (is.null(attr(idata[,i], "contrasts"))){ 353 contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2] 354 else icontrasts[1] 355 } 356 } 357 if (is.null(idesign)) stop("idesign (intra-subject design) missing.") 358 X.design <- model.matrix(idesign, data=idata) 359 if (check.imatrix) check.imatrix(X.design) 360 intercept <- has.intercept(X.design) 361 term.names <- term.names(idesign) 362 if (intercept) term.names <- c("(Intercept)", term.names) 363 which.terms <- match(iterms, term.names) 364 if (any(nas <- is.na(which.terms))){ 365 if (sum(nas) == 1) 366 stop('The term "', iterms[nas],'" is not in the intrasubject design.') 367 else stop("The following terms are not in the intrasubject design: ", 368 paste(iterms[nas], collapse=", "), ".") 369 } 370 select <- apply(outer(which.terms, attr(X.design, "assign") + intercept, "=="), 371 2, any) 372 P <- X.design[, select, drop=FALSE] 373 } 374 if (!is.null(P)){ 375 rownames(P) <- colnames(B) 376 SSPE <- t(P) %*% SSPE %*% P 377 B <- B %*% P 378 } 379 rank <- sum(eigen(SSPE, only.values=TRUE)$values >= sqrt(.Machine$double.eps)) 380 if (!singular.ok && rank < ncol(SSPE)) 381 stop("The error SSP matrix is apparently of deficient rank = ", 382 rank, " < ", ncol(SSPE)) 383 r <- ncol(B) 384 if (is.null(rhs)) rhs <- matrix(0, nrow(L), r) 385 rownames(rhs) <- rownames(L) 386 colnames(rhs) <- colnames(B) 387 q <- NROW(L) 388 if (verbose){ 389 cat("\nHypothesis matrix:\n") 390 print(L) 391 cat("\nRight-hand-side matrix:\n") 392 print(rhs) 393 cat("\nEstimated linear function (hypothesis.matrix %*% coef - rhs):\n") 394 print(drop(L %*% B - rhs)) 395 cat("\n") 396 } 397 SSPH <- t(L %*% B - rhs) %*% solve(L %*% V %*% t(L)) %*% (L %*% B - rhs) 398 rval <- list(SSPH=SSPH, SSPE=SSPE, df=q, r=r, df.residual=df.residual, P=P, 399 title=title, test=test, singular=rank < ncol(SSPE)) 400 class(rval) <- "linearHypothesis.mlm" 401 rval 402} 403 404#linearHypothesis.mlm <- function(model, hypothesis.matrix, rhs=NULL, SSPE, V, 405# test, idata, icontrasts=c("contr.sum", "contr.poly"), idesign, iterms, 406# check.imatrix=TRUE, P=NULL, title="", verbose=FALSE, ...){ 407# if (missing(test)) test <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy") 408# test <- match.arg(test, c("Pillai", "Wilks", "Hotelling-Lawley", "Roy"), 409# several.ok=TRUE) 410# df.residual <- df.residual(model) 411# if (missing (V)) V <- solve(crossprod(model.matrix(model))) 412# B <- coef(model) 413# if (is.character(hypothesis.matrix)) { 414# L <- makeHypothesis(rownames(B), hypothesis.matrix, rhs) 415# if (is.null(dim(L))) L <- t(L) 416# L <- L[, -NCOL(L), drop = FALSE] 417# rownames(L) <- hypothesis.matrix 418# } 419# else { 420# L <- if (is.null(dim(hypothesis.matrix))) t(hypothesis.matrix) 421# else hypothesis.matrix 422# } 423# if (missing(SSPE)) SSPE <- crossprod(residuals(model)) 424# if (missing(idata)) idata <- NULL 425# if (missing(idesign)) idesign <- NULL 426# if (!is.null(idata)){ 427# for (i in 1:length(idata)){ 428# if (is.null(attr(idata[,i], "contrasts"))){ 429# contrasts(idata[,i]) <- if (is.ordered(idata[,i])) icontrasts[2] 430# else icontrasts[1] 431# } 432# } 433# if (is.null(idesign)) stop("idesign (intra-subject design) missing.") 434# X.design <- model.matrix(idesign, data=idata) 435# if (check.imatrix) check.imatrix(X.design) 436# intercept <- has.intercept(X.design) 437# term.names <- term.names(idesign) 438# if (intercept) term.names <- c("(Intercept)", term.names) 439# which.terms <- match(iterms, term.names) 440# if (any(nas <- is.na(which.terms))){ 441# if (sum(nas) == 1) 442# stop('The term "', iterms[nas],'" is not in the intrasubject design.') 443# else stop("The following terms are not in the intrasubject design: ", 444# paste(iterms[nas], collapse=", "), ".") 445# } 446# select <- apply(outer(which.terms, attr(X.design, "assign") + intercept, "=="), 447# 2, any) 448# P <- X.design[, select, drop=FALSE] 449# } 450# if (!is.null(P)){ 451# rownames(P) <- colnames(B) 452# SSPE <- t(P) %*% SSPE %*% P 453# B <- B %*% P 454# } 455# rank <- sum(eigen(SSPE, only.values=TRUE)$values >= sqrt(.Machine$double.eps)) 456# if (rank < ncol(SSPE)) 457# stop("The error SSP matrix is apparently of deficient rank = ", 458# rank, " < ", ncol(SSPE)) 459# r <- ncol(B) 460# if (is.null(rhs)) rhs <- matrix(0, nrow(L), r) 461# rownames(rhs) <- rownames(L) 462# colnames(rhs) <- colnames(B) 463# q <- NROW(L) 464# if (verbose){ 465# cat("\nHypothesis matrix:\n") 466# print(L) 467# cat("\nRight-hand-side matrix:\n") 468# print(rhs) 469# cat("\nEstimated linear function (hypothesis.matrix %*% coef - rhs):\n") 470# print(drop(L %*% B - rhs)) 471# cat("\n") 472# } 473# SSPH <- t(L %*% B - rhs) %*% solve(L %*% V %*% t(L)) %*% (L %*% B - rhs) 474# rval <- list(SSPH=SSPH, SSPE=SSPE, df=q, r=r, df.residual=df.residual, P=P, 475# title=title, test=test) 476# class(rval) <- "linearHypothesis.mlm" 477# rval 478#} 479 480print.linearHypothesis.mlm <- function(x, SSP=TRUE, SSPE=SSP, 481 digits=getOption("digits"), ...){ 482 test <- x$test 483 if (!is.null(x$P) && SSP){ 484 P <- x$P 485 cat("\n Response transformation matrix:\n") 486 attr(P, "assign") <- NULL 487 attr(P, "contrasts") <- NULL 488 print(P, digits=digits) 489 } 490 if (SSP){ 491 cat("\nSum of squares and products for the hypothesis:\n") 492 print(x$SSPH, digits=digits) 493 } 494 if (SSPE){ 495 cat("\nSum of squares and products for error:\n") 496 print(x$SSPE, digits=digits) 497 } 498 if ((!is.null(x$singular)) && x$singular){ 499 warning("the error SSP matrix is singular; multivariate tests are unavailable") 500 return(invisible(x)) 501 } 502 SSPE.qr <- qr(x$SSPE) 503 # the following code is adapted from summary.manova 504 eigs <- Re(eigen(qr.coef(SSPE.qr, x$SSPH), symmetric = FALSE)$values) 505 tests <- matrix(NA, 4, 4) 506 rownames(tests) <- c("Pillai", "Wilks", "Hotelling-Lawley", "Roy") 507 if ("Pillai" %in% test) 508 tests[1, 1:4] <- Pillai(eigs, x$df, x$df.residual) 509 if ("Wilks" %in% test) 510 tests[2, 1:4] <- Wilks(eigs, x$df, x$df.residual) 511 if ("Hotelling-Lawley" %in% test) 512 tests[3, 1:4] <- HL(eigs, x$df, x$df.residual) 513 if ("Roy" %in% test) 514 tests[4, 1:4] <- Roy(eigs, x$df, x$df.residual) 515 tests <- na.omit(tests) 516 ok <- tests[, 2] >= 0 & tests[, 3] > 0 & tests[, 4] > 0 517 ok <- !is.na(ok) & ok 518 tests <- cbind(x$df, tests, pf(tests[ok, 2], tests[ok, 3], tests[ok, 4], 519 lower.tail = FALSE)) 520 colnames(tests) <- c("Df", "test stat", "approx F", "num Df", "den Df", "Pr(>F)") 521 tests <- structure(as.data.frame(tests), 522 heading = paste("\nMultivariate Test", 523 if (nrow(tests) > 1) "s", ": ", x$title, sep=""), 524 class = c("anova", "data.frame")) 525 print(tests, digits=digits) 526 invisible(x) 527} 528 529linearHypothesis.survreg <- function(model, hypothesis.matrix, rhs=NULL, 530 test=c("Chisq", "F"), vcov., verbose=FALSE, ...){ 531 suppress.vcov.msg <- missing(vcov.) 532 if (missing(vcov.)) { 533 vcov. <- vcov(model, complete=FALSE) 534 b <- coef(model) 535 if (length(b) != nrow(vcov.)){ 536 p <- which(rownames(vcov.) == "Log(scale)") 537 if (length(p) > 0) vcov. <- vcov.[-p, -p] 538 } 539 } 540 linearHypothesis.default(model, hypothesis.matrix, rhs, test, vcov., verbose=verbose, 541 suppress.vcov.msg = suppress.vcov.msg, ...) 542} 543 544linearHypothesis.polr <- function (model, hypothesis.matrix, rhs=NULL, vcov., verbose=FALSE, ...){ 545 suppress.vcov.msg <- missing(vcov.) 546 k <- length(coef(model)) 547# V <- vcov(model, complete=FALSE)[1:k, 1:k] 548 V <- getVcov(vcov., model, complete=FALSE)[1:k, 1:k] 549 linearHypothesis.default(model, hypothesis.matrix, rhs, vcov.=V, verbose=verbose, 550 suppress.vcov.msg = suppress.vcov.msg, ...) 551} 552 553coef.multinom <- function(object, ...){ 554 # the following local function is copied from nnet:::coef.multinom 555 coef.m <- function (object, ...) { 556 r <- length(object$vcoefnames) 557 if (length(object$lev) == 2L) { 558 coef <- object$wts[1L + (1L:r)] 559 names(coef) <- object$vcoefnames 560 } 561 else { 562 coef <- matrix(object$wts, nrow = object$n[3L], byrow = TRUE)[, 563 1L + (1L:r), drop = FALSE] 564 if (length(object$lev)) 565 dimnames(coef) <- list(object$lev, object$vcoefnames) 566 if (length(object$lab)) 567 dimnames(coef) <- list(object$lab, object$vcoefnames) 568 coef <- coef[-1L, , drop = FALSE] 569 } 570 coef 571 } 572 573 b <- coef.m(object, ...) 574 cn <- colnames(b) 575 rn <- rownames(b) 576 b <- as.vector(t(b)) 577 names(b) <- as.vector(outer(cn, rn, function(c, r) paste(r, c, sep=":"))) 578 b 579} 580 581## functions for mixed models 582 583linearHypothesis.merMod <- function(model, hypothesis.matrix, rhs=NULL, 584 vcov.=NULL, test=c("Chisq", "F"), 585 singular.ok=FALSE, verbose=FALSE, ...){ 586 linearHypothesis.mer(model=model, hypothesis.matrix=hypothesis.matrix, 587 vcov.=vcov., test=test, singular.ok=singular.ok, 588 verbose=verbose, ...) 589} 590 591linearHypothesis.mer <- function(model, hypothesis.matrix, rhs=NULL, 592 vcov.=NULL, test=c("Chisq", "F"), singular.ok=FALSE, verbose=FALSE, ...){ 593 test <- match.arg(test) 594 V <- as.matrix(if (is.null(vcov.))vcov(model, complete=FALSE) 595 else if (is.function(vcov.)) vcov.(model) else vcov.) 596 b <- fixef(model) 597 if (any(aliased <- is.na(b)) && !singular.ok) 598 stop("there are aliased coefficients in the model") 599 b <- b[!aliased] 600 if (is.character(hypothesis.matrix)) { 601 L <- makeHypothesis(names(b), hypothesis.matrix, rhs) 602 if (is.null(dim(L))) L <- t(L) 603 rhs <- L[, NCOL(L)] 604 L <- L[, -NCOL(L), drop = FALSE] 605 rownames(L) <- hypothesis.matrix 606 } 607 else { 608 L <- if (is.null(dim(hypothesis.matrix))) t(hypothesis.matrix) 609 else hypothesis.matrix 610 if (is.null(rhs)) rhs <- rep(0, nrow(L)) 611 } 612 q <- NROW(L) 613 if (verbose){ 614 cat("\nHypothesis matrix:\n") 615 print(L) 616 cat("\nRight-hand-side vector:\n") 617 print(rhs) 618 cat("\nEstimated linear function (hypothesis.matrix %*% coef - rhs)\n") 619 print(drop(L %*% b - rhs)) 620 cat("\n") 621 } 622 if (test == "Chisq"){ 623 df <- Inf 624 SSH <- as.vector(t(L %*% b - rhs) %*% solve(L %*% V %*% t(L)) %*% (L %*% b - rhs)) 625 } 626 else { 627 if (!requireNamespace("lme4")) stop("lme4 package is missing") 628 if (!lme4::isREML(model)) 629 stop("F test available only for linear mixed model fit by REML") 630 if (!all(rhs == 0)) warning("rhs of hypothesis ignored, set to 0") 631 res <- pbkrtest::KRmodcomp(model, L)$test 632 df <- res["Ftest", "ddf"] 633 F <- res["Ftest", "stat"] 634 p <- res["Ftest", "p.value"] 635 } 636 name <- try(formula(model), silent = TRUE) 637 if (inherits(name, "try-error")) name <- substitute(model) 638 title <- "Linear hypothesis test\n\nHypothesis:" 639 topnote <- paste("Model 1: restricted model","\n", "Model 2: ", 640 paste(deparse(name), collapse = "\n"), sep = "") 641 note <- if (is.null(vcov.)) "" 642 else "\nNote: Coefficient covariance matrix supplied.\n" 643 rval <- matrix(rep(NA, 8), ncol = 4) 644 if (test == "Chisq"){ 645 colnames(rval) <- c("Res.Df", "Df", "Chisq", paste("Pr(>Chisq)", sep = "")) 646 rownames(rval) <- 1:2 647 rval[,1] <- c(df+q, df) 648 p <- pchisq(SSH, q, lower.tail = FALSE) 649 rval[2, 2:4] <- c(q, SSH, p) 650 rval <- rval[,-1] 651 } 652 else{ 653 colnames(rval) <- c("Res.Df", "Df", "F", paste("Pr(>F)", sep = "")) 654 rownames(rval) <- 1:2 655 rval[,1] <- c(df+q, df) 656 rval[2, 2:4] <- c(q, F, p) 657 } 658 structure(as.data.frame(rval), 659 heading = c(title, printHypothesis(L, rhs, names(b)), "", topnote, note), 660 class = c("anova", "data.frame")) 661} 662 663linearHypothesis.lme <- function(model, hypothesis.matrix, rhs=NULL, 664 vcov.=NULL, singular.ok=FALSE, verbose=FALSE, ...){ 665 V <- as.matrix(if (is.null(vcov.))vcov(model, complete=FALSE) 666 else if (is.function(vcov.)) vcov.(model) else vcov.) 667 b <- fixef(model) 668 if (any(aliased <- is.na(b)) && !singular.ok) 669 stop("there are aliased coefficients in the model") 670 b <- b[!aliased] 671 if (is.character(hypothesis.matrix)) { 672 L <- makeHypothesis(names(b), hypothesis.matrix, rhs) 673 if (is.null(dim(L))) L <- t(L) 674 rhs <- L[, NCOL(L)] 675 L <- L[, -NCOL(L), drop = FALSE] 676 rownames(L) <- hypothesis.matrix 677 } 678 else { 679 L <- if (is.null(dim(hypothesis.matrix))) t(hypothesis.matrix) 680 else hypothesis.matrix 681 if (is.null(rhs)) rhs <- rep(0, nrow(L)) 682 } 683 q <- NROW(L) 684 if (verbose){ 685 cat("\nHypothesis matrix:\n") 686 print(L) 687 cat("\nRight-hand-side vector:\n") 688 print(rhs) 689 cat("\nEstimated linear function (hypothesis.matrix %*% coef - rhs)\n") 690 print(drop(L %*% b - rhs)) 691 cat("\n") 692 } 693 df <- Inf 694 SSH <- as.vector(t(L %*% b - rhs) %*% solve(L %*% V %*% t(L)) %*% (L %*% b - rhs)) 695 name <- try(formula(model), silent = TRUE) 696 if (inherits(name, "try-error")) name <- substitute(model) 697 title <- "Linear hypothesis test\n\nHypothesis:" 698 topnote <- paste("Model 1: restricted model","\n", "Model 2: ", 699 paste(deparse(name), collapse = "\n"), sep = "") 700 note <- if (is.null(vcov.)) "" 701 else "\nNote: Coefficient covariance matrix supplied.\n" 702 rval <- matrix(rep(NA, 8), ncol = 4) 703 colnames(rval) <- c("Res.Df", "Df", "Chisq", paste("Pr(>Chisq)", sep = "")) 704 rownames(rval) <- 1:2 705 rval[,1] <- c(df+q, df) 706 p <- pchisq(SSH, q, lower.tail = FALSE) 707 rval[2, 2:4] <- c(q, SSH, p) 708 rval <- rval[,-1] 709 structure(as.data.frame(rval), 710 heading = c(title, printHypothesis(L, rhs, names(b)), "", topnote, note), 711 class = c("anova", "data.frame")) 712} 713 714## for svyglm 715 716linearHypothesis.svyglm <- function(model, ...) linearHypothesis.default(model, ...) 717 718## for rlm 719 720df.residual.rlm <- function(object, ...){ 721 p <- length(coef(object)) 722 wt.method <- object$call$wt.method 723 if (!is.null(wt.method) && wt.method == "case") { 724 sum(object$weights) - p 725 } 726 else length(object$wresid) - p 727} 728 729linearHypothesis.rlm <- function(model, ...) linearHypothesis.default(model, test="F", ...) 730 731 732## matchCoefs 733 734matchCoefs <- function(model, pattern, ...) UseMethod("matchCoefs") 735 736matchCoefs.default <- function(model, pattern, coef.=coef, ...){ 737 names <- names(coef.(model)) 738 grep(pattern, names, value=TRUE) 739} 740 741matchCoefs.mer <- function(model, pattern, ...) NextMethod(coef.=fixef) 742 743matchCoefs.merMod <- function(model, pattern, ...) NextMethod(coef.=fixef) 744 745matchCoefs.lme <- function(model, pattern, ...) NextMethod(coef.=fixef) 746 747matchCoefs.mlm <- function(model, pattern, ...){ 748 names <- rownames(coef(model)) 749 grep(pattern, names, value=TRUE) 750} 751 752matchCoefs.lmList <- function(model, pattern, ...){ 753 names <- names(unlist(lapply(model, coef))) 754 grep(pattern, names, value=TRUE) 755} 756 757 758