Robust Multivariate Linear Models

The IRLS Algorithm

The IRLS algorithm for robust multivariate regression is shown in the figure below.

Flowchart for the iteratively reweighted least squares algorithm

The method proceeds as follows:

1. Initialization: Begin with an initial estimate of the regression coefficients, typically obtained from ordinary least squares (OLS).

2. Residual calculation: Compute the multivariate residuals for each observation: \[\mathbf{r}_i = \mathbf{y}_i - \mathbf{X}_i\hat{\boldsymbol{\beta}}\] where \(\mathbf{y}_i\) is the \(p \times 1\) response vector for observation \(i\), \(\mathbf{X}_i\) is the corresponding row of the design matrix, and \(\hat{\boldsymbol{\beta}}\) is the current estimate of the coefficient matrix.

3. Distance computation: Calculate the squared Mahalanobis distance of each residual vector from the origin: \[d_i^2 = \mathbf{r}_i^T \mathbf{S}^{-1} \mathbf{r}_i\] where \(\mathbf{S}\) is a robust estimate of the residual covariance matrix, computed using MASS::cov.trob().

4. Weight assignment: Assign weights to each observation based on their residual distances. Observations with larger distances receive smaller weights, effectively downweighting potential outliers: \[w_i = \rho(d_i^2)\] where \(\rho(\cdot)\) is a weight function that decreases as the distance increases.

5. Weighted Least Squares: Update the coefficient estimates using weighted least squares: \[\hat{\boldsymbol{\beta}}^{(new)} = (\mathbf{X}^T\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}\mathbf{Y}\] where \(\mathbf{W}\) is a diagonal matrix of weights.

6. Convergence check: Repeat steps 2-5 until convergence, typically assessed by monitoring changes in the coefficient estimates or weights between iterations.

Features of the Implementation

The robmlm() implementation incorporates several important features:

• Robust covariance estimation: The use of MASS::cov.trob() provides a robust estimate of the residual covariance matrix, which is crucial for computing meaningful Mahalanobis distances in the presence of outliers. This uses the multivariate \(t\) distribution that allows for longer tails. (There are other robust covariance estimators, such as Minimum Covariance Determinant (MCD) and Minimum Volume Ellipse (MVE), which have a high tolerance—breakdown-bound– for outliers. These might come to a future version of robmlm().)

• Inheritance Structure: The returned object inherits from both "mlm" and "lm" classes, ensuring compatibility with existing R infrastructure for linear models while adding robust-specific methods.

• Weight Preservation: The final weights are preserved in the fitted object, allowing users to identify influential observations and assess the impact of the robust fitting procedure.

• Diagnostic Capabilities: The plot.robmlm() method provides immediate visual feedback on the weighting scheme, plotting final weights against case numbers to highlight observations that were down-weighted during the fitting process.

• Visualization methods: Because a fitted robmlm object inherits the "mlm" class, it works directly with heplots’ existing diagnostic and hypothesis-test plots. distancePlot() plots the Mahalanobis distances of the model-matrix predictors against the distances of the residuals, so cases with both high leverage and a large residual — and so likely down-weighted during IRLS — stand out just as they would for a classical mlm. heplot() goes further: fitting the classical and robust models separately and overlaying the robust H and E ellipses on the classical model’s plot (add = TRUE) gives a direct geometric view of how down-weighting outliers reshapes the hypothesis and error variation — typically a visibly smaller E ellipse for the robust fit. See the example below.

Theoretical Properties

The IRLS approach offers several desirable theoretical properties:

• Breakdown Point: While not achieving the highest possible breakdown point, M-estimators like those implemented in IRLS can handle a reasonable proportion of outliers before completely breaking down.

• Influence Function: The bounded influence function of M-estimators ensures that no single observation can have unlimited impact on the final estimates.

• Asymptotic Efficiency: Under ideal conditions (no outliers, multivariate normality), robust M-estimators achieve high efficiency relative to classical least squares estimators.

• Equivariance: The robust estimates maintain appropriate equivariance properties under linear transformations of the data.

Classical MANOVA

We begin with the standard MANOVA model, and then examine some diagnostic plots.

