Ordinary Least Squares (OLS) Linear Regression

In plain language

OLS linear regression draws the best-fitting straight line through a set of data points by choosing a slope and intercept that minimize the sum of squared gaps between the actual outcome values and the line's predictions. Each coefficient tells you how much the average outcome changes for a one-unit increase in that predictor, holding all other predictors constant — so it is a natural tool for comparing two groups while accounting for differences in age, sickness, or other characteristics. One honest caveat: OLS assumes the relationship is a straight line, and for outcomes like raw healthcare costs — which are right-skewed and always positive — a different type of model often fits better.

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 to use

Use OLS when: the outcome is continuous and the conditional mean is the target quantity; the relationship between predictors and outcome is plausibly linear over the observed covariate range; the sample is large enough for the CLT to protect inference (roughly n ≥ 100 per arm as a heuristic); and heteroscedasticity-robust or cluster-robust SEs are applied. OLS is appropriate for continuous endpoints in both RCTs and observational studies where confounders are measured and included: blood pressure, HbA1c, health utility scores, length of stay, patient-reported outcome scales. It is also the estimator inside difference-in-differences and event-study designs where linearity on the untreated outcome trend is the identifying assumption.

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.

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

• Compute the sample means. x-bar = (1+2+3+4+5)/5 = 15/5 = 3. y-bar = (3+5+12+9+11)/5 = 40/5 = 8.

• Compute the numerator of the OLS slope: the five cross-products (x_i - x-bar)*(y_i - y-bar) give values 10, 3, 0, 1, and 6 respectively (deviations for patients 1-5). Their sum: 10+3+0+1+6 = 20.

• Compute the denominator: the five squared deviations (x_i - x-bar)^2 are 4, 1, 0, 1, 4. Their sum: 4+1+0+1+4 = 10.

• OLS slope: beta_1 = 20/10 = 2. Interpretation: each additional prior hospitalization is associated with $2K higher predicted annual cost on average.

• OLS intercept: beta_0 = 8 - 2*3 = 8-6 = 2. The model predicts $2K cost for a patient with zero prior hospitalizations (an extrapolation beyond the data range here, so interpret cautiously).

• Compute fitted values (y-hat = 2x + 2): for x=1 the fitted value is 21+2 = 4; for x=2 it is 22+2 = 6; for x=3 it is 23+2 = 8; for x=4 it is 24+2 = 10; for x=5 it is 25+2 = 12.

• Compute residuals (y - y-hat): patient 1 gives 3-4 = -1; patient 2 gives 5-6 = -1; patient 3 gives 12-8 = 4; patient 4 gives 9-10 = -1; patient 5 gives 11-12 = -1. Residual sum check: (-1)+(-1)+4+(-1)+(-1) = 0. OLS residuals always sum to zero when the model includes an intercept.

• Compute R-squared. SST (total sum of squares) = 25+9+16+1+9 = 60. SSE (residual sum of squares) = 1+1+16+1+1 = 20. SSR (model sum of squares) = 60-20 = 40. R-squared = 40/60 = 0.667. About two-thirds of the cost variation is explained by prior hospitalizations; one-third (patient 3's high cost) is unexplained by the linear fit.

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.

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")