Hypothesis Testing in Linear Regression

Wald tests, joint vs. individual significance, and nonlinear hypotheses

N. Muriel

Today’s plan

  • Block 1 — The general linear hypothesis \(R\beta = q\): Wald, F vs. \(\chi^2\)
  • Block 2 — Joint vs. individual significance
  • (15 min break)
  • Block 3 — Nonlinear hypotheses \(c(\beta) = q\): the delta method
  • Block 4 — Guided lab, all three together

Builds directly on last session’s lm() fundamentals.

Block 1 — The General Linear Hypothesis

From single tests to \(R\beta = q\)

A single-coefficient t-test — \(H_0: \beta_j = 0\) — is a special case of testing:

\[H_0: R\beta = q\]

  • \(R\): a \(q \times p\) matrix of restrictions
  • \(q\) (vector): the hypothesized values
  • Covers: one coefficient, several coefficients jointly, or linear combinations of coefficients

The data: NMES1988

Relevance: health utilization/expenditure data — the kind of variables that feed health insurance pricing.

data("NMES1988", package = "AER")
str(NMES1988)
'data.frame':   4406 obs. of  19 variables:
 $ visits   : int  5 1 13 16 3 17 9 3 1 0 ...
 $ nvisits  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ ovisits  : int  0 2 0 5 0 0 0 0 0 0 ...
 $ novisits : int  0 0 0 0 0 0 0 0 0 0 ...
 $ emergency: int  0 2 3 1 0 0 0 0 0 0 ...
 $ hospital : int  1 0 3 1 0 0 0 0 0 0 ...
 $ health   : Factor w/ 3 levels "poor","average",..: 2 2 1 1 2 1 2 2 2 2 ...
  ..- attr(*, "contrasts")= num [1:3, 1:2] 1 0 0 0 0 1
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "poor" "average" "excellent"
  .. .. ..$ : chr [1:2] "poor" "excellent"
 $ chronic  : int  2 2 4 2 2 5 0 0 0 0 ...
 $ adl      : Factor w/ 2 levels "normal","limited": 1 1 2 2 2 2 1 1 1 1 ...
 $ region   : Factor w/ 4 levels "northeast","midwest",..: 4 4 4 4 4 4 2 2 2 2 ...
  ..- attr(*, "contrasts")= num [1:4, 1:3] 1 0 0 0 0 1 0 0 0 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:4] "northeast" "midwest" "west" "other"
  .. .. ..$ : chr [1:3] "northeast" "midwest" "west"
 $ age      : num  6.9 7.4 6.6 7.6 7.9 6.6 7.5 8.7 7.3 7.8 ...
 $ afam     : Factor w/ 2 levels "no","yes": 2 1 2 1 1 1 1 1 1 1 ...
 $ gender   : Factor w/ 2 levels "female","male": 2 1 1 2 1 1 1 1 1 1 ...
 $ married  : Factor w/ 2 levels "no","yes": 2 2 1 2 2 1 1 1 1 1 ...
 $ school   : int  6 10 10 3 6 7 8 8 8 8 ...
 $ income   : num  2.881 2.748 0.653 0.659 0.659 ...
 $ employed : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...
 $ insurance: Factor w/ 2 levels "no","yes": 2 2 1 2 2 1 2 2 2 2 ...
 $ medicaid : Factor w/ 2 levels "no","yes": 1 1 2 1 1 2 1 1 1 1 ...

A baseline model

fit <- lm(visits ~ chronic + age + income + insurance + gender + married,
          data = NMES1988)
summary(fit)

Call:
lm(formula = visits ~ chronic + age + income + insurance + gender + 
    married, data = NMES1988)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.831  -3.769  -1.584   1.808  82.597 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.828464   1.237259   3.903 9.66e-05 ***
chronic       1.337049   0.072978  18.321  < 2e-16 ***
age          -0.245810   0.160580  -1.531   0.1259    
income        0.004867   0.034849   0.140   0.8889    
insuranceyes  1.350541   0.241377   5.595 2.34e-08 ***
gendermale   -0.533735   0.216784  -2.462   0.0139 *  
marriedyes   -0.259540   0.227981  -1.138   0.2550    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.496 on 4399 degrees of freedom
Multiple R-squared:  0.0776,    Adjusted R-squared:  0.07634 
F-statistic: 61.68 on 6 and 4399 DF,  p-value: < 2.2e-16