# Classical MANOVA model
pottery.mlm <- lm(cbind(Al, Fe, Mg, Ca, Na) ~ Site, data = Pottery)
Anova(pottery.mlm)
#> 
#> Type II MANOVA Tests: Pillai test statistic
#>      Df test stat approx F num Df den Df  Pr(>F)    
#> Site  3      1.55      4.3     15     60 2.4e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Chisquare QQ plot: As an initial check, a cqplot() of the model gives a \(\chi^2\) QQ plot of the residuals from the model, which would identify badly fit observations. None seem particularly large here.

cqplot(pottery.mlm, id.n = 5)

Chisquare QQ plot for the Pottery model.

Influence plot: An influence plot for this model, done using mvinfluence::influencePlot() plots the leverage multivariate hat-values for the predictors (site) against squared multivariate studentized residuals, using the size of the point symbol proportional to a multivariate generalization of Cook’s D statistic.

res <- influencePlot(pottery.mlm, id.n = 2)

Influence plot for the Pottery model.

res |>
  arrange(desc(CookD))
#>          H        Q   CookD       L       R
#> 25 0.20000 0.621766 0.68394 0.25000 0.77721
#> 18 0.20000 0.355067 0.39057 0.25000 0.44383
#> 11 0.07143 0.587152 0.23067 0.07692 0.63232
#> 16 0.50000 0.006379 0.01754 1.00000 0.01276
#> 15 0.50000 0.006379 0.01754 1.00000 0.01276

Because site is a factor, the hat-values are inversely proportional to the sample size. Points for Llanedryn (\(n=14\)) are in the left-most column, followed by AshleyRails and IsleThorns (each with \(n=5\)) and then Caldicot (\(n=14\)).

Here, case 25 stands out with the largest value of Cook’s D, followed by 18 and 11.

Robust MANOVA

Now, fit the robust model using robmlm(). Because this uses an iterative IRLS method, points that might not seem unusual in the initial model can become more noteworthy when the extreme observations are down-weighted in a subsequent iteration.

# Robust MANOVA model
pottery.rlm <- robmlm(cbind(Al, Fe, Mg, Ca, Na) ~ Site, data = Pottery)
Anova(pottery.rlm)
#> 
#> Type II MANOVA Tests: Pillai test statistic
#>      Df test stat approx F num Df den Df  Pr(>F)    
#> Site  3      1.98     6.55     15     51 1.7e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s compare the results. From the results of Anova() above, you can see that the \(F\) statistic for the robust model is greater than that for OLS, which suggests that possible outliers may have reduced the strength of evidence for differences among the sites.

It is useful to compare the coefficients of the two models, and for this, their relative difference is a useful metric. The simple function reldiff() expresses these as the signed percent of difference between the classical and robust estimates.

b.mlm <- coef(pottery.mlm)
b.rlm <- coef(pottery.rlm)

reldiff <- function(x, y, pct=TRUE) {
  res <- abs(x - y) / x
  if (pct) res <- 100 * res
  res
}

reldiff(b.mlm, b.rlm)
#>                     Al      Fe      Mg      Ca      Na
#> (Intercept)      3.616  20.480  2.5908  10.333  0.6925
#> SiteCaldicot   -11.145   7.934  0.4832   2.211 16.6192
#> SiteIsleThorns  64.510 136.881 25.1469 -21.792  1.9605
#> SiteLlanedyrn  -11.559   9.969  2.0078   6.241 22.7443

Among these, the coefficients for IsleThorns are quite different for most of the variables.

Examining the Robust Weights

The robust fitting procedure assigns weights to each observation based on their deviation from the fitted model. Observations that appear to be outliers receive lower weights. The plot() method for a "roblm" object gives an index plot of the weight values.

# Plot the weights from robust fitting
plot(pottery.rlm, col=Pottery$Site, segments=TRUE)
xloc <- c(7.5, 15.5, 19.5, 24)
text(xloc, rep(c(1.0, 1.05), length=5), 
     levels(Pottery$Site), pos =3, xpd = TRUE)

Weights from robust MANOVA fitting showing potential outliers

