1 2anova.glmrob <- function(object, ..., 3 test = c("Wald", "QD", "QDapprox")) 4{ 5 dotargs <- list(...) 6 if (!is.null(names(dotargs))) { 7 named <- (names(dotargs) != "") 8 if (any(named)) { 9 warning("the following arguments to 'anova.glmrob' are invalid and", 10 "dropped:\n", 11 pasteK(deparse(dotargs[named]))) 12 dotargs <- dotargs[!named] 13 } 14 } 15 is.glmrob <- vapply(dotargs, inherits, NA, what="glmrob") 16 if(!all(is.glmrob) || !inherits(object, "glmrob")) 17 stop("anova.glmrob() only works for 'glmrob' objects") 18 test <- match.arg(test) 19 if (length(dotargs) > 0) 20 anovaGlmrobList(c(list(object), dotargs), test=test) 21 else { 22 ## 23 ## "'Anova Table' for a single model object 24 stop("'Anova Table' for a single model object not yet implemented") 25 } 26} 27 28 29anovaGlmrobList <- function (object, test=NULL) 30{ 31 nmodels <- length(object) 32 stopifnot(nmodels >= 2) 33 responses <- as.character(lapply(object, 34 function(x) deparse(formula(x)[[2]]))) 35 if (!all(responses == responses[1])) 36 stop("Not the same response used in the fitted models") 37 nobs <- sapply(object, function(x) length(x$residuals)) 38 if (any(nobs != nobs[1])) 39 stop("models were not all fitted to the same size of dataset") 40 methods <- as.character(lapply(object, function(x) x$method)) 41 if(!all(methods == methods[1])) 42 stop("Not the same method used for fitting the models") 43 note <- paste("Models fitted by method '", methods[1], "'", sep="") 44 tccs <- sapply(object, function(x) length(x$tcc)) 45 if(!all(tccs == tccs[1])) 46 stop("Not the same tuning constant c used in the robust fits") 47 ## 48 tbl <- matrix(rep(NA, nmodels*4), ncol = 4) 49 tbl[1,1] <- nobs[1] - length(coef(object[[1]])) 50 for(k in 2:nmodels) 51 tbl[k,] <- anovaGlmrobPair(object[[k-1]], object[[k]], test=test) 52 53 ## return 54 dimnames(tbl) <- list(1:nmodels, 55 c("pseudoDf", "Test.Stat", "Df", "Pr(>chisq)")) 56 title <- switch(test, 57 Wald = "Robust Wald Test Table", 58 QD = "Robust Quasi-Deviance Table", 59 QDapprox = 60 "Robust Quasi-Deviance Table Based on a Quadratic Approximation", 61 "") 62 variables <- lapply(object, function(x) 63 paste(deparse(formula(x)), collapse = "\n")) 64 topnote <- paste("Model ", format(1:nmodels), ": ", variables, 65 sep = "", collapse = "\n") 66 structure(as.data.frame(tbl), heading = c(title, "", topnote, note,""), 67 class = c("anova", "data.frame")) 68} 69 70 71anovaGlmrobPair <- function(obj1, obj2, test) 72{ 73 if(length(coef(obj1)) < length(coef(obj2))){ 74 Sign <- 1 75 full.mfit <- obj2 76 reduced.mfit <- obj1 77 } 78 else { 79 Sign <- -1 80 full.mfit <- obj1 81 reduced.mfit <- obj2 82 } 83 X <- model.matrix(full.mfit) 84 asgn <- attr(X, "assign") 85 tt <- terms(full.mfit) 86 tt0 <- terms(reduced.mfit) 87 tl <- attr(tt, "term.labels") 88 tl0 <- attr(tt0, "term.labels") 89 numtl0 <- match(tl0 , tl, nomatch = -1) 90 if(attr(tt0, "intercept") == 1) numtl0 <- c(0, numtl0) 91 if(any(is.na(match(numtl0, unique(asgn))))) 92 stop("Models are not nested!") 93 mod0 <- seq(along = asgn)[!is.na(match(asgn, numtl0))] 94 if (length(asgn) == length(mod0)) 95 stop("Models are not strictly nested") 96 H0ind <- setdiff(seq(along = asgn), mod0) 97 H0coef <- coef(full.mfit)[H0ind] 98 df <- length(H0coef) 99 pp <- df + length(mod0) 100 101 if(test == "Wald") { 102 t.cov <- full.mfit$cov 103 t.chisq <- sum(H0coef * solve(t.cov[H0ind, H0ind], H0coef)) 104 statistic <- c(chisq = t.chisq) 105 } 106 else if(full.mfit$method=="Mqle" && (test == "QD" || test == "QDapprox")) { 107 matM <- full.mfit$matM 108 if(test == "QDapprox") { 109 ## Difference of robust quasi-deviances 110 ## via the asymptotically equivalent quadratic form 111 matM11 <- matM[mod0, mod0, drop=FALSE] 112 matM12 <- matM[mod0, H0ind, drop=FALSE] 113 matM22 <- matM[H0ind, H0ind, drop=FALSE] 114 matM22.1 <- matM22 - crossprod(matM12, solve(matM11, matM12)) 115 Dquasi.dev <- nrow(X) * c(H0coef %*% matM22.1 %*% H0coef) 116 } 117 else { 118 quasiDev <- switch(full.mfit$family$family, 119 poisson = glmrobMqleDiffQuasiDevPois, 120 binomial = glmrobMqleDiffQuasiDevB, 121 Gamma = glmrobMqleDiffQuasiDevGamma, 122 stop("This family is not implemented")) 123 124 ## note that qdev and qdev0 do depend on an incorrectly specified 125 ## lower limits in the integration. But this does't matter in 126 ## the following difference, because the difference does not 127 ## deepend on it! (Hence I could use the centered nui 128 ## (cnui= nui - Enui) in quasiDev as the function to be integrated. 129 Dquasi.dev <- quasiDev(mu = full.mfit$fitted.values, 130 mu0 = reduced.mfit$fitted.values, 131 y = full.mfit$y, ni = full.mfit$ni, 132 w.x = full.mfit$w.x, phi=full.mfit$dispersion, 133 tcc = full.mfit$tcc) 134 } 135 ## Asymptotic distribution: variance and weights of the sum of chi2 136 matQ <- full.mfit$matQ 137 matM11inv <- solve(matM[mod0,mod0]) 138 Mplus <- matrix(0, ncol = pp, nrow = pp) 139 Mplus[mod0, mod0] <- matM11inv 140 141 d.ev <- Re(eigen(matQ %*% (solve(matM)-Mplus), only.values=TRUE)$values) 142 d.ev <- d.ev[1:df] ## just the q (=df) lagest eigenvalues are needed 143 144 if(any(d.ev < 0)) warning("some eigenvalues are negative") 145 146 ## p-value: exact computation for q=1, approximated for q>1 (q=df) 147 statistic <- c(quasi.dev = Dquasi.dev/mean(d.ev)) 148 149 } else stop("non-implemented test method:", test, 150 "for fitting method", full.mfit$method) 151 152 ## return 153 c(nrow(X)-pp+df*(Sign<0), Sign*statistic, Sign*df, 154 pchisq(as.vector(statistic), df=df, lower.tail = FALSE)) 155} 156 157