Ordinary Least Squares (OLS) Linear Regression

What OLS is and how it works

Ordinary Least Squares (OLS) linear regression finds the line — or hyperplane in the multivariable case — that minimizes the sum of squared vertical distances between each observed outcome yᵢ and its model-predicted value ŷᵢ = β₀ + β₁x₁ᵢ + … + βₚxₚᵢ. The name "least squares" refers precisely to this criterion: minimize Σ(yᵢ − ŷᵢ)². For a single predictor, the closed-form solution is β₁ = Σ(xᵢ − x̄)(yᵢ − ȳ) / Σ(xᵢ − x̄)², which equals the sample covariance of x and y divided by the sample variance of x.

The intercept is β₀ = ȳ − β₁x̄. In the multivariable model, each coefficient β_j represents the expected change in the mean of y per one-unit increase in x_j, holding all other predictors fixed. This "holding other predictors fixed" is precisely what "covariate adjustment" means in observational research — it estimates the association within strata defined by the other covariates.

Assumptions: coefficient validity vs inference validity

A critical distinction that analysts frequently blur: the conditions for OLS coefficients to be unbiased are different from the conditions needed for valid standard errors and CIs.

For coefficient validity (unbiasedness), OLS requires: (1) the conditional mean E(y|X) is correctly specified as linear in X; and (2) the errors are mean-zero and uncorrelated with the predictors (exogeneity, i.e., no omitted confounders correlated with X). Violations of linearity or exogeneity bias the coefficients themselves and no standard error correction can fix this.

For inference validity (correct SEs and CIs), OLS additionally requires: (3) homoscedasticity — constant error variance across observations; and (4) independent errors across observations. The normality of residuals is a fifth classical assumption that matters far less in practice: at the sample sizes typical of RWE studies (thousands to millions of patients), the Central Limit Theorem ensures the sampling distribution of OLS coefficients is approximately normal regardless of the residual distribution. Non-normality of residuals is not a meaningful concern in large claims or EHR analyses (see Lumley et al. 2002).

Heteroscedasticity — non-constant error variance — is endemic in health data: costs are right-skewed, sicker patients have more variable utilization, and length of stay fans out with disease severity. Heteroscedasticity biases conventional (homoscedastic) SEs but does not bias the coefficients. The practical solution is to use heteroscedasticity-robust standard errors — specifically the HC3 variant due to MacKinnon and White — as the default in all OLS analyses of health data (Long and Ervin 2000). HC3 inflates the influence of high-leverage observations, providing better finite-sample coverage than HC0 or HC1. Conventional SEs should be reserved for sensitivity checks and clearly labeled as such.

Interpreting the output

Consider a multivariable OLS model of hospital length of stay (LOS, days) regressed on a binary treatment indicator and five covariates: age, sex, Charlson comorbidity score, index year, and payer type. The adjusted coefficient on treatment is β = −1.8 days (95% CI −3.1 to −0.5, p = 0.007), and the model R² = 0.32.

Formal interpretation. β = −1.8 is the model-estimated difference in mean LOS comparing treated to untreated patients at fixed values of age, sex, Charlson score, index year, and payer type — a confounding-adjusted association within those covariate strata. The 95% CI (−3.1 to −0.5) is a repeated-sampling interval: if the study were repeated many times under identical conditions, approximately 95% of such intervals would contain the true adjusted mean difference — it does NOT mean "there is a 95% probability the true effect is between −3.1 and −0.5" (that is a Bayesian credible interval). R² = 0.32 means 32% of the total variance in LOS is explained by the model collectively; it does not measure whether coefficients are causal, whether the model predicts well in a new population, or whether the unexplained variance (68%) is clinically negligible.

Practical interpretation. On average, treated patients in this cohort stayed about 1.8 fewer days than untreated patients with similar demographics, comorbidity burden, and payer mix — an adjusted association compatible with a true difference plausibly ranging from half a day to just over three days based on the confidence interval. If the included covariates adequately capture confounding, this is a reasonable proxy for the causal effect. If important confounders are unmeasured (e.g., functional status, disease severity not captured by Charlson), the estimate remains a biased adjusted association. R² = 0.32 tells us the model explains a moderate share of LOS variability — appropriate context for setting expectations about residual individual variation, but not a measure of causal completeness.

