08 — Robustness

Which conclusions survive when assumptions are relaxed?

Scientific question

In Chapter 07, the assumptions of our interaction model appeared reasonable.

That is useful, but it is not the end of the story.

Real analytical datasets can contain heteroscedasticity, unusual observations and non-normal residuals. Sometimes these departures are mild and have little practical consequence. Sometimes they substantially alter the uncertainty attached to an estimate.

The relevant question is therefore not:

“Can we find a method that does not require assumptions?”

Every statistical method makes assumptions.

The useful question is:

“How sensitive are our conclusions to the assumptions of the method we chose?”

This chapter explores that question by deliberately stressing the dataset and comparing different inferential approaches.

The purpose is not to replace the original analysis with a collection of alternative tests.

It is to understand what changes when the assumptions change.

The dataset

We return to the balanced dataset.

The original data remain untouched.

For the robustness exercises, we create separate modified datasets so that every analysis remains reproducible and the original observations are never overwritten.

library(here)
library(data.table)

dat <- fread(
  here("data", "dataset_v01_clean_balanced.csv")
)

dat[, Method := factor(
  Method,
  levels = c("Reference", "Method_A", "Method_B", "Method_C")
)]

dat[, Matrix := factor(
  Matrix,
  levels = c("Fill", "Soil", "Sediment")
)]

Robustness is not the same as robustness to everything

Before comparing methods, it is useful to be precise about what we mean by “robust”.

A confidence interval can be robust to moderate non-normality but not to dependence between observations.

A heteroscedasticity-consistent standard error can protect an inference against unequal variances, but it does not solve a badly misspecified mean structure.

A bootstrap can reduce reliance on a parametric sampling distribution, but it still depends on the resampling scheme being appropriate.

A permutation test can provide a useful alternative under certain null hypotheses, but exchangeability is itself an assumption.

There is therefore no single concept of “robustness”.

We always need to ask:

Robust to which violation, for which estimate or hypothesis, under which design?

Unequal variances: Welch’s t-test

A simple case is the comparison between Method B and Method C.

The classical two-sample t-test assumes a common variance:

\[ \sigma_B^2 = \sigma_C^2 \]

Welch’s t-test relaxes this assumption.

Instead of pooling the two sample variances into a common estimate, it uses the individual variances and modifies the degrees of freedom accordingly.

We can compare the two approaches on the original data.

bc <- dat[
  Method %in% c("Method_B", "Method_C")
]

t_classical <- t.test(
  RelativeBias_pct ~ Method,
  data = bc,
  var.equal = TRUE
)

t_welch <- t.test(
  RelativeBias_pct ~ Method,
  data = bc,
  var.equal = FALSE
)

cat("Classical t-test:\n")

Classical t-test:

print(t_classical)


    Two Sample t-test

data:  RelativeBias_pct by Method
t = 6.696, df = 118, p-value = 7.652e-10
alternative hypothesis: true difference in means between group Method_B and group Method_C is not equal to 0
95 percent confidence interval:
 3.206907 5.900260
sample estimates:
mean in group Method_B mean in group Method_C 
              1.952717              -2.600867

cat("\nWelch's t-test:\n")


Welch's t-test:

print(t_welch)


    Welch Two Sample t-test

data:  RelativeBias_pct by Method
t = 6.696, df = 117.91, p-value = 7.669e-10
alternative hypothesis: true difference in means between group Method_B and group Method_C is not equal to 0
95 percent confidence interval:
 3.206896 5.900270
sample estimates:
mean in group Method_B mean in group Method_C 
              1.952717              -2.600867

For the present data, the two analyses give very similar results.

That is what we would expect: there is no strong evidence that unequal variances are driving the B–C comparison.

The lesson is not that Welch’s test is “better”.

It is that when the equal-variance assumption is questionable, an inferential procedure that does not require it may be preferable.

A deliberately heteroscedastic dataset

We can make the issue more visible by creating a deliberately modified dataset.

