heplots (version 1.3-8)

# robmlm: Robust Fitting of 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

```robmlm(X, ...)# S3 method for formula
robmlm(formula, data, subset, weights, na.action,
model = TRUE, contrasts = NULL, ...)# S3 method for default
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 robmlm
print(x, ...)# S3 method for summary.robmlm
print(x, ...)# S3 method for 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

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

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

## References

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

## See Also

`rlm`, `cov.trob`

## Examples

```# NOT RUN {
##############
# 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, ylab=vlab, 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="")

# }
```