Student-t / robust Linear Model — R companion

Companion to student_t_lm_python.ipynb (from-scratch Geweke Gibbs) and student_t_lm_pymc.ipynb (NUTS). The R engines for the same model, on the same data:

engine	what it is
hett::tlm	Student-t linear model by maximum likelihood — estimates β, the scale, and ν (the dof). The frequentist-ML counterpart of the scale-mixture Bayesian model.
MASS::rlm	classic robust regression (M-estimation, Huber ψ) — down-weights outliers without assuming a specific error law.
lm	OLS — the non-robust baseline.

Same Section 1 synthetic data ($n=120$, $t(4)$ errors, $\beta=[2,0.8,-0.5]$) used by the Python/PyMC notebooks.

suppressMessages({library(hett); library(MASS); library(ggplot2)})
d <- read.csv('studentt_synth.csv')
ols <- lm(y ~ x1 + x2, d)
rob <- rlm(y ~ x1 + x2, d)
tl  <- suppressWarnings(tlm(y ~ x1 + x2, data = d, estDof = TRUE))
cmp <- rbind(OLS = coef(ols), 'rlm (robust M)' = coef(rob), 'tlm (Student-t ML)' = coef(tl$loc.fit))
colnames(cmp) <- c('intercept','x1','x2')
cat('coefficient estimates (true: 2.0, 0.8, -0.5):\n'); print(round(cmp, 3))
cat(sprintf('\ntlm estimated dof  nu = %.2f   (true nu = 4; Python Gibbs mode 4, PyMC median ~6.6)\n', tl$dof))

coefficient estimates (true: 2.0, 0.8, -0.5):

                   intercept    x1     x2
OLS                    2.061 0.745 -0.675
rlm (robust M)         2.044 0.787 -0.673
tlm (Student-t ML)     2.019 0.805 -0.652

tlm estimated dof  nu = 3.93   (true nu = 4; Python Gibbs mode 4, PyMC median ~6.6)

# Robustness on a 1-predictor contaminated demo (same data as the Python student_t_robust figure)
options(repr.plot.width = 9, repr.plot.height = 5.5)
r1 <- read.csv('studentt_robust1d.csv')
o1 <- lm(y ~ x, r1); rb1 <- rlm(y ~ x, r1); t1 <- suppressWarnings(tlm(y ~ x, data = r1, estDof = TRUE))
co <- coef(o1); cr <- coef(rb1); ct <- coef(t1$loc.fit); r1$absres <- abs(resid(o1))
ggplot(r1, aes(x, y)) +
  geom_point(aes(colour = absres), size = 2.6) +
  scale_colour_viridis_c(direction = -1, name = '|OLS resid|\n(outliers dark)') +
  geom_abline(aes(intercept = co[1], slope = co[2], linetype = 'OLS'),  colour = 'steelblue', linewidth = 1.2) +
  geom_abline(aes(intercept = cr[1], slope = cr[2], linetype = 'rlm'),  colour = 'darkgreen', linewidth = 1.2) +
  geom_abline(aes(intercept = ct[1], slope = ct[2], linetype = 'tlm'),  colour = 'firebrick', linewidth = 1.2) +
  geom_abline(aes(intercept = 1,     slope = 2,     linetype = 'truth'), colour = 'black') +
  scale_linetype_manual(name = 'fit', values = c(OLS='solid', rlm='solid', tlm='solid', truth='dashed')) +
  labs(title = sprintf('OLS slope %.2f (dragged) vs rlm %.2f vs tlm %.2f vs truth 2.00', co[2], cr[2], ct[2]),
       x = 'x', y = 'y') +
  theme_minimal(base_size = 12)

The one-predictor robustness demo on the same $t(4)$-error data (points shaded by |OLS residual|, outliers dark): OLS (blue) is dragged off the true slope of 2.00 by the outliers, while rlm (Huber M-estimation, green) and tlm (Student-t ML, red) both track the truth. It is the frequentist echo of the Bayesian robustness figure — and tlm additionally recovers ν ≈ 3.9, matching the true 4 and the Gibbs / PyMC posteriors.

# Nelson-Plosser Real GNP: Student-t ML trend fit (compare nu to the Bayesian notebooks)
npd <- read.csv('nelplo.csv'); npd <- npd[npd$year >= 1909 & npd$year <= 1970, ]
g <- data.frame(t = as.numeric(scale(npd$year)), y = log(npd$gnp.real)); g <- g[!is.na(g$y), ]
tg <- suppressWarnings(tlm(y ~ t, data = g, estDof = TRUE))
cat(sprintf('Real GNP (log, linear trend):  tlm nu = %.1f   trend slope = %.3f\n', tg$dof, coef(tg$loc.fit)[2]))
cat('  (heavy tails -> low nu, consistent with the Gibbs/PyMC posteriors for Real GNP)\n')

Real GNP (log, linear trend):  tlm nu = 1.6   trend slope = 0.097

  (heavy tails -> low nu, consistent with the Gibbs/PyMC posteriors for Real GNP)

Results

• tlm (Student-t ML) recovers ν ≈ 3.9 on the synthetic data — matching the true ν=4 and the Bayesian posteriors (Gibbs mode 4, PyMC median ~6.6). It is the maximum-likelihood twin of the scale-mixture Gibbs: same model, a point estimate of ν instead of a posterior.
• β estimates agree across tlm, rlm, and the Bayesian fits; OLS is pulled toward the outliers (visible in the 1-D robustness plot — the blue OLS line is dragged off the truth while the rlm/tlm lines track it).
• rlm gives similar robustness via M-estimation but without an explicit error model or a ν; tlm is preferable when you want the heavy-tailed generative model and its dof.
• Real GNP: tlm estimates a low ν (heavy tails), consistent with the Bayesian notebooks — three engines (Gibbs, NUTS, ML) agree the macro series are heavy-tailed.
• For a fully Bayesian R version (continuous ν with a prior, like the PyMC notebook), brms::brm(y ~ x, family = student()) via Stan is the route — omitted here to avoid the Stan toolchain install.