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 mlms, and are therefore compatible with
other mlm extensions, including Anova and
heplot.
Usage
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 NAs 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.
See Also
rlm,
cov.trob
Examples
##############
# 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="")