Adjustment, dummy coding, and the Table-2 fallacy

In a multivariable OLS model, each coefficient is estimated conditional on all other model terms. Continuous covariates are "held fixed" in the mathematical sense of partial derivatives. Binary exposures coded 0/1 (dummy coding) produce a coefficient equal to the adjusted mean difference. Categorical variables with k categories require k−1 dummy indicators with one reference category. Interactions between predictors model departure from additivity on the mean-difference scale and must be pre-specified.

The Table-2 fallacy is a high-value RWE warning: the coefficients on confounders in a multivariable regression model are NOT causal effects of those confounders and should not be reported as such. In an OLS model of cost adjusted for age, sex, comorbidity, and treatment, the coefficient on age is the estimated association of age with cost after conditioning on sex, comorbidity, and treatment — a different causal quantity from the total effect of age, often biased if any covariates are mediators of or colliders with the age-cost pathway. Only the coefficient on the pre-specified primary exposure has the intended adjustment interpretation. Report all other coefficients as associations with appropriate caveats, or suppress them from primary results tables.

Clustering: patients within providers and plans

When patients are nested within clusters — hospitals, medical practices, insurance plans, geographic regions — their outcomes are correlated within clusters even after controlling for patient-level covariates. Ignoring this correlation understates standard errors and inflates false-positive rates. The solution is cluster-robust standard errors (also called clustered SEs or the Huber-White sandwich with cluster correction), which require only that the clusters themselves are a random sample from a large population of clusters. A common rule of thumb is at least 30 clusters for reliable cluster-robust inference; with fewer clusters, bootstrap methods or the CR2 small-sample correction should be considered. Cluster-robust SEs are available across all major packages: statsmodels `cov_type="cluster"`, R `coeftest(..., vcov=vcovCL(..., cluster=...))`, SAS PROC SURVEYREG. For claims data with patients nested in health plans, or EHR data with patients nested in hospitals, cluster- robust SEs should be the default when OLS is used for comparative analyses.

OLS on cost and utilization outcomes: when it fails and when it can be defended

OLS applied directly to raw healthcare costs faces three structural problems:

1. right skew — the distribution has a long right tail driven by high-cost outliers, so the residuals are far from symmetric;
2. heteroscedasticity — variance of costs rises with the mean (sicker patients have higher mean and higher variance); and
3. mass at zero — many patients have zero cost in a window, creating a spike that no continuous linear model accommodates well. These problems motivate the gamma GLM with a log link (or a two-part model when zeros are structurally informative) as the modern standard for primary inference on mean cost differences — these models respect positivity of costs, accommodate skew through the variance function, and estimate mean ratios directly interpretable for payers.

However, OLS is not always wrong for cost data. In very large samples, the CLT ensures that the OLS estimator of the mean difference is consistent and approximately normally distributed even with skewed residuals — the coefficient is valid for inference on the mean difference (the policy-relevant quantity for budget-impact models). When the target estimand is a difference in means rather than a ratio, when zero-inflation is modest, and when robust SEs are used, OLS on raw or log-transformed costs can be a defensible sensitivity analysis or quick descriptive estimate. The principled choice remains a GLM for confirmatory analysis; OLS serves as a transparency check.

Pros, cons, and trade-offs

Pros: OLS coefficients are directly interpretable as differences in the conditional mean of y — the most intuitive scale for clinical and policy audiences. The method is computationally trivial, extends naturally to multivariable adjustment, dummy coding, and interactions, and has a closed-form solution. With robust SEs, OLS yields valid inference on means even under heteroscedasticity. In large RWE datasets, non-normality of residuals is irrelevant (CLT), so the classical criticism of OLS as "requiring normal data" is largely misplaced for large health studies. Predictions are on the original outcome scale, and coefficients have direct units (days, dollars, points on a score).

Cons: OLS assumes a linear conditional mean — if the true relationship is substantially nonlinear, coefficients are biased. R² is easily misread as a measure of causal explanatory power or predictive accuracy; it is neither. For skewed outcomes, OLS can produce negative predictions, suffers from heteroscedastic residuals that inflate conventional SEs, and has lower efficiency than a correctly specified GLM. Confounding bias is not addressed by OLS itself; it requires a well-designed covariate set or a causal identification strategy.

