Assessing the Linear Modelling Assumption of Homoscedasticity



Julia Piaskowski

March 3, 2026



“Residuals were inspected for normality and common variance visually using Q–Q plots and plots of residuals versus predicted values.”

Linear Model Review

\[Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\]

\(Y_i\) = dependent variable (there is only 1)

\(X_i\) = independent variable(s) (there may be many)

\(B_0\) = model intercept, the overall mean of \(Y\)

\(B_1\) = how \(Y\) changes with \(X\)

\(\epsilon_i\) = model residual, the gap between the predicted value for \(Y_i\) and its observed value

The Residual Is Important

\[ \epsilon_i \sim N(0, \sigma)\]

  • The residuals are normally distributed with a mean of zero and some standard deviation that we will estimate during the model-fitting process.

  • All model residuals are independent and identically distributed, “i.i.d.” – they are uncorrelated with each other and drawn from the same distribution (that has exactly one variance)

More on iid

  • For every observation, there is a residual: \(\{n_1, n_2,...,n_j\}\), \(\{e_1, e_2,...,e_j\}\)

  • imagine we have a data set of 3 observations and hence 3 residuals: \(\{e_1, e_2, e_3\}\)

  • there is a variance-covariance matrix for the residuals:

\[ \epsilon_i \sim N (0, \begin{bmatrix} \sigma & 0 & 0 \\ 0 & \sigma & 0 \\ 0 & 0 & \sigma \end{bmatrix} ) \]

  • all diagonal elements share a common variance (identical distribution), the off-diagonals are zero (independently distributed)

Why is iid Important?

  • ANOVA F-tests rely on this
  • The standard errors of the estimates and pairwise comparisons rely on this assumption
  • unequal error violations can result in higher type I error rate (false positives)
  • independence violations can results in a higher type II error rate (false negatives)

Example Analysis 1: Assumptions Met

  • White pine study, 4 pollen parents crossed with 7 female trees (4 reps), epicotyl length was assessed
data("hanover.whitepine", package = "agridat")
hist(hanover.whitepine$length, main = NULL, xlab = "epicotyl length", cex.lab = cex_lab, col = "#65B362")

Example Analysis 1: Assumptions Met

library(lme4)
m1 <- lmer(length ~ male*female + (1|rep), data = hanover.whitepine)
plot(m1, xlab = "predicted values", ylab = "Pearson residuals",  main = "Standard Residual Plot")

Example Analysis 1: Assumptions Met

plot(rstudent(m1), col = "blue3", main = "Plot by Data Set Order", ylab = "studentized residuals", xlab = "spreadsheet row"); abline(h = 0)

Example Analysis 1: Assumptions Met

plot(m1, resid(., scaled=TRUE) ~fitted(.) | rep, xlab = "predicted values", ylab = "Pearson residuals",  main = "Standard Residual Plot by Rep", abline = 0)

Example Analysis 1: Assumptions Met

plot(m1, male ~ resid(.,type = "response", scaled=TRUE), main = "Residuals by a Variable", xlab = "standardized Residuals")

Example Analysis 1: Assumptions Met

plot(m1, female ~ resid(.,type = "response", scaled=TRUE), main = "Residuals by a Variable", xlab = "standardized Residuals")

What if the iid assumptions are not met?

In the “general” linear model days, when a premium was placed on the i.i.d. error paradigm, if we did reject \(H_0\), it would set off a minor crisis in the form of a hunt for a variance stabilizing transformation. In contemporary modeling, we simply proceed with inference on estimable functions using the unequal variance model.


– Walt Stroup, Marina Ptukhiuna and Julie Garai (2024), Generalized Linear Mixed Models, 2nd Ed, Section 7.2.3.1

Assumptions Not Met

Heteroscedasticity [1]

  • When an independent factor in the study is causing heteroscedasticity
  • multi-environment agronomic study: 38 plant varieties, 9 locations, 4 reps each (RCBD)
var_ex1 <- read.csv(here::here("static", "tutorials", "data", "MET_trial_variance.csv")) |> 
  mutate(block = as.character(block))
head(var_ex1, n = 4)
    site block    variety    yield
1 site_7     4 variety_23       NA
2 site_7     1 variety_22 2.881003
3 site_7     2 variety_23 4.321505
4 site_7     4  variety_3 4.321505

Assumptions Not Met

Heteroscedasticity [1]

hist(var_ex1$yield, main=NA, col = "orange", ylab=NULL, xlab = "seed yield", cex.lab = cex_lab, cex.axis = cex_axis, breaks = 15)

Assumptions Not Met

Heteroscedasticity [1]