A joint restriction

\(H_0\): the coefficients on insurance and married are both zero.

linearHypothesis(fit, c("insuranceyes = 0", "marriedyes = 0"))

Linear hypothesis test:
insuranceyes = 0
marriedyes = 0

Model 1: restricted model
Model 2: visits ~ chronic + age + income + insurance + gender + married

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   4401 186960                                  
2   4399 185635  2    1325.3 15.703 1.602e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The Wald statistic

\[W = (R\hat\beta - q)' \left[R\, \widehat{\mathrm{Var}}(\hat\beta)\, R'\right]^{-1} (R\hat\beta - q)\]

  • \(\sigma^2\) estimated (our usual case) → \(W / q \sim F(q, n-p)\)
  • \(\sigma^2\) known, or asymptotically → \(W \sim \chi^2(q)\)

car::linearHypothesis() computes the same Wald statistic either way — the test= argument only changes the reference distribution.

F vs. Chi-squared, explicitly

linearHypothesis(fit, c("insuranceyes = 0", "marriedyes = 0"), test = "F")

Linear hypothesis test:
insuranceyes = 0
marriedyes = 0

Model 1: restricted model
Model 2: visits ~ chronic + age + income + insurance + gender + married

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   4401 186960                                  
2   4399 185635  2    1325.3 15.703 1.602e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
linearHypothesis(fit, c("insuranceyes = 0", "marriedyes = 0"), test = "Chisq")

Linear hypothesis test:
insuranceyes = 0
marriedyes = 0

Model 1: restricted model
Model 2: visits ~ chronic + age + income + insurance + gender + married

  Res.Df    RSS Df Sum of Sq  Chisq Pr(>Chisq)    
1   4401 186960                                   
2   4399 185635  2    1325.3 31.406  1.515e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Building \(R\) and \(W\) by hand

R <- matrix(0, nrow = 1, ncol = length(coef(fit)))
colnames(R) <- names(coef(fit))
R[1, "age"] <- 1
q <- 0

beta_hat <- coef(fit)
V <- vcov(fit)
W <- t(R %*% beta_hat - q) %*% solve(R %*% V %*% t(R)) %*% (R %*% beta_hat - q)
W
         [,1]
[1,] 2.343233

Comparing to F and \(\chi^2\) p-values

pf(W, df1 = 1, df2 = fit$df.residual, lower.tail = FALSE)
          [,1]
[1,] 0.1259001
pchisq(W, df = 1, lower.tail = FALSE)
          [,1]
[1,] 0.1258282
linearHypothesis(fit, "age = 0")

Linear hypothesis test:
age = 0

Model 1: restricted model
Model 2: visits ~ chronic + age + income + insurance + gender + married

  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1   4400 185734                           
2   4399 185635  1    98.883 2.3432 0.1259
summary(fit)$coefficients["age", ]
  Estimate Std. Error    t value   Pr(>|t|) 
-0.2458097  0.1605799 -1.5307621  0.1259001

For one restriction, \(t^2 = F\) — this single test is exactly what summary() already gave us.

Convergence as \(n\) grows

q_restr <- 2
n_vals <- c(30, 100, 1000, 100000)
crit_F     <- qf(0.95, df1 = q_restr, df2 = n_vals - 7) * q_restr
crit_Chisq <- qchisq(0.95, df = q_restr)
data.frame(n = n_vals, F_based = crit_F, Chisq_based = crit_Chisq)
      n  F_based Chisq_based
1 3e+01 6.844264    5.991465
2 1e+02 6.188675    5.991465
3 1e+03 6.009576    5.991465
4 1e+05 5.991644    5.991465

F-based critical values shrink toward the \(\chi^2\)-based one as \(n\) grows — the two tests agree asymptotically.

Block 2 — Joint vs. Individual Significance

The idea

Correlated regressors can inflate standard errors enough that each coefficient looks individually insignificant, while the group is jointly informative.

