glmrob: 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 only for discrete response distributions, i.e. family = binomial or = poisson.

Usage

glmrob(formula, family, data, weights, subset, na.action,
       start = NULL, offset, method = "Mqle",
       weights.on.x = c("none", "hat", "robCov", "covMcd"), control = NULL,
       model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, trace = FALSE, ...)

Arguments

formula

a formula, i.e., a symbolic description of the model to be fit (cf. glm or lm).

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 function or the result of a call to a family function. (See

data

an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which glmrob is called.

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 NAs. The default is set by the na.action setting in options. The factory-fresh

start

starting values for the parameters in the linear predictor.

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

character string (can be abbreviated) specifying how points (potential outliers) in x-space are downweighted. If "hat", weights on the design of the form $\sqrt{1-h_{ii}}$ are used, where $h_{ii}$ are the diagonal elements of

control

a list of parameters for controlling the fitting process. See the documentation for glmrobMqle.control for details.

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 contrasts.arg of model.matrix.default.

trace

logical (or integer) indicating if intermediate results should be printed; defaults to FALSE.

...

arguments passed to glmrobMqle.control when control is NULL (as per default).

Value

glmrob returns an object of class "glmrob" and is also inheriting from glm. The function summary (i.e., summary.glmrob) can be used to obtain or print a summary of the results. The generic accessor functions coefficients, effects, fitted.values and residuals 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:
coefficientsa named vector of coefficients
residualsthe Pearson residuals
fitted.valuesthe fitted mean values, obtained by transforming the linear predictors by the inverse of the link function.
w.rrobustness weights for each observations; i.e., residuals $\times$ w.r equals the psi-function of the Preason's residuals.
w.xweights used to down-weight observations based on the position of the observation in the design space.
covthe estimated asymptotic covariance matrix of the estimated coefficients.
tccthe tuning constant c in Huber's psi-function.
familythe family object used.
linear.predictorsthe linear fit on link scale.
devianceup to a constant, minus twice the maximized log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero.
null.devianceThe deviance for the null model, comparable with deviance. The null model will include the offset, and an intercept if there is one in the model
iterthe number of iterations used by the influence algorithm.
convergedlogical. Was the IWLS algorithm judged to have converged?
callthe matched call.
formulathe formula supplied.
termsthe terms object used.
datathe data argument.
offsetthe offset vector used.
controlthe value of the control argument used.
methodthe 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.
If a binomial is specified by giving a two-column response ???

Details

method="Mqle" fits a generalized linear model using Mallows or Huber type robust estimators, as described in Cantoni and Ronchetti (2001). In contrast to the implementation described in Cantoni (2004), the pure influence algorithm is implemented. Currently no other method is implemented. weights.on.x= "robCov" makes sense if all explanatory variables are continuous.

References

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, http://www.jstatsoft.org

Examples

Run this code

## Binomial response --------------
data(carrots)

Cfit1 <- glm(cbind(success, total-success) ~ logdose + block,
             data = carrots, family = binomial)
summary(Cfit1)

Rfit1 <- glmrob(cbind(success, total-success) ~ 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))
## Note the problems with  tcc <= 3 %% FIXME algorithm ???
coef(Vfit2) # c = 3.5 ==> not much different from classical



## 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))
)

### -------- 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")

Run the code above in your browser using DataLab