Here we increase the dispersion of Method B while keeping its mean approximately unchanged.

The purpose is pedagogical: we are not claiming that this is a realistic model of a particular analytical failure. We are creating a controlled stress test.

dat_hetero <- copy(dat)

idx_b <- which(dat_hetero$Method == "Method_B")

mu_b <- mean(dat_hetero$RelativeBias_pct[idx_b])
sd_b <- sd(dat_hetero$RelativeBias_pct[idx_b])

dat_hetero[idx_b, RelativeBias_pct :=  mu_b +
  3 * (RelativeBias_pct - mu_b)]

We can inspect the group-specific standard deviations.

dat_hetero[
  ,
  .(
    Mean = mean(RelativeBias_pct),
    SD = sd(RelativeBias_pct),
    N = .N
  ),
  by = Method
]

      Method       Mean        SD     N
      <fctr>      <num>     <num> <int>
1: Reference  -9.522000  3.226779    60
2:  Method_A -17.627167  4.060234    60
3:  Method_B   1.952717 11.325641    60
4:  Method_C  -2.600867  3.673626    60

The variance structure has now been deliberately altered.

This gives us a useful stress test: does our inference depend strongly on the equal-variance assumption?

Classical versus Welch

bc_hetero <- dat_hetero[
  Method %in% c("Method_B", "Method_C")
]

t_classical_hetero <- t.test(
  RelativeBias_pct ~ Method,
  data = bc_hetero,
  var.equal = TRUE
)

t_welch_hetero <- t.test(
  RelativeBias_pct ~ Method,
  data = bc_hetero,
  var.equal = FALSE
)

cat("Classical t-test:\n")

Classical t-test:

print(t_classical_hetero)


    Two Sample t-test

data:  RelativeBias_pct by Method
t = 2.9624, df = 118, p-value = 0.003692
alternative hypothesis: true difference in means between group Method_B and group Method_C is not equal to 0
95 percent confidence interval:
 1.509652 7.597515
sample estimates:
mean in group Method_B mean in group Method_C 
              1.952717              -2.600867

cat("\nWelch's t-test:\n")


Welch's t-test:

print(t_welch_hetero)


    Welch Two Sample t-test

data:  RelativeBias_pct by Method
t = 2.9624, df = 71.279, p-value = 0.004146
alternative hypothesis: true difference in means between group Method_B and group Method_C is not equal to 0
95 percent confidence interval:
 1.488846 7.618321
sample estimates:
mean in group Method_B mean in group Method_C 
              1.952717              -2.600867

The exact numerical difference depends on the data-generating perturbation, but the principle is more important than the particular p-value.

The classical test answers the comparison under a common-variance model.

Welch’s test allows the two groups to have different variances.

If the scientific question is simply whether the group means differ, and unequal variances are plausible, Welch’s procedure is often a safer default for a two-group comparison.

Robust standard errors in the linear model

The same principle extends to the linear model.

Suppose the mean structure is correctly specified but the residual variance is not constant.

The ordinary least-squares coefficient estimates can still be useful, but their conventional standard errors may be unreliable.

One approach is to retain the model coefficients while replacing the classical covariance estimator with a heteroscedasticity-consistent estimator.

We use HC3, which is a commonly used small-sample adjustment.

library(sandwich)
library(lmtest)

fit_classical <- lm(
  RelativeBias_pct ~ Method,
  data = dat
)

cat("Classical standard errors:\n")

Classical standard errors:

print(coeftest(fit_classical))


t test of coefficients:

               Estimate Std. Error t value  Pr(>|t|)    
(Intercept)    -9.52200    0.47717 -19.955 < 2.2e-16 ***
MethodMethod_A -8.10517    0.67482 -12.011 < 2.2e-16 ***
MethodMethod_B 11.47472    0.67482  17.004 < 2.2e-16 ***
MethodMethod_C  6.92113    0.67482  10.256 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nHC3 robust standard errors:\n")


HC3 robust standard errors:

