robmlm: Robust Fitting of Multivariate Linear Models In heplots: Visualizing Hypothesis Tests in Multivariate Linear Models

Description

Fit a multivariate linear model by robust regression using a simple M estimator.

These S3 methods are designed to provide a specification of a class of robust methods which extend `mlm`s, and are therefore compatible with other `mlm` extensions, including `Anova` and `heplot`.

Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20``` ```robmlm(X, ...) ## S3 method for class 'formula' robmlm(formula, data, subset, weights, na.action, model = TRUE, contrasts = NULL, ...) ## Default S3 method: robmlm(X, Y, w, P = 2 * pnorm(4.685, lower.tail = FALSE), tune, max.iter = 100, psi = psi.bisquare, tol = 1e-06, initialize, verbose = FALSE, ...) ## S3 method for class 'robmlm' print(x, ...) ## S3 method for class 'summary.robmlm' print(x, ...) ## S3 method for class 'robmlm' summary(object, ...) ```

Arguments

 `formula` a formula of the form `cbind(y1, y2, ...) ~ x1 + x2 + ...`. `data` a data frame from which variables specified in `formula` are preferentially to be taken. `subset` An index vector specifying the cases to be used in fitting. `weights` a vector of prior weights for each case. `na.action` A function to specify the action to be taken if `NA`s are found. The 'factory-fresh' default action in R is `na.omit`, and can be changed by `options``(na.action=)`. `model` should the model frame be returned in the object? `contrasts` optional contrast specifications; see `lm` for details. `...` other arguments, passed down. In particular relevant control arguments can be passed to the to the `robmlm.default` method. `X` for the default method, a model matrix, including the constant (if present) `Y` for the default method, a response matrix `w` prior weights `P` two-tail probability, to find cutoff quantile for chisq (tuning constant); default is set for bisquare weight function `tune` tuning constant (if given directly) `max.iter` maximum number of iterations `psi` robustness weight function; `psi.bisquare` is the default `tol` convergence tolerance, maximum relative change in coefficients `initialize` modeling function to find start values for coefficients, equation-by-equation; if absent WLS (`lm.wfit`) is used `verbose` show iteration history? (`TRUE` or `FALSE`) `x` a `robmlm` object `object` a `robmlm` object

Details

Fitting is done by iterated re-weighted least squares (IWLS), using weights based on the Mahalanobis squared distances of the current residuals from the origin, and a scaling (covariance) matrix calculated by `cov.trob`. The design of these methods were loosely modeled on `rlm`.

An internal `vcov.mlm` function is an extension of the standard `vcov` method providing for observation weights.

Value

An object of class `"robmlm"` inheriting from `c("mlm", "lm")`.

This means that the returned `"robmlm"` contains all the components of `"mlm"` objects described for `lm`, plus the following:

 `weights ` final observation weights `iterations ` number of iterations `converged ` logical: did the IWLS process converge?

The generic accessor functions `coefficients`, `effects`, `fitted.values` and `residuals` extract various useful features of the value returned by `robmlm`.

Author(s)

John Fox; packaged by Michael Friendly

References

A. Marazzi (1993) Algorithms, Routines and S Functions for Robust Statistics. Wadsworth & Brooks/Cole.

`rlm`, `cov.trob`
 ``` 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``` ```############## # Skulls data # make shorter labels for epochs and nicer variable labels in heplots Skulls\$epoch <- factor(Skulls\$epoch, labels=sub("c","",levels(Skulls\$epoch))) # variable labels vlab <- c("maxBreadth", "basibHeight", "basialLength", "nasalHeight") # fit manova model, classically and robustly sk.mod <- lm(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls) sk.rmod <- robmlm(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls) # standard mlm methods apply here coefficients(sk.rmod) # index plot of weights plot(sk.rmod\$weights, type="h", xlab="Case Index", ylab="Robust mlm weight", col="gray") points(sk.rmod\$weights, pch=16, col=Skulls\$epoch) axis(side=1, at=15+seq(0,120,30), labels=levels(Skulls\$epoch), tick=FALSE, cex.axis=1) # heplots to see effect of robmlm vs. mlm heplot(sk.mod, hypotheses=list(Lin="epoch.L", Quad="epoch.Q"), xlab=vlab[1], ylab=vlab[2], cex=1.25, lty=1) heplot(sk.rmod, hypotheses=list(Lin="epoch.L", Quad="epoch.Q"), add=TRUE, error.ellipse=TRUE, lwd=c(2,2), lty=c(2,2), term.labels=FALSE, hyp.labels=FALSE, err.label="") ############## # Pottery data pottery.mod <- lm(cbind(Al,Fe,Mg,Ca,Na)~Site, data=Pottery) pottery.rmod <- robmlm(cbind(Al,Fe,Mg,Ca,Na)~Site, data=Pottery) Anova(pottery.mod) Anova(pottery.rmod) # index plot of weights plot(pottery.rmod\$weights, type="h") points(pottery.rmod\$weights, pch=16, col=Pottery\$Site) # heplots to see effect of robmlm vs. mlm heplot(pottery.mod, cex=1.3, lty=1) heplot(pottery.rmod, add=TRUE, error.ellipse=TRUE, lwd=c(2,2), lty=c(2,2), term.labels=FALSE, err.label="") ############### # Prestige data # treat women and prestige as response variables for this example prestige.mod <- lm(cbind(women, prestige) ~ income + education + type, data=Prestige) prestige.rmod <- robmlm(cbind(women, prestige) ~ income + education + type, data=Prestige) coef(prestige.mod) coef(prestige.rmod) # how much do coefficients change? round(coef(prestige.mod) - coef(prestige.rmod),3) # pretty plot of case weights plot(prestige.rmod\$weights, type="h", xlab="Case Index", ylab="Robust mlm weight", col="gray") points(prestige.rmod\$weights, pch=16, col=Prestige\$type) legend(0, 0.7, levels(Prestige\$type), pch=16, col=palette()[1:3], bg="white") heplot(prestige.mod, cex=1.4, lty=1) heplot(prestige.rmod, add=TRUE, error.ellipse=TRUE, lwd=c(2,2), lty=c(2,2), term.labels=FALSE, err.label="") ```