boxplot(yield ~ site, data = var_ex1, col = "orange", xlab = "location")

Assumptions Not Met

Heteroscedasticity [1]

library(nlme)
m2 <- lme(yield ~ site:variety + variety, random = ~ 1 |site/block, 
          na.action = na.exclude, data = var_ex1)
plot(m2)

Assumptions Not Met

Heteroscedasticity [1]

plot(m2, site  ~ resid(.,type = "response", scaled=TRUE))

Assumptions Not Met

Heteroscedasticity [1]

Solution: model variance for each site with varIdent(), an ‘nlme’ function

m2_new <- lme(yield ~ site:variety + variety, random = ~ 1 |site/block, 
              na.action = na.exclude, data = var_ex1, 
              weights = varIdent(form = ~1|site))
plot(m2_new)

Assumptions Not Met

Heteroscedasticity [1]

library(emmeans)
emmeans(m2, ~ variety*site, at = list(variety = "variety_1", site = c("site_1", "site_8", "site_7")))
 variety   site   emmean  SE df lower.CL upper.CL
 variety_1 site_1   1713 126  8     1423     2002
 variety_1 site_8   1327 126  8     1038     1617
 variety_1 site_7    147 126  8     -142      437

Degrees-of-freedom method: containment 
Confidence level used: 0.95 
emmeans(m2_new, ~ variety*site, at = list(variety = "variety_1", site = c("site_1", "site_8", "site_7")))
 variety   site   emmean    SE df lower.CL upper.CL
 variety_1 site_1   1713 186.0  8   1283.8     2142
 variety_1 site_8   1327 124.0  8   1040.6     1614
 variety_1 site_7    147  98.8  8    -80.5      375

Degrees-of-freedom method: containment 
Confidence level used: 0.95

Assumptions Not Met

Heteroscedasticity [1]

anova(m2, type = "marginal") 
             numDF denDF   F-value p-value
(Intercept)      1   977 185.95512  <.0001
variety         37   977  28.28370  <.0001
site:variety   304   977   7.80575  <.0001
anova(m2_new, type = "marginal") 
             numDF denDF  F-value p-value
(Intercept)      1   977 84.79352  <.0001
variety         37   977 10.75473  <.0001
site:variety   304   977 10.19020  <.0001

Assumptions Not Met

Heteroscedasticity [2]

  • When the mean and variance are linked
  • single agronomic study: 38 plant varieties evaluated in one trial, 4 reps each (RCBD)
var_ex2 <- read.csv(here::here("static", "tutorials", "data", "single_trial_variance.csv")) |> 
  dplyr::mutate(block = as.character(block))

head(var_ex2, n = 4)
  block    variety    yield
1     1  variety_1 127.4420
2     1 variety_10 528.7488
3     1 variety_11 409.4414
4     1 variety_12 173.5381

Assumptions Not Met

Heteroscedasticity [2]

hist(var_ex2$yield, col = "#F88379", main = NULL, xlab = "seed yield",  cex.lab = cex_lab, cex.axis = cex_axis)

Assumptions Not Met

Heteroscedasticity [2]

m3 <- lme(yield ~ variety, random = ~ 1 |block, 
          na.action = na.exclude, data = var_ex2)
plot(m3, main = "Classic Funnel Shape", abline = 0)

Assumptions Not Met

Heteroscedasticity [2]

Solution: model variance as a function of the model fitted values using varPower(), an ‘nlme’ function

m3_new <- lme(yield ~ variety, random = ~ 1 |block, na.action = na.exclude, data = var_ex2,
              weights = varPower())
plot(m3_new)

Assumptions Not Met

Heteroscedasticity [2]

emmeans(m3, ~ variety, at = list(variety = c("variety_22", "variety_3", "variety_26"))) 
 variety    emmean   SE df lower.CL upper.CL
 variety_22   28.5 59.8  3   -161.9      219
 variety_3   168.1 59.8  3    -22.3      359
 variety_26  725.3 59.8  3    534.9      916

Degrees-of-freedom method: containment 
Confidence level used: 0.95 
emmeans(m3_new, ~ variety, at = list(variety = c("variety_22", "variety_3", "variety_26")))
 variety    emmean    SE df lower.CL upper.CL
 variety_22   28.5  6.43  3      8.0     48.9
 variety_3   168.1 27.70  3     79.9    256.3
 variety_26  725.3 92.30  3    431.5   1019.1

Degrees-of-freedom method: containment 
Confidence level used: 0.95

Assumptions Not Met

Heteroscedasticity [2]

anova(m3, type = "marginal") 
            numDF denDF F-value p-value