print(
  coeftest(
    fit_classical,
    vcov = vcovHC(fit_classical, type = "HC3")
  )
)


t test of coefficients:

               Estimate Std. Error t value  Pr(>|t|)    
(Intercept)    -9.52200    0.42009 -22.666 < 2.2e-16 ***
MethodMethod_A -8.10517    0.67520 -12.004 < 2.2e-16 ***
MethodMethod_B 11.47472    0.64656  17.747 < 2.2e-16 ***
MethodMethod_C  6.92113    0.63656  10.873 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The important point is what has and has not changed.

The fitted coefficients have not changed.

The model’s estimated mean structure has not changed.

What changes is our estimate of the uncertainty associated with those coefficients.

This is fundamentally different from transforming the response or fitting a different mean model.

Robust standard errors therefore address a relatively specific problem:

The mean model is retained, but the classical variance calculation is relaxed.

They do not make an incorrect model correct.

Bootstrap: relaxing the parametric sampling distribution

Another strategy is the bootstrap.

Suppose we are interested in the B − C difference in mean relative bias:

\[ \Delta = \mu_B - \mu_C \]

The classical confidence interval is derived from the sampling distribution implied by the linear model.

The bootstrap instead approximates the sampling distribution by repeatedly resampling observations from the observed data.

For this comparison, the appropriate resampling scheme is to resample within Method.

We preserve the group structure while allowing the observations within each group to vary.

set.seed(42)

b_values <- dat[
  Method == "Method_B",
  RelativeBias_pct
]

c_values <- dat[
  Method == "Method_C",
  RelativeBias_pct
]

observed_diff <-
  mean(b_values) -
  mean(c_values)

n_boot <- 10000

boot_diff <- numeric(n_boot)

for (i in seq_len(n_boot)) {

  boot_b <- sample(
    b_values,
    length(b_values),
    replace = TRUE
  )

  boot_c <- sample(
    c_values,
    length(c_values),
    replace = TRUE
  )

  boot_diff[i] <-
    mean(boot_b) -
    mean(boot_c)
}

ci_boot <- quantile(
  boot_diff,
  c(0.025, 0.975)
)

cat(
  "Observed B - C difference:",
  round(observed_diff, 3),
  "percentage points\n"
)

Observed B - C difference: 4.554 percentage points

cat(
  "Bootstrap 95% CI:",
  round(ci_boot[1], 3),
  "to",
  round(ci_boot[2], 3),
  "percentage points\n"
)

Bootstrap 95% CI: 3.217 to 5.879 percentage points

The bootstrap interval is expected to be close to the classical interval for this dataset.

That agreement is useful evidence that the B − C estimate is not particularly sensitive to the normal-theory calculation.

But the bootstrap is not assumption-free.

It assumes that the observed sample is a reasonable representation of the population from which we resample.

It also assumes that our resampling scheme respects the experimental structure.

If the observations are correlated because several measurements come from the same sample, independently resampling individual observations would not reproduce the relevant sampling process.

The bootstrap is therefore not magic.

It moves the assumptions to a different place.

A deliberately extreme observation

We now examine a different kind of violation.

Suppose one measurement for Method B is unexpectedly extreme.

This is deliberately constructed as a single observation, rather than shifting an entire group.

dat_outlier <- copy(dat)

set.seed(42)

candidate_idx <- which(
  dat_outlier$Method == "Method_B" &
  dat_outlier$Matrix == "Fill"
)

outlier_idx <- sample(candidate_idx, 1)

dat_outlier[
  outlier_idx,
  RelativeBias_pct := RelativeBias_pct + 50
]

cat(
  "Modified observation:",
  outlier_idx,
  "\n"
)

Modified observation: 67

We can compare the original and contaminated Method B means.