Trade-offs vs GLMs: The gamma GLM with log link is preferred for cost and utilization outcomes because it respects positivity and accommodates heteroscedasticity through the variance function. OLS is preferred when the effect scale is an additive mean difference, when interactions on the additive scale are the estimand, or when direct coefficient interpretability in original units is required.

When NOT to use

Do not use OLS when:

• the outcome is binary — use logistic regression or modified Poisson for risk differences
• the outcome is a count — use Poisson or negative binomial regression
• the outcome is time-to-event with censoring — use Cox proportional-hazards or parametric survival models
• the outcome is bounded or ordinal — use ordinal regression
• the sample is small (n < 30 per arm) with clearly non-normal residuals and no CLT protection
• or raw cost/utilization is the outcome and a mean-ratio estimand is needed — use a gamma GLM or two-part model for confirmatory analysis.

Do not use OLS and report causal effects without a credible identification design — an adjusted regression coefficient is a confounding- adjusted association, not a causal effect, unless the study design (RCT, valid instrumental variable, sharp regression discontinuity) provides identification.

The Table-2 fallacy warning applies: confounders' coefficients are not their causal effects and should not be narrated as such.

For heavy-tailed distributions at small n, robust SEs reduce but do not eliminate the risk of poorly calibrated inference.

Decision diagram

flowchart TD
  Q[Is the outcome<br/>continuous?] --> No[No → use logistic, Poisson,<br/>Cox, or ordinal regression]
  Q --> Yes[Yes]
  Yes --> Cost[Is the outcome<br/>raw cost / LOS?]
  Cost --> CostYes["Consider gamma GLM<br/>as primary analysis<br/>(see gamma-distribution)"]
  Cost --> CostNo[No — lab values, utility<br/>scores, log-cost, etc.]
  CostNo --> ClusterQ[Are patients nested<br/>in providers / plans?]
  CostYes --> CostOLS["OLS as robustness check<br/>(large n only, HC3 SEs,<br/>label as approximation)"]
  ClusterQ --> Cluster["Yes → cluster-robust SEs<br/>(≥30 clusters)"]
  ClusterQ --> NoCluster["No → HC3 robust SEs<br/>(default for health data)"]
  Cluster --> FitOLS[Fit multivariable OLS<br/>with cluster-robust SEs]
  NoCluster --> FitOLS
  FitOLS --> Diagnose["Diagnose: residual plots,<br/>leverage, influential obs<br/>(see regression-diagnostics)"]
  Diagnose --> Report["Report: coefficient + 95% CI<br/>in original units; label as<br/>adjusted association unless<br/>causal design is established"]

Worked example

Scenario

A health outcomes researcher fits a simple OLS regression to explore the relationship between number of prior hospitalizations in the 12 months before an index date (x, ranging from 1 to 5) and total annual healthcare cost in the following 12 months (y, in thousands of dollars) across 5 patients from a registry. The analyst works through the exact least-squares arithmetic by hand to understand where the slope and intercept come from and how R-squared is computed from the residuals.

Dataset

Five-patient registry excerpt. x = prior hospitalizations (count); y = annual cost ($K). Patient 3 has unusually high cost (12) relative to the linear trend, which will produce a larger residual and drive the R-squared below 1.

Steps

Result

beta_0 = 2 (intercept, $2K baseline cost), beta_1 = 2 (each additional hospitalization adds $2K predicted cost). Fitted values: 4, 6, 8, 10, 12. Residuals: -1, -1, 4, -1, -1. Residual sum check: (-1)+(-1)+4+(-1)+(-1) = 0. SST = 25+9+16+1+9 = 60, SSE = 1+1+16+1+1 = 20, SSR = 60-20 = 40. R-squared = 40/60 = 0.667.

Trade-offs

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm

# ── Toy dataset (matches beginner-layer worked example) ──────────────────────────────
data_toy = pd.DataFrame({
    "prior_hosp": [1, 2, 3, 4, 5],
    "annual_cost_k": [3, 5, 12, 9, 11],
})
X_toy = sm.add_constant(data_toy[["prior_hosp"]])   # adds intercept column
y_toy = data_toy["annual_cost_k"]

