anova.glmrob.RdCompute an analysis of robust quasi-deviance table for one or more
generalized linear models fitted by glmrob.
objects of class glmrob, typically
the result of a call to glmrob.
a character string specifying the test statistic to be
used. (Partially) matching one of "Wald", "QD" or
"QDapprox". See Details.
Specifying a single object gives a sequential analysis of robust quasi-deviance table for that fit. That is, the reductions in the robust residual quasi-deviance as each term of the formula is added in turn are given in as the rows of a table. (Currently not yet implemented.)
If more than one object is specified, the table has a row for the residual quasi-degrees of freedom (However, this information is never used in the asymptotic tests). For all but the first model, the change in degrees of freedom and robust quasi-deviance is also given. (This only makes statistical sense if the models are nested.) It is conventional to list the models from smallest to largest, but this is up to the user.
In addition, the table will contain test statistics and P values
comparing the reduction in robust quasi-deviance for the model on the
row to that on top of it. For all robust fitting methods, the
“Wald”-type test between two models can be applied (test
= "Wald").
When using Mallows or Huber type robust estimators
(method="Mqle" in glmrob), then there are
additional test methods. One is the robust quasi-deviance test
(test = "QD"), as described by Cantoni and Ronchetti (2001).
The asymptotic distribution is approximated by a chi-square
distibution. Another test (test = "QDapprox") is based on a
quadratic approximation of the robust quasi-deviance test
statistic. Its asymptotic distribution is chi-square (see the reference).
The comparison between two or more models by anova.glmrob
will only be valid if they are fitted to the same dataset and by the same
robust fitting method using the same tuning constant \(c\) (tcc in
glmrob).
Basically, an object of class anova inheriting from class
data.frame.
E. Cantoni and E. Ronchetti (2001) Robust Inference for Generalized Linear Models. JASA 96 (455), 1022–1030.
E.Cantoni (2004) Analysis of Robust Quasi-deviances for Generalized Linear Models. Journal of Statistical Software 10, https://www.jstatsoft.org/article/view/v010i04
## Binomial response -----------
data(carrots)
Cfit2 <- glmrob(cbind(success, total-success) ~ logdose + block,
family=binomial, data=carrots, method="Mqle",
control=glmrobMqle.control(tcc=1.2))
summary(Cfit2)
#>
#> Call: glmrob(formula = cbind(success, total - success) ~ logdose + block, family = binomial, data = carrots, method = "Mqle", control = glmrobMqle.control(tcc = 1.2))
#>
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.3883 0.6923 3.450 0.000561 ***
#> logdose -2.0491 0.3685 -5.561 2.68e-08 ***
#> blockB2 0.2351 0.2122 1.108 0.267828
#> blockB3 -0.4496 0.2409 -1.866 0.061989 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> Robustness weights w.r * w.x:
#> 15 weights are ~= 1. The remaining 9 ones are
#> 2 5 6 7 13 14 21 22 23
#> 0.7756 0.7026 0.6751 0.9295 0.8536 0.2626 0.8337 0.9051 0.9009
#>
#> Number of observations: 24
#> Fitted by method ‘Mqle’ (in 9 iterations)
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> No deviance values available
#> Algorithmic parameters:
#> acc tcc
#> 0.0001 1.2000
#> maxit
#> 50
#> test.acc
#> "coef"
#>
Cfit4 <- glmrob(cbind(success, total-success) ~ logdose * block,
family=binomial, data=carrots, method="Mqle",
control=glmrobMqle.control(tcc=1.2))
anova(Cfit2, Cfit4, test="Wald")
#> Robust Wald Test Table
#>
#> Model 1: cbind(success, total - success) ~ logdose + block
#> Model 2: cbind(success, total - success) ~ logdose * block
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 20
#> 2 18 2.4895 2 0.288
anova(Cfit2, Cfit4, test="QD")
#> Robust Quasi-Deviance Table
#>
#> Model 1: cbind(success, total - success) ~ logdose + block
#> Model 2: cbind(success, total - success) ~ logdose * block
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 20
#> 2 18 2.8412 2 0.2416
anova(Cfit2, Cfit4, test="QDapprox")
#> Robust Quasi-Deviance Table Based on a Quadratic Approximation
#>
#> Model 1: cbind(success, total - success) ~ logdose + block
#> Model 2: cbind(success, total - success) ~ logdose * block
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 20
#> 2 18 2.5022 2 0.2862
## Poisson response ------------
data(epilepsy)
Efit2 <- glmrob(Ysum ~ Age10 + Base4*Trt, family=poisson, data=epilepsy,
method="Mqle", control=glmrobMqle.control(tcc=1.2,maxit=100))
summary(Efit2)
#>
#> Call: glmrob(formula = Ysum ~ Age10 + Base4 * Trt, family = poisson, data = epilepsy, method = "Mqle", control = glmrobMqle.control(tcc = 1.2, maxit = 100))
#>
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.036768 0.154168 13.211 < 2e-16 ***
#> Age10 0.158434 0.047444 3.339 0.000840 ***
#> Base4 0.085132 0.004174 20.395 < 2e-16 ***
#> Trtprogabide -0.323886 0.087421 -3.705 0.000211 ***
#> Base4:Trtprogabide 0.011842 0.004967 2.384 0.017124 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> Robustness weights w.r * w.x:
#> 26 weights are ~= 1. The remaining 33 ones are summarized as
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.07328 0.30750 0.50730 0.49220 0.68940 0.97240
#>
#> Number of observations: 59
#> Fitted by method ‘Mqle’ (in 14 iterations)
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> No deviance values available
#> Algorithmic parameters:
#> acc tcc
#> 0.0001 1.2000
#> maxit
#> 100
#> test.acc
#> "coef"
#>
Efit3 <- glmrob(Ysum ~ Age10 + Base4 + Trt, family=poisson, data=epilepsy,
method="Mqle", control=glmrobMqle.control(tcc=1.2,maxit=100))
anova(Efit3, Efit2, test = "Wald")
#> Robust Wald Test Table
#>
#> Model 1: Ysum ~ Age10 + Base4 + Trt
#> Model 2: Ysum ~ Age10 + Base4 * Trt
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 55
#> 2 54 5.6836 1 0.01712 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(Efit3, Efit2, test = "QD")
#> Robust Quasi-Deviance Table
#>
#> Model 1: Ysum ~ Age10 + Base4 + Trt
#> Model 2: Ysum ~ Age10 + Base4 * Trt
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 55
#> 2 54 2.9691 1 0.08487 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## trivial intercept-only-model:
E0 <- update(Efit3, . ~ 1)
anova(E0, Efit3, Efit2, test = "QDapprox")
#> Robust Quasi-Deviance Table Based on a Quadratic Approximation
#>
#> Model 1: Ysum ~ 1
#> Model 2: Ysum ~ Age10 + Base4 + Trt
#> Model 3: Ysum ~ Age10 + Base4 * Trt
#> Models fitted by method 'Mqle'
#>
#> pseudoDf Test.Stat Df Pr(>chisq)
#> 1 58
#> 2 55 2006.26 3 < 2e-16 ***
#> 3 54 5.68 1 0.01712 *
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1