The weight plot reveals which observations were considered potentially problematic during the robust fitting process. Observations with weights substantially below 1.0 were down-weighted, indicating they may be outliers in the multivariate space defined by the five chemical measurements.

Distance Plot

distancePlot() gives a complementary view of the same idea. It plots the Mahalanobis distances of the Site predictors (the model matrix X) against the distances of the residuals from the robust fit, labeling cases beyond the cutoff lines (dashed, level = 0.975 by default). This plot method, due to Rousseeuw & van Zomeren (1990, 1991) and Rousseeuw et al. (2004), helps you see WHY observations are downweighted, by combining information on regression outliers (unusual in the Y responses) and leverage points (unusual in the Xs).

distancePlot(pottery.rlm)

Distance plot for the robust Pottery model: Mahalanobis distances of the Site predictors vs. distances of the residuals.

Case 25 – already flagged with the largest Cook’s D in the classical influence plot above – again stands out, here with an unusually large residual distance. Cases 15 and 16 stand out instead on the predictor (Site) side, reflecting their group’s small sample size rather than an unusual response pattern.

Hypothesis-Error (HE) plots

HE plots provide a visual representation of the multivariate hypothesis test by showing the relationship between hypothesis and error variation in a reduced dimensional space. We’ll create HE plots for the first two variables (Al and Fe) and compare the classical and robust fits.

# Classical HE plot for Al and Fe
heplot(pottery.mlm, variables = c("Al", "Fe"), 
       main = "Classical vs Robust MANOVA: Al vs Fe",
       col = c("blue", "blue"), fill = TRUE, fill.alpha = 0.2)

# Overlay robust HE plot
heplot(pottery.rlm, variables = c("Al", "Fe"), 
       add = TRUE, error.ellipse = TRUE,
       col = c("red", "red"), 
       fill = TRUE, fill.alpha = 0.2, lty = 2)

# Add legend
legend("topright", 
       legend = c("Classical", "Robust"), 
       col = c("blue", "red"), 
       lty = c(1, 2), 
       fill = c("blue", "red")
       )

HE plot comparing classical (blue) and robust (red) MANOVA for Al vs Fe

The H and E ellipses have approximately the same shape and orientation for the classical and robust models, indicating that the pattern of differences among the group means is quite similar for both models. However, the E ellipse for the robust model is noticeably smaller. This goes into the greater \(F\) statistic for the robust model.

For a more complete interpretation of robust model, a scatterplot matrix of HE plots is shown below, using the pairs() method for a MLM.

pairs(pottery.rlm, 
      fill=TRUE, fill.alpha = 0.1)

Pairwise HE plots for all response variables in the robust model pottery.rlm

Interpretation and Discussion

Several important insights emerge from this analysis of the Pottery data

• Model evaluation: The comparison of MANOVA results shows …

• Outlier detection: The weight plot identifies specific observations that deviate from the typical pattern. In archaeological studies, such outliers might represent pottery from different time periods, contaminated samples, or measurement errors that warrant further investigation.

• Geometric interpretation: The HE plot comparison reveals how outliers affect the geometric representation of the hypothesis test. The hypothesis ellipse (representing the Site effect) and error ellipse (representing within-group variation) may differ substantially between classical and robust fits when influential outliers are present. The usual effect is that of seeing hypothesis effects more strongly when potential outliers are down-weighted.

• Practical implications: For archaeologists studying historical artifacts, the robust analysis provides more reliable conclusions about site or other factors differences by reducing the influence of potentially problematic observations. This is particularly important when making inferences about ancient trade routes, cultural practices, or chronological relationships based on chemical composition data. In some circumstances, the identification of such unusual observations can lead to insight about conventional theory.

The robust approach demonstrates its value by providing a more stable analysis that is less susceptible to the influence of outlying observations, while still maintaining high efficiency when the data conform to standard assumptions.

Classical and Robust Fits

To get started, fit the classical OLS multivariate regression and the robust version and run car::Anova() on each

pulp.mod <- lm(cbind(Y1, Y2, Y3, Y4) ~ X1 + X2 + X3 + X4, data = pulpfiber)
pulp.rlm <- robmlm(cbind(Y1, Y2, Y3, Y4) ~ X1 + X2 + X3 + X4, data = pulpfiber)