Not a contradiction — exactly what multicollinearity does to \(t\)-tests vs. \(F\)-tests.

A model with a correlated group

fit2 <- lm(visits ~ chronic + health + hospital + adl + age + income,
           data = NMES1988)
summary(fit2)

Call:
lm(formula = visits ~ chronic + health + hospital + adl + age + 
    income, data = NMES1988)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.870  -3.696  -1.458   1.865  73.624 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      7.04775    1.18368   5.954 2.82e-09 ***
chronic          0.94134    0.07831  12.021  < 2e-16 ***
healthpoor       1.40013    0.32247   4.342 1.44e-05 ***
healthexcellent -1.21563    0.36560  -3.325 0.000891 ***
hospital         1.62283    0.13402  12.109  < 2e-16 ***
adllimited       0.38400    0.26915   1.427 0.153741    
age             -0.46612    0.15999  -2.914 0.003592 ** 
income           0.03382    0.03314   1.020 0.307595    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.376 on 4398 degrees of freedom
Multiple R-squared:  0.1116,    Adjusted R-squared:  0.1102 
F-statistic: 78.94 on 7 and 4398 DF,  p-value: < 2.2e-16

Look closely at hospital and adl — how convincing do they look alone?

Checking the correlation

cor(data.frame(chronic = NMES1988$chronic,
               hospital = NMES1988$hospital,
               adl = as.numeric(NMES1988$adl)))
           chronic  hospital       adl
chronic  1.0000000 0.2335242 0.2553067
hospital 0.2335242 1.0000000 0.1591745
adl      0.2553067 0.1591745 1.0000000
vif(fit2)
             GVIF Df GVIF^(1/(2*Df))
chronic  1.210303  1        1.100138
health   1.271045  2        1.061794
hospital 1.084296  1        1.041295
adl      1.275171  1        1.129235
age      1.112730  1        1.054860
income   1.018153  1        1.009036

Finding the right coefficient name

Factor dummies are named by their levels, not a guess from the variable name — always check first:

names(coef(fit2))
[1] "(Intercept)"     "chronic"         "healthpoor"      "healthexcellent"
[5] "hospital"        "adllimited"      "age"             "income"         
adl_coef <- grep("^adl", names(coef(fit2)), value = TRUE)
adl_coef
[1] "adllimited"

The joint test

linearHypothesis(fit2, c("hospital = 0", paste0(adl_coef, " = 0")))

Linear hypothesis test:
hospital = 0
adllimited = 0

Model 1: restricted model
Model 2: visits ~ chronic + health + hospital + adl + age + income

  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   4400 184944                                  
2   4398 178789  2    6154.5 75.697 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

If individual \(p\)-values were high but this rejects \(H_0\): individually insignificant, jointly significant.

The same story via anova()

fit2_restricted <- lm(visits ~ chronic + age + income, data = NMES1988)
anova(fit2_restricted, fit2)
Analysis of Variance Table

Model 1: visits ~ chronic + age + income
Model 2: visits ~ chronic + health + hospital + adl + age + income
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   4402 187337                                  
2   4398 178789  4    8548.1 52.568 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Same \(F\)-statistic logic, framed as nested model comparison.

Actuarial takeaway

Don’t drop a correlated rating variable from a tariff model just because its own \(p\)-value is high.

Test the whole correlated block jointly before deciding it adds nothing.

Block 3 — Nonlinear Hypotheses

Motivation

Not every natural hypothesis is linear in \(\beta\).

Claim frequency often follows a U-shape in policyholder age:

\[E[\text{claims}] = \beta_0 + \beta_1\,\text{age} + \beta_2\,\text{age}^2\]

Risk-minimizing age: \(\text{age}^* = -\beta_1 / (2\beta_2)\)

\(H_0: \text{age}^* = 45\) is nonlinear in \((\beta_1, \beta_2)\) — linearHypothesis() doesn’t apply directly.

Simulated data, known turning point

