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