comparison <- rbind(
  dat[
    ,
    .(
      Mean = mean(RelativeBias_pct),
      SD = sd(RelativeBias_pct)
    ),
    by = Method
  ][
    Method == "Method_B"
  ][
    ,
    Dataset := "Original"
  ],
  dat_outlier[
    ,
    .(
      Mean = mean(RelativeBias_pct),
      SD = sd(RelativeBias_pct)
    ),
    by = Method
  ][
    Method == "Method_B"
  ][
    ,
    Dataset := "With one extreme observation"
  ]
)

comparison[]

     Method     Mean       SD                      Dataset
     <fctr>    <num>    <num>                       <char>
1: Method_B 1.952717 3.775214                     Original
2: Method_B 2.786050 7.596881 With one extreme observation

A single observation can have a surprisingly large effect on a sample mean, especially when the quantity of interest is itself a mean.

This is not necessarily evidence that the observation should be removed.

The first question is scientific:

Is the observation plausible?

If it corresponds to a transcription error, instrument malfunction or demonstrable analytical failure, correction or exclusion may be justified according to a predefined protocol.

If it is a legitimate but unusual sample, removing it simply because it changes the answer would introduce bias.

Does the conclusion change?

We can compare the one-way models fitted to the original and contaminated data.

fit_clean <- lm(
  RelativeBias_pct ~ Method,
  data = dat
)

fit_outlier <- lm(
  RelativeBias_pct ~ Method,
  data = dat_outlier
)

cat("Original model:\n")

Original model:

print(anova(fit_clean))

Analysis of Variance Table

Response: RelativeBias_pct
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Method      3 13127.4  4375.8  320.31 < 2.2e-16 ***
Residuals 236  3224.1    13.7                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nModel with one extreme observation:\n")


Model with one extreme observation:

print(anova(fit_outlier))

Analysis of Variance Table

Response: RelativeBias_pct
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Method      3 14048.9  4683.0  190.94 < 2.2e-16 ***
Residuals 236  5788.2    24.5                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The global conclusion that the methods differ may remain unchanged.

But the estimated Method B mean and some of the pairwise differences can move substantially.

This distinction is important.

A global hypothesis test can remain significant while the estimated magnitude of a particular effect changes materially.

Therefore:

Robustness of a p-value is not the same thing as robustness of an effect estimate.

In analytical chemistry, the second question can be more important than the first.

If a method difference moves from 4% to 7% depending on one observation, the practical interpretation may change even if the global ANOVA remains highly significant.

Permutation test: an alternative when distributional assumptions are questionable

Permutation tests provide another way to construct a null distribution.

For the B − C comparison, the null hypothesis is that the two groups are exchangeable with respect to Method.

Under that null, the Method labels can be permuted while keeping the observed response values fixed.

We can therefore compare the observed mean difference with the distribution obtained after repeatedly shuffling the labels.

set.seed(42)

bc <- dat[
  Method %in% c("Method_B", "Method_C")
]

y <- bc$RelativeBias_pct
group <- bc$Method

observed_diff <- mean(
  y[group == "Method_B"]
) - mean(
  y[group == "Method_C"]
)

n_perm <- 10000

perm_diff <- numeric(n_perm)

for (i in seq_len(n_perm)) {

  perm_group <- sample(group)

  perm_diff[i] <- mean(
    y[perm_group == "Method_B"]
  ) - mean(
    y[perm_group == "Method_C"]
  )
}

p_perm <- mean(
  abs(perm_diff) >= abs(observed_diff)
)

cat(
  "Observed B - C difference:",
  round(observed_diff, 3),
  "percentage points\n"
)

Observed B - C difference: 4.554 percentage points

cat(
  "Permutation p-value:",
  format.pval(p_perm),
  "\n"
)

Permutation p-value: < 2.22e-16

The permutation distribution provides an empirical reference distribution for the contrast under the null hypothesis.

This can be attractive when the normal approximation is questionable.

But again, there is an assumption hiding underneath:

The observations must be exchangeable under the null hypothesis.

If the observations have a hierarchical structure, temporal dependence, paired measurements or other constraints, arbitrary permutation of individual observations can be invalid.