(Intercept)     1   111 33.3863  <.0001
variety        37   111  9.6935  <.0001
anova(m3_new, type = "marginal") 
            numDF denDF  F-value p-value
(Intercept)     1   111 47.47258  <.0001
variety        37   111 29.91474  <.0001


Note: LSD (least significant difference) calculations rely on a single variance term, so LSD cannot be used on these models (and is ill defined for mixed models)

Hypothesis Testing

  • Tests for unequal variance by groups

    • Bartlett’s test: bartlett.test()
    • Levene’s test: car::leveneTest(..., center = "mean")
    • Brown-Forsythe test: car::leveneTest(..., center = "median")

  • Test for unequal variance for regression: lmtest::bptest()

  • Log likelihood ratio tests between models 🤩

Hypothesis Testing


  • Barlett’s test for unequal variance (\(H_0\): group variances are equal)
bartlett.test(yield ~ site, data = var_ex1)

    Bartlett test of homogeneity of variances

data:  yield by site
Bartlett's K-squared = 567.26, df = 8, p-value < 2.2e-16

Hypothesis Testing


- Log likelihood ratio test (\(H_0\): the models are equivalent)
- “Better” model: lower AIC, lower BIC, higher log likelihood
- Test to determine if change in log likelihood is meaningful

anova(m2, m2_new)
       Model  df      AIC      BIC    logLik   Test  L.Ratio p-value
m2         1 345 14814.25 16511.55 -7062.128                        
m2_new     2 353 14489.23 16225.88 -6891.615 1 vs 2 341.0259  <.0001

SAS Options

For evaluating homoscedasticity


proc glm data=var_ex1;
  class site variety rep;
  model yield = variety|site rep(site);
  means site / hovtest=bartlett bf levene;
run;


Proc GLM Documentation

SAS Options

For modelling heteroscedasticity


proc mixed data=var_ex1;
  class site variety rep;
  model yield = site*variety + variety;
  random rep(site);
  repeated / group = site;
run;


Proc Mixed Documentation

Intractable iid Violations

Lognormal Models

Almost normal, just take the (natural) log

Earthquake acceleration

Lognormal Models

m4 <- gls(accel ~ Quake, data = Earthquake)
plot(m4, abline = 0)

Lognormal Models

m4_new <- gls(log(accel) ~ Quake, data = Earthquake)
plot(m4_new, abline = 0)

Lognormal Models

anova(m4) 
Denom. DF: 159 
            numDF   F-value p-value
(Intercept)     1 278.13410  <.0001
Quake          22   4.50921  <.0001
anova(m4_new) 
Denom. DF: 159 
            numDF   F-value p-value
(Intercept)     1 1232.4445  <.0001
Quake          22    6.8964  <.0001

Lognormal Models

emmeans(m4, "Quake", at = list(Quake = c("20", "10", "5")))
 Quake emmean     SE  df lower.CL upper.CL
 20    0.1630 0.0312 159   0.1014    0.225
 10    0.0300 0.1250 159  -0.2164    0.276
 5     0.0275 0.0376 159  -0.0468    0.102

Degrees-of-freedom method: df.error 
Confidence level used: 0.95 
emmeans(m4_new, "Quake", type = "response", at = list(Quake = c("20", "10", "5")))
 Quake response     SE  df lower.CL upper.CL
 20      0.1439 0.0335 159  0.09086   0.2281
 10      0.0300 0.0280 159  0.00476   0.1890
 5       0.0156 0.0044 159  0.00898   0.0272

Degrees-of-freedom method: df.error 
Confidence level used: 0.95 
Intervals are back-transformed from the log scale

Other Generalized Models

  • Poisson: the mean equals the variance
  • negative binomial, binomial (logistic), gamma, geometric, beta, …(more): the variance is a mathematical function of the mean
  • common deviations: count data (unless even the lowest counts are very high), percent data that might be subject to bounding of the values, binomial outcomes, multinomial outcomes, highly skewed data
  • A generalized linear model is needed; these use another distribution to describe the behavior of a dependent variable (not the residuals). Learning these could constitute a semester-long course.

Final Thoughts

  • In standard linear models, there is an ‘iid’ assumption (independent and identically distributed): we assume that model residuals shares a single common variance and are uncorrelated with each other

  • Violating these assumptions can impact the power of hypothesis testing and the width of confidence intervals

  • Failure to address these violations invalidates the analysis and inference from it

  • We can model variance by categorical variable or by a covariate using tools in the R package ‘nlme’

  • Some variables by their nature violate the iid assumption and required generalized linear models to analyze

Additional Resources