set.seed(2026)
n <- 2000
age <- runif(n, 18, 80)
true_b0 <- 2; true_b1 <- -0.09; true_b2 <- 0.001  # turning point = 45
claims <- true_b0 + true_b1*age + true_b2*age^2 + rnorm(n, sd = 0.5)
sim_data <- data.frame(age, claims)

The simulated pattern

ggplot(sim_data, aes(age, claims)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = TRUE)) +
  theme_minimal()

Fitting the quadratic model

fit3 <- lm(claims ~ age + I(age^2), data = sim_data)
summary(fit3)

Call:
lm(formula = claims ~ age + I(age^2), data = sim_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.84386 -0.32160  0.00378  0.32748  1.40402 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.1506803  0.0869191   24.74   <2e-16 ***
age         -0.0968886  0.0038447  -25.20   <2e-16 ***
I(age^2)     0.0010642  0.0000387   27.50   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4917 on 1997 degrees of freedom
Multiple R-squared:  0.3114,    Adjusted R-squared:  0.3107 
F-statistic: 451.5 on 2 and 1997 DF,  p-value: < 2.2e-16

The delta method

\[\sqrt{n}\big(c(\hat\beta) - c(\beta)\big) \;\xrightarrow{d}\; N\!\left(0,\; \nabla c(\beta)'\, \Sigma\, \nabla c(\beta)\right)\]

car::deltaMethod() computes \(\nabla c\) (numerically or symbolically) and the resulting SE for you.

Estimating the turning point

deltaMethod(fit3, "-age / (2*`I(age^2)`)")
                      Estimate       SE    2.5 % 97.5 %
-age/(2 * `I(age^2)`) 45.52018  0.31701 44.89886 46.142

Estimate, delta-method SE, and a confidence interval — does it contain 45, the true value?

Testing \(H_0: \text{age}^* = 45\)

dm <- deltaMethod(fit3, "-age / (2*`I(age^2)`)")
est <- dm$Estimate
se  <- dm$SE
wald_stat <- ((est - 45) / se)^2
wald_stat
[1] 2.692654
pchisq(wald_stat, df = 1, lower.tail = FALSE)
[1] 0.1008118
pf(wald_stat, df1 = 1, df2 = fit3$df.residual, lower.tail = FALSE)
[1] 0.1009693

Alternatives

  • msm::deltamethod() — same idea, you supply the gradient formula explicitly (good for seeing the mechanics once)
  • Bootstrap — non-parametric fallback when \(c(\cdot)\) is messy or you distrust the normal approximation
boot_turning_point <- function(data, n_boot = 500) {
  est <- numeric(n_boot)
  for (i in seq_len(n_boot)) {
    idx <- sample(nrow(data), replace = TRUE)
    m <- lm(claims ~ age + I(age^2), data = data[idx, ])
    est[i] <- -coef(m)["age"] / (2 * coef(m)["I(age^2)"])
  }
  est
}
boot_est <- boot_turning_point(sim_data)
quantile(boot_est, c(0.025, 0.975))
    2.5%    97.5% 
44.81702 46.08511

Block 4 — Guided Lab

Your turn: NMES1988, all three ideas together

  1. Fit visits ~ chronic + age + I(age^2) + income + insurance + medicaid + gender. Report individual significance.
  2. Joint test: insurance and medicaid jointly zero — both test = "F" and test = "Chisq". Do the \(p\)-values agree?
  3. Find two correlated predictors (cor()/vif()) and test their joint significance. Individually insignificant but jointly significant?
  4. Nonlinear hypothesis: use deltaMethod() on the age/age^2 terms to estimate the visits-minimizing (or maximizing) age. Test it against a reference value.
  5. Write 3–4 sentences: how would an insurer use this in age-based rating factor design?

Wrap-up

  • \(R\beta = q\) unifies t-tests, F-tests, and arbitrary linear restrictions
  • test = "F" vs. "Chisq" — same Wald statistic, different reference distribution, converge as \(n \to \infty\)
  • Joint significance \(\neq\) sum of individual significances — multicollinearity is the usual culprit
  • Nonlinear hypotheses need the delta method (or bootstrap) — car::deltaMethod() does the work
  • All three show up directly in actuarial practice: tariff variable selection, rating factor justification, risk-curve turning points

Questions?