# ── Fit OLS ──────────────────────────────────────────────────────────────────────────
model_toy = sm.OLS(y_toy, X_toy).fit()
print("=== Toy dataset (n=5): conventional SEs ===")
print(model_toy.summary2())
# Intercept = 2.0, slope (prior_hosp) = 2.0, R-squared = 0.667

# ── HC3 robust standard errors (default for health data) ─────────────────────────────
model_hc3 = model_toy.get_robustcov_results(cov_type="HC3")
print("\n=== Toy dataset: HC3 robust SEs ===")
print(model_hc3.summary2())
# In a 5-obs example HC3 CIs will be very wide; robust SEs matter in larger samples.

# ── Larger simulated RWE dataset: OLS with HC3 SEs ──────────────────────────────────
rng = np.random.default_rng(42)
n = 2000
df = pd.DataFrame({
    "treated":      rng.binomial(1, 0.5, n),
    "age":          rng.normal(65, 12, n),
    "charlson":     rng.poisson(2, n),
    "female":       rng.binomial(1, 0.52, n),
})
# Simulate LOS outcome (days): true adjusted treatment effect = -2.0 days
df["los"] = (
    8.0
    - 2.0 * df["treated"]
    + 0.08 * (df["age"] - 65)
    + 0.4 * df["charlson"]
    - 0.5 * df["female"]
    + rng.normal(0, 3.5, n)   # heteroscedastic-ish residuals
)

covariates = ["treated", "age", "charlson", "female"]
X = sm.add_constant(df[covariates])
y = df["los"]

# Fit with HC3 robust SEs — single call via cov_type argument
ols_hc3 = sm.OLS(y, X).fit(cov_type="HC3")

print("\n=== Multivariable OLS, HC3 robust SEs (n=2000) ===")
print(ols_hc3.summary2())
# Focus on the 'treated' coefficient: should be near -2.0
coef_trt = ols_hc3.params["treated"]
ci_trt   = ols_hc3.conf_int().loc["treated"]
print(f"\nAdjusted treatment effect: {coef_trt:.2f} days")
print(f"95% CI (HC3): [{ci_trt[0]:.2f}, {ci_trt[1]:.2f}]")
print(f"Model R-squared: {ols_hc3.rsquared:.3f}")

# ── Cluster-robust SEs (patients nested in 50 hypothetical plans) ───────────────────
df["plan_id"] = rng.integers(1, 51, n)   # 50 insurance plans
ols_cluster = sm.OLS(y, X).fit(
    cov_type="cluster",
    cov_kwds={"groups": df["plan_id"]}
)
ci_trt_cl = ols_cluster.conf_int().loc["treated"]
print(f"\nCluster-robust 95% CI (clustered by plan): [{ci_trt_cl[0]:.2f}, {ci_trt_cl[1]:.2f}]")
print("Note: cluster-robust SEs are wider than HC3 when within-plan correlation is positive.")

library(sandwich)
library(lmtest)

# ── Toy dataset (matches beginner-layer worked example) ──────────────────────────────
toy <- data.frame(
  prior_hosp    = 1:5,
  annual_cost_k = c(3, 5, 12, 9, 11)
)
m_toy <- lm(annual_cost_k ~ prior_hosp, data = toy)
cat("=== Toy dataset (n=5): conventional SEs ===\n")
print(summary(m_toy))
# Intercept = 2, slope (prior_hosp) = 2, R-squared = 0.6667

# ── HC3 robust SEs via sandwich + lmtest ─────────────────────────────────────────────
cat("\n=== Toy dataset: HC3 robust SEs ===\n")
print(coeftest(m_toy, vcov = vcovHC(m_toy, type = "HC3")))

# ── Larger simulated RWE dataset ─────────────────────────────────────────────────────
set.seed(42)
n <- 2000
df <- data.frame(
  treated  = rbinom(n, 1, 0.5),
  age      = rnorm(n, 65, 12),
  charlson = rpois(n, 2),
  female   = rbinom(n, 1, 0.52),
  plan_id  = sample(1:50, n, replace = TRUE)
)
df$los <- (
  8.0
  - 2.0 * df$treated
  + 0.08 * (df$age - 65)
  + 0.4 * df$charlson
  - 0.5 * df$female
  + rnorm(n, 0, 3.5)
)

