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