Anova(pulp.mod)
#> 
#> Type II MANOVA Tests: Pillai test statistic
#>    Df test stat approx F num Df den Df  Pr(>F)    
#> X1  1     0.157     2.51      4     54   0.053 .  
#> X2  1     0.404     9.16      4     54 1.0e-05 ***
#> X3  1     0.211     3.62      4     54   0.011 *  
#> X4  1     0.695    30.80      4     54 2.3e-13 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(pulp.rlm)
#> 
#> Type II MANOVA Tests: Pillai test statistic
#>    Df test stat approx F num Df den Df  Pr(>F)    
#> X1  1     0.357      6.9      4     50 0.00016 ***
#> X2  1     0.534     14.3      4     50 7.2e-08 ***
#> X3  1     0.368      7.3      4     50 0.00011 ***
#> X4  1     0.764     40.4      4     50 4.4e-15 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The contrast here is more dramatic than for Pottery:

• In the classical fit, X2 and X4 stand out clearly (\(p < .0001\)), X3 is only marginally significant (\(p = .011\)), and X1 falls just short of conventional significance (\(p = .053\)).
• In the robust fit, by contrast, all four fiber properties are highly significant predictors of the paper properties (\(p < .001\) in every case). This suggests that a handful of atypical samples in the classical fit are masking real predictive relationships involving X1 and X3.

Robust Weights

plot(pulp.rlm, segments = TRUE)

Weights from robust fitting of the pulp fiber data

Four cases – 51, 52, 56 and 61 – are driven to (essentially) zero weight; two more (22, 44) are markedly reduced (0.43 and 0.69); every other sample retains weight above 0.7. The weight plot alone doesn’t say why these six are atypical; for that we turn to the distance plot.

Distance Plot

A distancePlot() involves two independent choices: which residuals define the Y (vertical) axis – from pulp.mod (OLS) or pulp.rlm (robust) – and which covariance estimator computes the distances themselves, via method. In this example, only the second choice matters much: with the default method = "classical", the flagged cases are identical whether the residuals come from pulp.mod or pulp.rlm – the same eight cases (46, 51, 52, 56, 57, 58, 60, 61) either way. The main difference is seen in the covariance estimate, not the fit.

The reason is a masking effect: classical covariance for the X1-X4 fiber measurements is computed from all 62 samples, including the very outliers it is meant to help detect, so it is itself distorted by them and understates how unusual those samples really are. Rousseeuw et al. (2004) address exactly this for the pulp fiber data by using the MCD covariance estimator instead, which is resistant to the outliers it’s estimating around. The panels below compare the two, both using the robust-fit residuals from pulp.rlm:

op <- par(mfrow = c(1, 2))
distancePlot(pulp.rlm, method = "classical", main = "Classical distances")
distancePlot(pulp.rlm, method = "mcd", main = "MCD distances")

Distance plots for the robust pulp fiber model, comparing classical (left) and MCD-based (right) covariance estimates for the fiber (X) predictor distances.

par(op)

Under classical covariance (left), cases 60 and 62 – which turn out to have the two largest predictor distances of any sample once MCD is used – don’t stand out at all: their unusual fiber measurements are masked by a covariance matrix that treats them as ordinary. Under MCD (right), the picture sharpens considerably, flagging 21 cases rather than 8. Most of these retain weight above 0.8 in pulp.rlm – unusual in X, but well fit by the model – so we focus on the few that also matter for the robust fit:

• Cases 51, 52 and 56 have unremarkable fiber measurements (small distance on X) but very large residual distances: vertical outliers, poorly predicted by the fitted relationship regardless of their (ordinary) predictor values. These are exactly the three cases driven to zero weight above.

• Case 61 is extreme on both axes: a bad leverage point, unusual in its fiber properties and poorly fit by the model – correctly downweighted to zero along with the vertical outliers.

• Cases 60 and 62 have the two largest predictor distances of any sample, by a wide margin, yet small residual distances: good leverage points. They are unusual fiber samples, but ones the fitted relationship predicts well, so robmlm() correctly leaves them at near-full weight (0.98 each) rather than downweighting them.

