Robust Fitting of Generalized Linear Models
Description
glmrob
is used to fit generalized linear models by robust
methods. The models are specified by giving a symbolic description of
the linear predictor and a description of the error distribution.
Currently, robust methods are implemented for family =
binomial
, = poisson
, = Gamma
and = gaussian
.
Usage
1 2 3 4 
Arguments
formula 
a 
family 
a description of the error distribution and link function to
be used in the model. This can be a character string naming a
family function, a family 
data 
an optional data frame containing the variables
in the model. If not found in 
weights 
an optional vector of weights to be used in the fitting process. 
subset 
an optional vector specifying a subset of observations to be used in the fitting process. 
na.action 
a function which indicates what should happen
when the data contain 
start 
starting values for the parameters in the linear
predictor. Note that specifying 
offset 
this can be used to specify an a priori known component to be included in the linear predictor during fitting. 
method 
a character string specifying the robust fitting method. The details of method specification are given below. 
weights.on.x 
a character string (can be abbreviated), a If If it is a 
control 
a list of parameters for controlling the fitting process.
See the documentation for 
model 
a logical value indicating whether model frame should be included as a component of the returned value. 
x, y 
logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value. 
contrasts 
an optional list. See the 
trace.lev 
logical (or integer) indicating if intermediate results
should be printed; defaults to 
... 
arguments passed to 
Details
method="model.frame"
returns the model.frame()
,
the same as glm()
.
method="Mqle"
fits a generalized linear
model using Mallows or Huber type robust estimators, as described in
Cantoni and Ronchetti (2001) and Cantoni and Ronchetti (2006). In
contrast to the implementation
described in Cantoni (2004), the pure influence algorithm is
implemented.
method="WBY"
and method="BY"
,
available for logistic regression (family = binomial
) only, call
BYlogreg(*, initwml= . )
for the (weighted) BiancoYohai
estimator, where initwml
is true for "WBY"
, and false
for "BY"
.
method="MT"
, currently only implemented for family = poisson
,
computes an “[M]Estimator based on [T]ransformation”,
by Valdora and Yohai (2013).
weights.on.x= "robCov"
makes sense if all explanatory variables
are continuous.
In the cases,where weights.on.x
is "covMcd"
or
"robCov"
, or list with a “robCov” function, the
mahalanobis distances D^2
are computed with respect to the
covariance (location and scatter) estimate, and the weights are
1/sqrt(1+ pmax.int(0, 8*(D2  p)/sqrt(2*p)))
,
where D2 = D^2
and p = ncol(X)
.
Value
glmrob
returns an object of class "glmrob"
and is also
inheriting from glm
.
The summary
method, see summary.glmrob
, can
be used to obtain or print a summary of the results.
The generic accessor functions coefficients
,
effects
, fitted.values
and residuals
(see
residuals.glmrob
) can be used to extract various useful
features of the value returned by glmrob()
.
An object of class "glmrob"
is a list with at least the
following components:
coefficients 
a named vector of coefficients 
residuals 
the working residuals, that is the (robustly “huberized”) residuals in the final iteration of the IWLS fit. 
fitted.values 
the fitted mean values, obtained by transforming the linear predictors by the inverse of the link function. 
w.r 
robustness weights for each observations; i.e.,

