Answer

Do robust standard errors fix floor and ceiling effects?

Inspired by a question on Cross Validated · November 8, 2026

The short answer

Only partly. Robust (sandwich) standard errors give honest intervals and p values when score variance differs between groups, which a ceiling or floor usually causes. They leave the estimates untouched, so you still get the difference in observed means, which a scale limit can squeeze well below the difference in the underlying trait. If the observed mean is your target, use HC3. If not, change the model.

What robust standard errors change, and what they leave alone

A least-squares fit gives you coefficients and standard errors. The classic standard error assumes the residual variance is the same for everyone. A heteroscedasticity-consistent standard error (also called robust or sandwich) drops that assumption and uses each observation's own squared residual instead. The coefficients are calculated exactly as before. Only the uncertainty around them changes.

So this is a repair for one specific problem. It does not move an estimate, it does not make a straight line fit any better, and it does not handle clustered or repeated observations, which need cluster-robust methods or a mixed model. In a model with only a grouping factor the link to a familiar test is exact: the version known as HC2 reproduces the standard error of the Welch t test. A robust standard error is essentially Welch's idea carried over to models with covariates.

The two things a scale limit does

It distorts the variances. A group close to the limit gets squashed, so its scores vary less than those of a group far from the limit, even when the underlying spread is identical. That is heteroscedasticity, and robust standard errors deal with it.
It changes what the difference means. A person at 100 might be just at that level or far beyond it, yet both count as 100. When one group sits near the limit, the gap in observed means is usually smaller than the gap in the underlying trait. Robust standard errors do nothing about this.

Whether the second point matters depends on your question. If the top of the scale really is the top, the mean of the observed score may be precisely what you want to compare. If the score is a stand-in for a trait that continues past the end of the scale, it is not.

See it in R and Python

Here two groups (25 and 1,500 people) have the same underlying spread, with latent means 35 points apart, and scores are cut off at 0 and 100. The script writes out the classic, HC1 and HC3 standard errors in full, fits one data set, then repeats the fit 3,000 times to see how often each 95% interval covers the true difference in observed means. Figures quoted come from one seeded run of the R code. Python's random numbers differ, so its simulated figures differ slightly but show the same pattern (the exact 27.84 target is identical).

set.seed(677562)
n <- c(small = 25, large = 1500)
group <- factor(rep(names(n), n), levels = c("small", "large"))

# Both groups have the SAME latent spread; only the 0-100 scale differs
draw <- function() {
  latent <- rnorm(sum(n), ifelse(group == "small", 60, 95), 25)
  pmin(pmax(latent, 0), 100)
}

# Classic, HC1 and HC3 standard errors for one coefficient of any lm fit
se_types <- function(fit, term = 2) {
  X <- model.matrix(fit); e <- resid(fit); h <- hatvalues(fit)
  N <- nrow(X); p <- ncol(X); bread <- solve(crossprod(X))
  sandwich <- function(w) (bread %*% crossprod(X * w) %*% bread)[term, term]
  c(classic = sqrt(sum(e^2) / (N - p) * bread[term, term]),
    HC1 = sqrt(N / (N - p) * sandwich(e)),
    HC3 = sqrt(sandwich(e / (1 - h))))
}

# Exact mean of a normal score clipped to 0-100
clipped_mean <- function(mu, s = 25) {
  a <- (0 - mu) / s; b <- (100 - mu) / s
  mu * (pnorm(b) - pnorm(a)) + s * (dnorm(a) - dnorm(b)) + 100 * (1 - pnorm(b))
}
truth <- clipped_mean(95) - clipped_mean(60)
round(c(latent_difference = 35, observed_difference = truth), 2)

# One data set
y <- draw()
fit <- lm(y ~ group)
round(rbind(sd = tapply(y, group, sd), at_100 = tapply(y == 100, group, mean)), 2)
se <- se_types(fit); b <- unname(coef(fit)[2]); df <- fit$df.residual
round(c(estimate = b, se), 2)
round(c(t_HC3 = b / se[["HC3"]], df = df,
        lower = b - qt(.975, df) * se[["HC3"]],
        upper = b + qt(.975, df) * se[["HC3"]]), 2)

# Repeat: how often does a 95% interval cover the observed-score difference?
res <- t(replicate(3000, {
  y <- draw(); fit <- lm(y ~ group); b <- unname(coef(fit)[2])
  c(estimate = b, abs(b - truth) <= qt(.975, fit$df.residual) * se_types(fit))
}))
round(colMeans(res), 3)

Python

import numpy as np
from scipy import stats

rng = np.random.default_rng(677562)
n_small, n_large = 25, 1500
large = np.r_[np.zeros(n_small), np.ones(n_large)]
X = np.column_stack([np.ones(n_small + n_large), large])   # reference: small group

# Both groups have the SAME latent spread; only the 0-100 scale differs
def draw():
    latent = rng.normal(np.where(large == 0, 60, 95), 25)
    return np.clip(latent, 0, 100)

# Classic, HC1 and HC3 standard errors for one coefficient of an OLS fit
def fit_ols(X, y, term=1):
    N, p = X.shape
    bread = np.linalg.inv(X.T @ X)
    b = bread @ X.T @ y
    e = y - X @ b
    h = np.sum((X @ bread) * X, axis=1)
    def sandwich(w):
        Xw = X * w[:, None]
        return (bread @ (Xw.T @ Xw) @ bread)[term, term]
    se = {"classic": np.sqrt(e @ e / (N - p) * bread[term, term]),
          "HC1": np.sqrt(N / (N - p) * sandwich(e)),
          "HC3": np.sqrt(sandwich(e / (1 - h)))}
    return b[term], se, N - p