• Case 22 sits in between, at weight 0.43: moderately unusual on both axes without being fully downweighted – a reminder that robmlm()’s weighting is continuous, not the all-or-nothing split the labels above suggest.

Multivariate Influence Plots

What we can learn from a distancePlot() can be compared to alternative displays that generalize the traditional influence plots of residuals vs. leverage, which are extended to the MLM case in the mvinfluence package.

influencePlot() for an "mlm" object draws a bubble plot, with circle areas proportional to a generalized Cook’s distance, using one of three types of axes. The default, type = "stres", plots (squared) Studentized residuals directly against the hat values (leverage). type = "LR" instead plots transformed residual and leverage components, R and L, chosen so that contours of constant Cook’s distance fall on diagonal lines – useful for seeing at a glance whether a case’s influence comes mainly from an unusual response, an unusual predictor pattern, or both.

influencePlot(pulp.mod, id.n = 5)

Default (type = 'stres') multivariate influence plot for the OLS pulp fiber model.

#>          H       Q  CookD       L       R
#> 46 0.24366 0.05983 0.1662 0.32216 0.07911
#> 51 0.06929 0.61222 0.4836 0.07445 0.65779
#> 52 0.10600 0.41406 0.5004 0.11857 0.46316
#> 56 0.12652 0.69124 0.9970 0.14485 0.79137
#> 57 0.20519 0.04397 0.1029 0.25817 0.05532
#> 58 0.21948 0.05280 0.1321 0.28119 0.06765
#> 59 0.18955 0.18764 0.4055 0.23389 0.23153
#> 60 0.49040 0.09928 0.5551 0.96231 0.19483
#> 61 0.41128 0.38188 1.7904 0.69859 0.64865

influencePlot(pulp.mod, id.n = 5, type = "LR")

type = 'LR' multivariate influence plot for the same model.

#>          H       Q  CookD       L       R
#> 46 0.24366 0.05983 0.1662 0.32216 0.07911
#> 51 0.06929 0.61222 0.4836 0.07445 0.65779
#> 52 0.10600 0.41406 0.5004 0.11857 0.46316
#> 56 0.12652 0.69124 0.9970 0.14485 0.79137
#> 57 0.20519 0.04397 0.1029 0.25817 0.05532
#> 58 0.21948 0.05280 0.1321 0.28119 0.06765
#> 59 0.18955 0.18764 0.4055 0.23389 0.23153
#> 60 0.49040 0.09928 0.5551 0.96231 0.19483
#> 61 0.41128 0.38188 1.7904 0.69859 0.64865

Both types flag the same nine cases here (46, 51, 52, 56, 57, 58, 59, 60, 61) – switching type reorganizes how the same influence information is displayed, not which cases stand out. Case 61 has the largest bubble in both (generalized Cook’s D = 1.79), consistent with its bad-leverage status above, and case 60 has by far the largest leverage component L, consistent with it being the more extreme of the two good-leverage cases. Case 62, though, doesn’t appear here at all: influencePlot() uses the ordinary (non-robust) hat values from pulp.mod, under which case 62’s leverage ranks only 11th of 62 – the same classical-covariance masking that hid it from the distance plot above.

Interpretation and Discussion

• Model evaluation: The gain from 2 to 4 significant predictors under the robust fit shows how a small number of outlying samples can mask genuine predictive relationships in classical multivariate regression – a stronger contrast than in the Pottery example.

• Outlier detection: The weight plot and distance plot together tell a more nuanced story than “outlier” vs. “not”: some flagged cases (51, 52, 56, 61) are poorly fit and rightly downweighted, while others (60, 62) are merely unusual in their predictors and rightly are not. Distinguishing good from bad leverage points this way is exactly the diagnostic problem distancePlot() is designed to address.

• Practical implications: For process data like this – where fiber properties are measured to predict paper quality – correctly identifying which atypical samples to downweight (vertical/bad-leverage cases) versus keep (good-leverage cases) matters for building a model that generalizes to new, unusual-but-plausible fiber batches, rather than one that has simply been fit to accommodate a few mismeasured samples.