w.x 
weights used to downweight observations based on the position of the observation in the design space. 
dispersion 
robust estimation of dispersion paramter if appropriate 
cov 
the estimated asymptotic covariance matrix of the estimated coefficients. 
tcc 
the tuning constant c in Huber's psifunction. 
family 
the 
linear.predictors 
the linear fit on link scale. 
deviance 
NULL; Exists because of compatipility reasons. 
iter 
the number of iterations used by the influence algorithm. 
converged 
logical. Was the IWLS algorithm judged to have converged? 
call 
the matched call. 
formula 
the formula supplied. 
terms 
the 
data 
the 
offset 
the offset vector used. 
control 
the value of the 
method 
the name of the robust fitter function used. 
contrasts 
(where relevant) the contrasts used. 
xlevels 
(where relevant) a record of the levels of the factors used in fitting. 
Author(s)
Andreas Ruckstuhl ("Mqle") and Martin Maechler
References
Eva Cantoni and Elvezio Ronchetti (2001) Robust Inference for Generalized Linear Models. JASA 96 (455), 1022–1030.
Eva Cantoni (2004) Analysis of Robust Quasideviances for Generalized Linear Models. Journal of Statistical Software, 10, http://www.jstatsoft.org/v10/i04
Eva Cantoni and Elvezio Ronchetti (2006) A robust approach for skewed and heavytailed outcomes in the analysis of health care expenditures. Journal of Health Economics 25, 198–213.
S. Heritier, E. Cantoni, S. Copt, M.P. VictoriaFeser (2009) Robust Methods in Biostatistics. Wiley Series in Probability and Statistics.
Marina Valdora and Víctor J. Yohai (2013) Robust estimators for Generalized Linear Models. In progress.
See Also
predict.glmrob
for prediction;
glmrobMqle.control
Examples
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  ## Binomial response 
data(carrots)
Cfit1 < glm(cbind(success, totalsuccess) ~ logdose + block,
data = carrots, family = binomial)
summary(Cfit1)
Rfit1 < glmrob(cbind(success, totalsuccess) ~ logdose + block,
family = binomial, data = carrots, method= "Mqle",
control= glmrobMqle.control(tcc=1.2))
summary(Rfit1)
Rfit2 < glmrob(success/total ~ logdose + block, weights = total,
family = binomial, data = carrots, method= "Mqle",
control= glmrobMqle.control(tcc=1.2))
coef(Rfit2) ## The same as Rfit1
## Binary response 
data(vaso)
Vfit1 < glm(Y ~ log(Volume) + log(Rate), family=binomial, data=vaso)
coef(Vfit1)
Vfit2 < glmrob(Y ~ log(Volume) + log(Rate), family=binomial, data=vaso,
method="Mqle", control = glmrobMqle.control(tcc=3.5))
coef(Vfit2) # c = 3.5 ==> not much different from classical
## Note the problems with tcc <= 3 %% FIXME algorithm ???
Vfit3 < glmrob(Y ~ log(Volume) + log(Rate), family=binomial, data=vaso,
method= "BY")
coef(Vfit3)## note that results differ much.
## That's not unreasonable however, see Kuensch et al.(1989), p.465
## Poisson response 
data(epilepsy)
Efit1 < glm(Ysum ~ Age10 + Base4*Trt, family=poisson, data=epilepsy)
summary(Efit1)
Efit2 < glmrob(Ysum ~ Age10 + Base4*Trt, family = poisson,
data = epilepsy, method= "Mqle",
control = glmrobMqle.control(tcc= 1.2))
summary(Efit2)
## 'x' weighting:
(Efit3 < glmrob(Ysum ~ Age10 + Base4*Trt, family = poisson,
data = epilepsy, method= "Mqle", weights.on.x = "hat",
control = glmrobMqle.control(tcc= 1.2)))
try( # gives singular cov matrix: 'Trt' is binary factor >
# affine equivariance and subsampling are problematic
Efit4 < glmrob(Ysum ~ Age10 + Base4*Trt, family = poisson,
data = epilepsy, method= "Mqle", weights.on.x = "covMcd",
control = glmrobMqle.control(tcc=1.2, maxit=100))
)
##> See example(possumDiv) for another Poissonregression
###  Gamma family  data from example(glm) 
clotting < data.frame(
u = c(5,10,15,20,30,40,60,80,100),
lot1 = c(118,58,42,35,27,25,21,19,18),
lot2 = c(69,35,26,21,18,16,13,12,12))
summary(cl < glm (lot1 ~ log(u), data=clotting, family=Gamma))
summary(ro < glmrob(lot1 ~ log(u), data=clotting, family=Gamma))
clotM5.high < within(clotting, { lot1[5] < 60 })
op < par(mfrow=2:1, mgp = c(1.6, 0.8, 0), mar = c(3,3:1))
plot( lot1 ~ log(u), data=clotM5.high)
plot(1/lot1 ~ log(u), data=clotM5.high)
par(op)
## Obviously, there the first observation is an outlier with respect to both
## representations!
cl5.high < glm (lot1 ~ log(u), data=clotM5.high, family=Gamma)
ro5.high < glmrob(lot1 ~ log(u), data=clotM5.high, family=Gamma)
with(ro5.high, cbind(w.x, w.r))## the 5th obs. is downweighted heavily!
plot(1/lot1 ~ log(u), data=clotM5.high)
abline(cl5.high, lty=2, col="red")
abline(ro5.high, lwd=2, col="blue") ## result is ok (but not "perfect")