The correct permutation scheme must reproduce the structure implied by the null hypothesis.

Permutation tests therefore do not eliminate assumptions.

They make a different set of assumptions.

Comparing the approaches

We have now used several methods to interrogate the same scientific problem.

Approach	Main concern addressed	What changes?
Classical t-test	Standard parametric inference	Assumes equal variances
Welch t-test	Unequal variances Standard error and degrees of freedom
HC3 standard errors	Heteroscedasticity Estimated uncertainty
Bootstrap	Parametric sampling distribution	Sampling distribution of the estimate
Permutation test	Null distribution under exchangeability	Reference distribution
Influence analysis	Sensitivity to individual observations	Stability of estimates

These approaches should not be thought of as interchangeable competitors.

Each one answers a somewhat different methodological concern.

What did we learn from the stress tests?

The original dataset was deliberately well behaved, to gives us a reference point.

When we perturb the data, however, different aspects of the analysis respond differently.

A variance violation primarily affects uncertainty estimates.

An extreme observation can substantially affect estimated means and contrasts.

Non-normality may have little practical effect in one design and a major effect in another.

And none of these procedures can compensate for a violation of independence caused by an incorrect experimental-unit definition.

This is why robustness analysis should be targeted.

There is little value in running five alternative methods simply to see whether they all produce \(p < 0.05\).

A more informative workflow is:

identify a plausible failure mode;
determine which inferential quantity it threatens;
choose an appropriate alternative or sensitivity analysis;
compare the resulting estimates and uncertainty;
decide whether the scientific conclusion changes.

Statistical significance versus practical stability

Our five percentage-point threshold provides a useful perspective.

Suppose a method difference is estimated as:

\[ \hat{\Delta} = 4.6\% \]

with a narrow confidence interval.

That result may be statistically very convincing, but it remains below our predefined threshold for chemical relevance.

Now suppose a single influential observation changes the estimate to 7%.

The issue is no longer simply whether the p-value remains below 0.05.

The practical conclusion has changed.

This is why sensitivity analysis should report effect estimates and uncertainty, not only significance decisions.

For analytical chemistry, the question is often:

Would the conclusion about method suitability change under a reasonable alternative analysis?

That is a much more useful definition of robustness than “the p-value stayed significant”.

What robustness cannot fix

There are limits to all of these approaches.

None of them can compensate for:

a wrong experimental unit;
pseudoreplication;
unrecognised dependence between measurements;
systematic measurement bias;
a missing scientifically important predictor;
an incorrect functional form;
confounding that the design cannot distinguish.

If aliquots from the same physical sample are treated as independent when the scientific unit is actually the sample, no robust standard error or bootstrap applied at the wrong level can rescue the analysis.

The resampling or variance estimation procedure must respect the structure of the experiment.

This is one reason why statistical robustness ultimately depends on understanding the analytical measurement process.

Take-home message

Robustness is not about finding a different answer. It is about understanding how much the answer depends on the assumptions behind the analysis.

A good robustness analysis does not produce a list of alternative p-values.

It identifies plausible ways the model could fail and asks whether those failures materially change the scientific conclusion.

For our dataset, the original inference is reassuringly stable under several reasonable checks.

But the most important lesson is broader: robust methods do not remove the need to understand the data-generating process.

They simply make some parts of our inference less dependent on particular assumptions.

And when the experimental structure itself changes — for example, when we have several correlated responses measured on the same sample — the problem is no longer just about making a univariate ANOVA more robust.

It becomes a different statistical question.

Next: MANOVA

So far, we have analysed one response variable at a time: relative bias.

But analytical chemistry rarely produces only one useful measurement.

A method may have a systematic bias and a precision characteristic that matter simultaneously. These responses can also be correlated.

Analysing them separately creates a new problem: multiple hypothesis tests and loss of information about their joint structure.

The next question is therefore:

“Do the analytical methods differ when bias and precision are considered simultaneously?”

This leads us to MANOVA.