# Exact mean of a normal score clipped to 0-100
def clipped_mean(mu, s=25):
    a, b = (0 - mu) / s, (100 - mu) / s
    return (mu * (stats.norm.cdf(b) - stats.norm.cdf(a))
            + s * (stats.norm.pdf(a) - stats.norm.pdf(b))
            + 100 * (1 - stats.norm.cdf(b)))

truth = clipped_mean(95) - clipped_mean(60)
print("latent difference 35, observed difference", round(truth, 2))

# One data set
y = draw()
for label, g in [("small", 0), ("large", 1)]:
    print(label, "sd", round(y[large == g].std(ddof=1), 2),
          "at 100", round(np.mean(y[large == g] == 100), 2))
b, se, df = fit_ols(X, y)
print("estimate", round(b, 2), {k: round(float(v), 2) for k, v in se.items()})
margin = stats.t.ppf(.975, df) * se["HC3"]
print("t_HC3", round(b / se["HC3"], 2), "df", df,
      "95% CI", round(b - margin, 2), round(b + margin, 2))

# Repeat: how often does a 95% interval cover the observed-score difference?
reps = 3000
estimates = np.empty(reps)
hits = {"classic": 0, "HC1": 0, "HC3": 0}
for i in range(reps):
    b, se, df = fit_ols(X, draw())
    estimates[i] = b
    for k in hits:
        hits[k] += abs(b - truth) <= stats.t.ppf(.975, df) * se[k]
print("mean estimate", round(estimates.mean(), 3),
      {k: round(v / reps, 3) for k, v in hits.items()})

Output from one seeded run of the R code (same numbers on a second run):
  latent_difference observed_difference 
              35.00               27.84 
       small large
sd     22.27 16.45
at_100  0.00  0.42
estimate  classic      HC1      HC3 
   24.96     3.34     4.39     4.57 
  t_HC3      df   lower   upper 
   5.47 1523.00   16.00   33.92 
estimate  classic      HC1      HC3 
  27.885    0.824    0.934    0.943

The cap turns a 35-point latent gap into an exact observed gap of 27.84. In the one data set, 42% of the large group sits at 100 and nobody in the small group does, so a shared latent SD of 25 becomes observed SDs of 16.45 and 22.27. The ceiling created the unequal variances. The estimate is 24.96 whichever standard error you choose, and the standard error is 3.34 (classic), 4.39 (HC1) or 4.57 (HC3).

Over 3,000 data sets, the average estimate is 27.885, right on the observed gap and far from 35. Classic 95% intervals cover the observed gap only 82.4% of the time, because the small group is the more variable one and the pooled variance is dominated by the large group. HC1 reaches 93.4% and HC3 94.3%. The robust intervals are now honest, but they are honest about the compressed difference.

Which robust version to use

HC0 is White's (1980) original. It runs too small in modest samples.
HC1 multiplies by N/(N − p). It is what Stata's robust option gives you.
HC3 inflates each squared residual by 1/(1 − h)², where h is the observation's leverage. Long and Ervin (2000) recommended it for samples under 250. In R, use sandwich::vcovHC(fit, type = "HC3") with lmtest::coeftest().

With unequal groups, the total N can mislead you. In a model with only group indicators, each member of a group of size n has leverage 1/n, so what counts is the size of your smallest group. HC3 did well with 25 people in the example; with only a handful in a group, any of these large-sample corrections can be unreliable.

When to change the model instead

Ordinal regression (MASS::polr, ordinal::clm) uses only the order of the scores, so a pile at the top is simply the highest category. You get an odds ratio rather than a difference in points.
Tobit regression treats 100 as "at least 100" and estimates the latent gap. It relies on a normal latent variable with constant variance, which is hard to defend if the ceiling comes from how the items were built, such as a few yes/no questions.
Fractional logit divides the score by its maximum and fits a quasi-binomial model with robust standard errors. Exact 0s and 1s are allowed, and fitted means stay inside the scale.
A two-part analysis models reaching the limit with logistic regression and then models scores below it, which suits data where being at the limit is meaningfully different.
The bootstrap replaces the standard error, not the target. It still estimates the difference in observed means. Resample within groups so that group sizes stay fixed.

If you run the same analysis on many correlated outcomes, one method across all of them makes results comparable, and robust least squares is a defensible choice when the observed mean is meaningful. Add an ordinal model as a sensitivity check for the outcomes with the heaviest piles, and decide in advance how to handle multiple comparisons. The DASS tools page has calculators for the simpler group comparisons.

How to report robust standard errors in APA style (7th edition)

Name the version, and give the estimate, its robust standard error and the interval. Using the single data set from the example:

"Mean scores were higher in the large group than in the small group (b = 24.96, HC3 SE = 4.57), t(1523) = 5.47, p < .001, 95% CI [16.00, 33.92]."

In the method section, say why: "Because ceiling effects produced unequal residual variances, we used heteroscedasticity-consistent standard errors (HC3; Long & Ervin, 2000)." If a limit may have shrunk the difference, say that in the discussion too, since the robust interval does not account for it.

Related tools and guides

Working with your own data?

General answers only go so far. Send us your situation and we'll reply by email within two business days.

Ask your question

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.