m <- lm(los ~ treated + age + charlson + female, data = df)
cat("\n=== Multivariable OLS, conventional SEs (n=2000) ===\n")
print(summary(m))

# HC3 robust SEs — recommended default for health data
cat("\n=== HC3 robust SEs ===\n")
print(coeftest(m, vcov = vcovHC(m, type = "HC3")))

# Confidence interval with HC3 SEs
ci_hc3 <- coefci(m, vcov = vcovHC(m, type = "HC3"))
cat(sprintf(
  "\nAdjusted treatment effect: %.2f days, 95%% CI (HC3): [%.2f, %.2f]\n",
  coef(m)["treated"], ci_hc3["treated", 1], ci_hc3["treated", 2]
))

# ── Cluster-robust SEs (patients nested in 50 insurance plans) ───────────────────────
cat("\n=== Cluster-robust SEs (clustered by plan_id) ===\n")
print(coeftest(m, vcov = vcovCL(m, cluster = ~plan_id, data = df)))
# Note: requires sandwich >= 2.5-1 for the cluster formula interface
cat("Rule of thumb: reliable with >= 30 clusters; 50 plans is adequate here.\n")

/* ── Create motivating dataset (toy: n=5) ─────────────────────────────────────── */
data work.toy;
  input prior_hosp annual_cost_k;
  datalines;
1  3
2  5
3 12
4  9
5 11
;
run;

/* ── PROC REG: OLS with ACOV (heteroscedasticity-consistent SEs, HC0 variant) ─── */
proc reg data=work.toy;
  model annual_cost_k = prior_hosp / acov;
  /* ACOV: "Asymptotic Covariance" — produces White (HC0) robust SEs.
     Use the ACOV rows in the output for robust inference.
     Intercept = 2, slope (prior_hosp) = 2, R-squared = 0.6667              */
run;

/* ── Create larger simulated RWE dataset ─────────────────────────────────────── */
data work.rwe;
  call streaminit(42);
  do i = 1 to 2000;
    treated  = rand("bernoulli", 0.5);
    age      = rand("normal", 65, 12);
    charlson = rand("poisson", 2);
    female   = rand("bernoulli", 0.52);
    plan_id  = ceil(rand("uniform") * 50);   /* 50 insurance plans */
    los = 8.0
          - 2.0 * treated
          + 0.08 * (age - 65)
          + 0.4 * charlson
          - 0.5 * female
          + rand("normal", 0, 3.5);
    output;
  end;
run;

/* ── PROC REG: multivariable OLS with HC0 robust SEs ──────────────────────────── */
proc reg data=work.rwe;
  model los = treated age charlson female / acov;
  /* Focus on ACOV row for the 'treated' coefficient.
     ACOV in PROC REG implements HC0 (White 1980). For HC3, use PROC IML
     or the %HCCME macro available from SAS support.                          */
run;

/* ── PROC GLM: equivalent OLS fit (useful for CLASS variables and LS means) ──── */
proc glm data=work.rwe;
  model los = treated age charlson female;
  /* PROC GLM does not have a built-in ACOV option; combine with PROC IML
     or PROC SURVEYREG (below) for robust SEs.                                */
run;

/* ── PROC SURVEYREG: cluster-robust SEs (the recommended approach in SAS) ────── */
proc surveyreg data=work.rwe;
  cluster plan_id;              /* cluster-robust SEs at the plan level         */
  model los = treated age charlson female;
  /* PROC SURVEYREG implements the Huber-White sandwich estimator with
     cluster correction. Use this procedure whenever patients are nested
     within plans, hospitals, or practices.
     Rule of thumb: >= 30 clusters for reliable inference; 50 plans here.     */
run;

Citations

1. [1]Schneider A, Hommel G, Blettner M. Linear regression analysis. Deutsches Arzteblatt International. 2010;107(44):776-782.
2. [2]Lumley T, Diehr P, Emerson S, Chen L. The importance of the normality assumption in large public health data sets. Annual Review of Public Health. 2002;23:151-169.

1. [3]Long JS, Ervin LH. Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician. 2000;54(3):217-224.

1. [4]White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838.
2. [5]Manning WG, Mullahy J. Estimating log models: to transform or not to transform? Journal of Health Economics. 2001;20(4):461-494.