29  Misc. Multivariate Topics

---

The remaining multivariate topics did not fit cleanly into the preceding chapters.

This chapter asks four questions.

• How do we predict, rather than describe or explain?
  • Prediction interval: this extends the interval for one new observation to each value of \(x\).
  • Cross-validation: this reuses the out-of-sample score from Simple Regression to choose a bandwidth.
• Which observations can we use to estimate a time series at \(t\)?
  • Filtering and smoothing: filtering uses only current and past observations, while smoothing also uses future ones.
• What does lm() compute with several explanatory variables?
  • Matrix OLS: this is the formula \(\hat{\mathbf{b}}=(\hat{\textbf{X}}'\hat{\textbf{X}})^{-1}\hat{\textbf{X}}'\hat{y}\). The same algebra shows when OLS is unbiased.
• How do we avoid false discoveries when we run many tests?
  • Family-wise error rate: this is the chance of at least one false positive across \(K\) tests.
  • Bonferroni and Holm corrections: these make each test stricter to hold that rate down.
  • Post-hoc pairwise tests: these find which groups differ after a significant ANOVA.

29.1 Prediction

Describe vs. Explain vs. Predict

Understanding whether we aim to describe, explain, or predict is central to empirical analysis. The three distinct purposes are to accurately describe empirical patterns, explain why the empirical patterns exist, and predict what empirical patterns will exist.

Three Objectives of Empirical Analysis

Objective	Core Question	Goal	Methods
Describe	What is happening?	Characterize patterns & facts	Summary stats, visualizations, correlations
Explain	Why is it happening?	Establish causal relationships & mechanisms	Theoretical models; experimentation; counterfactual reasoning
Predict	What will happen?	Anticipate outcomes; support decisions that affect the future	Machine learning, treatment effect forecasting, policy simulations

Here is an example for minimum wages.

• Describe: Studies document that following minimum-wage increases, overall low-wage employment may look roughly stable in the short run. But disaggregated data often show larger employment declines over longer horizons, especially among youths and racial minorities.
• Explain: Studies investigate economic mechanisms. One question is whether lower productivity populations are more vulnerable and employers adjust along hiring margins (fewer openings, higher required skills). Another is whether effects are larger in sectors like retail/food service, where minimum wages bite hardest.
• Predict: A policy simulation that raises the wage floor incorporates subgroup-specific elasticities to forecast different employment losses (e.g., unemployment for teenage or black workers and income gains for stayers)

The distinctions matter, as they help you recognize that predictive accuracy or good data visualization does not equal causal insight (Shmueli 2010). You can then better align statistical tools with questions. Try remembering this: Describe reality, Explain causes, Predict futures

NoteMust Know

Consider the relationship between class size and student test scores. Write one sentence for each: (a) a descriptive finding, (b) an explanatory hypothesis, and (c) a predictive statement. Which of the three requires the strongest assumptions?

Prediction Intervals

A confidence interval for the conditional expectation \(m_{Y}(x)=\mathbb{E}[Y_{i}\mid X_{i}=x]\) does not by itself tell us how far a new observation at \(x\) might land.

ImportantKey Definition

A prediction interval for \(Y_{i}\) at \(X_{i} = x\) is a range that captures a stated probability (e.g., 95%) of new observations at that \(x\).

Prediction intervals are useful when the question is about individuals rather than averages. We ask what wage a new worker with 12 years of schooling might actually earn, not what the average wage of such workers is. Because they also include the residual variability of \(Y_{i}\) around the mean, they are generally wider than a confidence interval for \(m_{Y}(x)\) at the same level.

For a nice overview of different types of intervals, see https://www.jstor.org/stable/2685212. For an in-depth view, see “Statistical Intervals: A Guide for Practitioners and Researchers” or “Statistical Tolerance Regions: Theory, Applications, and Computation”. See https://robjhyndman.com/hyndsight/intervals/ for constructing intervals for future observations in a time-series context. See Davison and Hinkley, chapters 5 and 6 (also Efron and Tibshirani, or Wehrens et al.)

## MAIN EXAMPLE

## Data
# 3294 US workers, with hourly wages in 1980 dollars
data(Wages1, package='Ecdat')
xy_hat <- Wages1[, c('school', 'wage')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage
xy_hat <- xy_hat[order(xy_hat[, 'x_hat']), ]


# Design Points
x_values <- unique(xy_hat[, 'x_hat'])
# Regression
reg_lo <- loess(y_hat ~ x_hat, data=xy_hat, span=.8)
y_pred_lo <- predict(reg_lo, newdata=data.frame(x_hat=x_values))

# Bootstrap Residuals
n <- nrow(xy_hat)
res_lo <- lapply(1:399, function(i){
    xy_boot <- xy_hat[sample(n, replace=TRUE), ]
    reg_b <- loess(y_hat ~ x_hat, data=xy_boot, span=.8)
    resids_b <- resid(reg_b)
    cbind(e_hat=resids_b, x_hat=xy_boot[, 'x_hat'])
})
res_lo <- do.call(rbind, res_lo)

## Smooth residuals
res_fun <- function(x0, h){
    # Include only residuals within h distance to x0
    ki0 <- abs(res_lo[, 'x_hat']-x0)/h
    ki0 <- dunif(ki0)/h
    e_local_hat <- res_lo[ki0 != 0, 'e_hat'] #subset estimates
    ## Local quantiles
    res_i <- quantile(e_local_hat, probs=c(.025, .975), na.rm=TRUE)
    res_i
}
res_ci <- sapply(x_values, res_fun, h=2)

## Prediction Interval
pi_lo1 <- y_pred_lo + res_ci[1, ]
pi_lo2 <- y_pred_lo + res_ci[2, ]

#Plot
plot(y_hat ~ x_hat, pch=16, col=grey(0, .1),
    data=xy_hat, ylim=c(0, 40), main=NA,
    xlab='Years of schooling', ylab='Hourly wage (1980 dollars)')
lines(x_values, y_pred_lo,
    col=rgb(0, .7, .7, .8), lwd=2, type='o')
polygon( c(x_values, rev(x_values)),
    c(pi_lo1, rev(pi_lo2)),
    col=rgb(0, .7, .7, .3), border=NA)

TipTest Yourself

Consider also this example below with a linear model. Expand on it by using the Normal Approximation for the residuals (hint: use qnorm).

reg <- lm(y_hat ~ x_hat, data=xy_hat)

# Bootstrap Prediction Interval
boot_resids <- lapply(1:399, function(b){
    b_id <- sample( nrow(xy_hat), replace=TRUE)
    xy_boot <- xy_hat[b_id, ]
    reg_b <- lm(y_hat ~ x_hat, data=xy_boot)
    e_b_hat <- resid(reg_b)
    x_boot <- reg_b[['model']][['x_hat']]
    res_b <- cbind(e_b_hat, x_boot)
})
boot_resids <- as.data.frame(do.call(rbind, boot_resids))

# Errors do not vary with x_hat
q_e_hat <- quantile(boot_resids[, 'e_b_hat'], probs=c(.025, .975))

# PI
y_fit <- predict(reg)
boot_pi <- cbind(y_fit + q_e_hat[1], y_fit + q_e_hat[2])

# Plot Bootstrap PI
plot(y_hat ~ x_hat, data=xy_hat, pch=16, col=grey(0, .1), main=NA,
    ylim=c(0, 40), xlab='Years of schooling', ylab='Hourly wage (1980 dollars)')
title('Prediction Intervals', font.main=1)
x_hat <- xy_hat[, 'x_hat']
lines(x_hat, y_fit, col=rgb(1, 0, 0, .8), lwd=2)
polygon( c(x_hat, rev(x_hat)),
    c(boot_pi[, 1], rev(boot_pi[, 2])),
    col=rgb(1, 0, 0, .3), border=NA)

# Normal PI (Make by hand)
#pi <- predict(reg, interval='prediction')
#lines( x_hat, pi[, 'lwr'], lty=2)
#lines( x_hat, pi[, 'upr'], lty=2)

TipTest Yourself

Here is the same prediction interval for monthly earnings by years of education.

## Data
# Monthly earnings of 935 US working men in 1980
data(wage2, package='wooldridge')
xy2_hat <- wage2[, c('educ', 'wage')]
colnames(xy2_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=monthly earnings
xy2_hat <- xy2_hat[order(xy2_hat[, 'x_hat']), ]


# Design Points
x_values <- unique(xy2_hat[, 'x_hat'])
# Regression
reg_lo <- loess(y_hat ~ x_hat, data=xy2_hat, span=.8)
y_pred_lo <- predict(reg_lo, newdata=data.frame(x_hat=x_values))

# Bootstrap Residuals
n2 <- nrow(xy2_hat)
res_lo <- lapply(1:399, function(i){
    xy2_boot <- xy2_hat[sample(n2, replace=TRUE), ]
    reg_b <- loess(y_hat ~ x_hat, data=xy2_boot, span=.8)
    resids_b <- resid(reg_b)
    cbind(e_hat=resids_b, x_hat=xy2_boot[, 'x_hat'])
})
res_lo <- do.call(rbind, res_lo)

## Smooth residuals
res_fun <- function(x0, h){
    # Include only residuals within h distance to x0
    ki0 <- abs(res_lo[, 'x_hat']-x0)/h
    ki0 <- dunif(ki0)/h
    e_local_hat <- res_lo[ki0 != 0, 'e_hat'] #subset estimates
    ## Local quantiles
    res_i <- quantile(e_local_hat, probs=c(.025, .975), na.rm=TRUE)
    res_i
}
res_ci <- sapply(x_values, res_fun, h=2)

## Prediction Interval
pi_lo1 <- y_pred_lo + res_ci[1, ]
pi_lo2 <- y_pred_lo + res_ci[2, ]

#Plot
plot(y_hat ~ x_hat, pch=16, col=grey(0, .25),
    data=xy2_hat, ylim=c(0, 3100), main=NA, xlab='educ', ylab='wage')
lines(x_values, y_pred_lo,
    col=rgb(0, .7, .7, .8), lwd=2, type='o')
polygon( c(x_values, rev(x_values)),
    c(pi_lo1, rev(pi_lo2)),
    col=rgb(0, .7, .7, .3), border=NA)

Cross Validation

A model that fits the training data perfectly can still predict new observations badly, so we need a score that rewards out-of-sample accuracy rather than in-sample fit.

ImportantKey Definition

Leave-one-out cross-validation (LOO-CV) fits the model \(n\) times, each time leaving one observation out and predicting it from the rest, then chooses the bandwidth \(h\) that minimizes the average squared prediction error, \[\min_{h} \quad \frac{1}{n} \sum_{i=1}^{n} \left[ \hat{y}_{i} - \hat{y}^{\mathrm{pred}(-i)}(\hat{x}_{i},h) \right]^2,\] where \(\hat{y}^{\mathrm{pred}(-i)}(\hat{x}_{i},h)\) is the model’s predicted value at \(\hat{x}_{i}\) from a fit that excluded \(\hat{x}_{i}\) and \(h\) is the bandwidth (e.g., bin size in a regressogram).

Leave-one-out cross-validation is useful as the simplest computational way to choose bandwidths, and it directly addresses an issue that plagues observational studies in the social sciences: your model explains everything and predicts nothing. The score it minimizes is an out-of-sample average, so a too-flexible model that fits in-sample noise will be penalized. Minimizing prediction error is not objectively “best”, but it is a sensible choice when you have no better criterion.

In practice, compare models with cross-validation (below). The preferred model is the one with best out-of-sample performance, not the one with the most flexible in-sample fit (Stone 1974; Hastie et al. 2009; Mullainathan and Spiess 2017). The same idea selects bin counts and bandwidths in Local Relationships.

## MAIN EXAMPLE
data(Wages1, package='Ecdat')
# Crossvalidated bandwidth for regression
xy_hat <- Wages1[, c('school', 'wage')]
colnames(xy_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=hourly wage

## Grid Search
# The grid starts at 0.75.
# With a span of 0.5, only workers with exactly 12 years of schooling get weight at 12 years.
BWS <- c(0.75, 1, 1.5, 2, 3)
BWS_CV <- sapply(BWS, function(h){
    e_bw_hat <- sapply(1:nrow(xy_hat), function(i){
        llls <- loess(y_hat ~ x_hat, data=xy_hat[-i, ], span=h,
            degree=1)
        y_pred_i <- predict(llls, newdata=xy_hat[i, ])
        e_i_hat <-  (y_pred_i- xy_hat[i, 'y_hat'])
        return(e_i_hat)
    })
    # LOESS does not extrapolate, so it cannot predict the only worker with 3 years of schooling
    return( mean(e_bw_hat^2, na.rm=TRUE) )
})

## Plot MSE
par(mfrow=c(1, 2))
plot(BWS, BWS_CV, ylab='CV error', pch=3,
    xlab='bandwidth (h)')
# The error is smallest at the smallest bandwidth.
# It rises toward that of a linear regression as the bandwidth grows.
# So a more flexible fit predicts hourly wages slightly better.

## Plot Optimal Predictions
h_star <- BWS[which.min(BWS_CV)]
llls <- loess(y_hat ~ x_hat, data=xy_hat, span=h_star,
    degree=1)
plot(xy_hat, pch=16, col=grey(0, .1),
    xlab='Years of schooling', ylab='Hourly wage (1980 dollars)', main=NA)

pred_id <- order(xy_hat[, 'x_hat'])
lines(xy_hat[pred_id, 'x_hat'], predict(llls)[pred_id],
    col=rgb(0, .7, .7, .8), lwd=2)

TipTest Yourself

Here is the same bandwidth search for monthly earnings by years of education.

# Monthly earnings of 935 US working men in 1980
data(wage2, package='wooldridge')
# Crossvalidated bandwidth for regression
xy2_hat <- wage2[, c('educ', 'wage')]
colnames(xy2_hat) <- c('x_hat', 'y_hat') # x_hat=years of schooling, y_hat=monthly earnings

## Grid Search
BWS <- c(0.5, 0.75, 1, 1.5, 2, 3)
BWS_CV <- sapply(BWS, function(h){
    e_bw_hat <- sapply(1:nrow(xy2_hat), function(i){
        llls <- loess(y_hat ~ x_hat, data=xy2_hat[-i, ], span=h,
            degree=1, surface='direct')
        y_pred_i <- predict(llls, newdata=xy2_hat[i, ])
        e_i_hat <-  (y_pred_i- xy2_hat[i, 'y_hat'])
        return(e_i_hat)
    })
    return( mean(e_bw_hat^2) )
})

## Plot MSE
par(mfrow=c(1, 2))
plot(BWS, BWS_CV, ylab='CV error', pch=3,
    xlab='bandwidth (h)')
# The error falls as the bandwidth grows.
# It flattens toward that of a linear regression.
# So more flexibility does not improve predictions of monthly earnings.

## Plot Optimal Predictions
h_star <- BWS[which.min(BWS_CV)]
llls <- loess(y_hat ~ x_hat, data=xy2_hat, span=h_star,
    degree=1, surface='direct')
plot(xy2_hat, pch=16, col=grey(0, .25),
    xlab='X', ylab='Predictions', main=NA)

pred_id <- order(xy2_hat[, 'x_hat'])
lines(xy2_hat[pred_id, 'x_hat'], predict(llls)[pred_id],
    col=rgb(0, .7, .7, .8), lwd=2)

NoteMust Know

Leave-one-out cross-validation is just an average of out-of-sample errors. Suppose \(n=4\) and, for one bandwidth \(h\), leaving out each point in turn gives the predictions \(\hat{y}^{\mathrm{pred}(-i)}(\hat{x}_{i})\) below.

Leave-One-Out Prediction Errors

\(i\)	\(\hat{y}_{i}\)	\(\hat{y}^{\mathrm{pred}(-i)}(\hat{x}_{i})\)	\(\hat{y}_{i} - \hat{y}^{\mathrm{pred}(-i)}(\hat{x}_{i})\)
1	10	9	1
2	12	14	-2
3	15	14	1
4	11	11	0

The cross-validation score is the mean squared error: \((1^2 + (-2)^2 + 1^2 + 0^2)/4 = 6/4 = 1.5\). Repeating this for several bandwidths, the one with the smallest score is preferred.

29.2 Filtering

Filtering

In some cases we want to estimate the level of a time series at \(t\) rather than forecast forward, and the relevant question is which neighbors we are allowed to use.

ImportantKey Definition

Filtering estimates the value of a time series at \(t\) using only current and past observations, \(\mathbb{E}[Y_{t} \mid X_{t}, X_{t-1}]\); smoothing uses past, current, and future observations, \(\mathbb{E}[Y_{t} \mid X_{t-1}, X_{t}, X_{t+1}]\).

Filtering is useful for forecasting, or any time the future is unknown at time \(t\). Smoothing gives lower-variance estimates, so use it to analyze the past when you have the full series. The two are otherwise built from the same kernel weights, so the only practical difference is which observations the weights are allowed to touch.

One example of filtering is Exponential Filtering (sometimes confusingly referred to as “Exponential Smoothing”) which weights only current and previous observations using an exponential kernel.

## MAIN EXAMPLE
# Time series data
## Underlying Trend
trend_sim <- cumsum(rnorm(500, 0, 1))
## Observed Datapoints
x_sim <- trend_sim + runif(length(trend_sim), -10, 10)
xy_sim <- data.frame(t=seq_along(x_sim), trend_sim=trend_sim, x_sim=x_sim)

## Asymmetric kernel (filtering: uses current and past only)
bw1 <- c(2/3, 1/3)
s1 <- filter(x_sim, bw1/sum(bw1), sides=1)

## Symmetric kernel (smoothing: uses past, current, and future)
bw2 <- c(1/6, 2/3, 1/6)
s2 <- filter(x_sim, bw2/sum(bw2), sides=2)

## Plot the series with both smoothers
plot(x_sim ~ t, data=xy_sim, pch=16, col=grey(0, .25),
    main=NA, xlab='Time', ylab='x')
lines(xy_sim$t, s1, col=rgb(1, 0, 0, .8), lwd=2)
lines(xy_sim$t, s2, col=rgb(0, 0, 1, .8), lwd=2)
legend('topleft', lwd=2, bty='n',
    col=c(rgb(1, 0, 0, .8), rgb(0, 0, 1, .8)),
    legend=c('Asymmetric (filtering)', 'Symmetric (smoothing)'))

There are several cross-validation procedures for filtering time series data (Opsomer et al. 2001). One is called time series cross-validation (TSCV), which is useful for temporally dependent data (Hart 1991; Haerdle and Vieu 1992; Bergmeir et al. 2018). Unlike ordinary cross-validation, TSCV only ever predicts a point from its past, never its future, which respects the time ordering of the data. We use it below to choose the window of a one-sided moving average.

# Time series cross-validation for a one-sided moving average
# For each window h, predict each point from the previous h points
x_sim <- xy_sim$x_sim
ma_windows <- 2:20
tscv_mse <- sapply(ma_windows, function(h){
    ## Weight the previous h points equally, 0 weight on the current point
    ma_coefs <- c(0, rep(1/h, h))
    x_fit <- filter(x_sim, ma_coefs, sides=1)
    mean((x_sim - x_fit)^2, na.rm=TRUE)
})

## CV-optimal window
h_star <- ma_windows[which.min(tscv_mse)]

par(mfrow=c(1, 2))
## Plot CV error against window size
plot(ma_windows, tscv_mse, type='o', pch=16, main=NA,
    xlab='Moving-average window (h)', ylab='TSCV mean squared error')
abline(v=h_star, col=rgb(1, 0, 0, .8), lwd=2)

## Plot the series with the CV-optimal filter
x_smooth <- filter(x_sim, c(0, rep(1/h_star, h_star)), sides=1)
plot(x_sim ~ t, data=xy_sim, pch=16, col=grey(0, .25), main=NA,
    xlab='Time', ylab='x')
lines(xy_sim$t, x_smooth, col=rgb(1, 0, 0, .8), lwd=2)
title(paste0('One-sided MA, h = ', h_star), font.main=1)

## See also https://robjhyndman.com/hyndsight/tscvexample/

29.3 Inference

Arrays

Arrays are a generalization of matrices to multiple dimensions. They are a very efficient way to store well-formatted numeric data, and are often used in spatial econometrics and time series (often in the form of “data cubes”).

# data square (matrix)
array(data = seq(1, 24), dim=c(3, 8))
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    1    4    7   10   13   16   19   22
## [2,]    2    5    8   11   14   17   20   23
## [3,]    3    6    9   12   15   18   21   24

# data cube
a <- array(data = seq(1, 24), dim=c(3, 2, 4))
a
## , , 1
## 
##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    7   10
## [2,]    8   11
## [3,]    9   12
## 
## , , 3
## 
##      [,1] [,2]
## [1,]   13   16
## [2,]   14   17
## [3,]   15   18
## 
## , , 4
## 
##      [,1] [,2]
## [1,]   19   22
## [2,]   20   23
## [3,]   21   24

a[1, , , drop = FALSE]  # Row 1
#a[, 1, , drop = FALSE]  # Column 1
#a[, , 1, drop = FALSE]  # Layer 1

a[ 1, 1,  ]  # Row 1, column 1
#a[ 1,  , 1]  # Row 1, 'layer' 1
#a[  , 1, 1]  # Column 1, 'layer' 1
a[1 , 1, 1]  # Row 1, column 1, 'layer' 1

Apply extends to arrays

apply(a, 1, mean)    # Row means
## [1] 11.5 12.5 13.5
apply(a, 2, mean)    # Column means
## [1] 11 14
apply(a, 3, mean)    # 'Layer' means
## [1]  3.5  9.5 15.5 21.5
apply(a, c(1, 2), mean)  # Row/Column combination 
##      [,1] [,2]
## [1,]   10   13
## [2,]   11   14
## [3,]   12   15

The outer function produces an array by multiplying every element of one vector with every element of another.

x <- c(1, 2, 3)
x_outer <- outer(x, x) # x %o% x
x_outer
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    4    6
## [3,]    3    6    9
is.array(x_outer) # Matrices are arrays
## [1] TRUE

x_mat <- matrix( seq(6, 1), 2, 3)
outer(x_mat, x)
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]    6    4    2
## [2,]    5    3    1
## 
## , , 2
## 
##      [,1] [,2] [,3]
## [1,]   12    8    4
## [2,]   10    6    2
## 
## , , 3
## 
##      [,1] [,2] [,3]
## [1,]   18   12    6
## [2,]   15    9    3
# outer(x_mat, matrix(x))
# outer(x_mat, t(x))
# outer(x_outer, x_mat)

Matrices are the two-dimensional case, and they carry the regression algebra that follows.

Matrix Calculations

First, note that you can use matrix algebra in R

x_mat1 <- matrix( seq(2, 7), 2, 3)
x_mat1
##      [,1] [,2] [,3]
## [1,]    2    4    6
## [2,]    3    5    7

x_mat2 <- matrix( seq(4, -1), 2, 3)
x_mat2
##      [,1] [,2] [,3]
## [1,]    4    2    0
## [2,]    3    1   -1

x_mat1*x_mat2 #NOT classical matrix multiplication
##      [,1] [,2] [,3]
## [1,]    8    8    0
## [2,]    9    5   -7
x_mat1^x_mat2
##      [,1] [,2]      [,3]
## [1,]   16   16 1.0000000
## [2,]   27    5 0.1428571

tcrossprod(x_mat1, x_mat2) #x_mat1 %*% t(x_mat2)
##      [,1] [,2]
## [1,]   16    4
## [2,]   22    7

crossprod(x_mat1, x_mat2)
##      [,1] [,2] [,3]
## [1,]   17    7   -3
## [2,]   31   13   -5
## [3,]   45   19   -7

TipTest Yourself

Matrix multiplication %*% follows the rule that element \((i,j)\) of the product equals the dot product of row \(i\) of the first matrix and column \(j\) of the second. For crossprod(A, B), R computes \(A'B\) efficiently without explicitly transposing. This is useful in OLS where we need \((\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}\) repeatedly.

Second, note that you can apply functions to matrices

sum_squared <- function(x1, x2) {
    y <- (x1 + x2)^2
    return(y)
}
sum_squared(x_mat1, x_mat2)
##      [,1] [,2] [,3]
## [1,]   36   36   36
## [2,]   36   36   36


solve(x_mat1[1:2, 1:2])
##      [,1] [,2]
## [1,] -2.5    2
## [2,]  1.5   -1

Third, note that we can conduct OLS regressions efficiently using matrix algebra. Grouping the trial coefficients as a vector \(\mathbf{b}=(b_{0} ~~ b_{1} ~~... ~~ b_{K})\), we want to find coefficients that minimize the sum of squared errors. The linear model has a fairly simple solution for \(\hat{\mathbf{b}}\) if you know linear algebra. Denoting \(\hat{\mathbf{X}}_{i} = [1~~ \hat{x}_{i1} ~~...~~ \hat{x}_{iK}]\) as a row vector, we can write the model as \(\hat{y}_{i} = \hat{\mathbf{X}}_{i}\mathbf{b} + \hat{e}_{i}\). The residuals \(\hat{e}_{i}\) depend on the trial coefficients but are computed from the observed sample. We can then write the model in matrix form \[\begin{aligned} \hat{y} &= \hat{\textbf{X}}\mathbf{b} + \hat{\mathbf{e}} \\ \hat{y} &= \begin{pmatrix} \hat{y}_{1} \\ \vdots \\ \hat{y}_{n} \end{pmatrix} \quad \hat{\textbf{X}} = \begin{pmatrix} 1 & \hat{x}_{11} & ... & \hat{x}_{1K} \\ & \vdots & & \\ 1 & \hat{x}_{n1} & ... & \hat{x}_{nK} \end{pmatrix} \end{aligned}\] Minimizing the squared errors: \(\min_{\mathbf{b}} \sum_{i} \hat{e}_{i}^2 = \min_{\mathbf{b}} (\hat{\mathbf{e}}' \hat{\mathbf{e}})\), yields coefficient estimates and predictions \[\begin{aligned} \hat{\mathbf{b}} &= (\hat{\textbf{X}}'\hat{\textbf{X}})^{-1}\hat{\textbf{X}}'\hat{y} \\ \hat{y}^{\mathrm{fit}} &= \hat{\textbf{X}} \hat{\mathbf{b}} \\ \hat{\mathbf{e}} &= \hat{y} - \hat{y}^{\mathrm{fit}} \end{aligned}\]

## MAIN EXAMPLE
# Multivariate data: all columns of Wages1, not only schooling and wage as above
data(Wages1, package='Ecdat')
xy_hat <- Wages1
# Three-level schooling factor, as in Comparing Multiple Groups
xy_hat$school_g <- cut(xy_hat$school, breaks=c(0, 11, 12, 16), labels=c('<12', '12', '>12'))

# Manually compute b_hat = (X'X)^{-1} X'y_hat step by step.
y_hat <- xy_hat[, 'wage']
X <- xy_hat[, c('school', 'exper')]
X <- as.matrix(cbind(1, X))

XtX  <- t(X) %*% X         # K+1 by K+1 cross-product
XtXi <- solve(XtX)         # invert
XtY  <- t(X) %*% y_hat          # K+1 by 1 cross-product with y_hat
b_hat <- XtXi %*% XtY        # K+1 coefficient vector
c(b_hat)
## [1] -2.4766831  0.5992289  0.1572607

# Confirms the same answer as lm()
coef(lm(wage ~ school + exper, data=xy_hat))
## (Intercept)      school       exper 
##  -2.4766831   0.5992289   0.1572607

Bias

With matrix notation, we can establish that OLS is an unbiased estimator for additively separable and linear relationships. Assume \[\begin{aligned} Y_{i} = X_{i} \beta + \epsilon_{i} &\quad \mathbb{E}[\epsilon_{i} \mid X_{i} ]=0, \end{aligned}\] then, stacking the \(Y_{i}\) into the vector \(\mathbf{Y}\) and the rows \(X_{i}\) into the matrix \(\mathbf{X}\), notice that \[\begin{aligned} \mathbf{B} &= (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'(\mathbf{X}\beta + \epsilon) = \beta + (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\epsilon \\ \mathbb{E}\left[ \mathbf{B} \right] &= \mathbb{E}\left[ (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} \right] = \beta + (\mathbf{X}'\mathbf{X})^{-1}\mathbb{E}\left[ \mathbf{X}'\epsilon \right] = \beta \end{aligned}\]

If we do not know that the data generating process is additively separable and linear, then we do not know how to interpret \(\hat{\mathbf{b}}\). There is a very mathematical expression for it, which I do not cover here. If we run a local linear regression analysis, on the other hand, we get the estimator \[\begin{aligned} \beta(\mathbf{x}) &= [\mathbf{X}'\mathbf{K}(\mathbf{x})\mathbf{X}]^{-1} \mathbf{X}'\mathbf{K}(\mathbf{x})\mathbf{Y} \\ \mathbf{K}(\mathbf{x}) &= \begin{pmatrix} k\left(\mathbf{X}_{1}, \mathbf{x}, \mathbf{H}\right) & ... & 0\\ \vdots & & \\ 0 & ... & k\left(\mathbf{X}_{n}, \mathbf{x}, \mathbf{H}\right) \end{pmatrix}, \end{aligned}\] which is unbiased far more often.

The weight matrix \(\mathbf{K}(\mathbf{x})\) is \(n \times n\) and diagonal, holding one weight per observation rather than one per variable. Each weight multiplies the univariate kernels across the \(K\) variables, \[\begin{aligned} k\left(\mathbf{X}_{i}, \mathbf{x}, \mathbf{H}\right) &= \prod_{k=1}^{K} k\left(\frac{X_{ik}-x_{k}}{h_{k}}\right), \end{aligned}\] where \(h_{k}\) is the \(k\)th diagonal element of the bandwidth matrix \(\mathbf{H}\) from Local Relationships. So observation \(i\) gets a large weight only when it is close to \(\mathbf{x}\) in every variable, and \(h_{k}\) sets how close it has to be in variable \(k\).

29.4 Multiple Testing

When you test one hypothesis at the \(5\%\) level, there is a \(5\%\) chance of a false positive. When you test many hypotheses, those chances accumulate. This is the statistical engine behind much of Data Scientism: search long enough and something will look significant.

Family-Wise Error Rate

A single \(5\%\) test promises a \(5\%\) false-positive rate, but a researcher rarely stops at one test. The relevant question is the chance of any false positive across the whole batch.

ImportantKey Definition

The family-wise error rate (FWER) is the probability of at least one false positive across a family of \(K\) tests. If the tests are independent and each has false-positive rate \(\alpha\), \[\text{FWER} = 1-(1-\alpha)^{K}.\]

The FWER is useful for any procedure that runs many hypothesis tests. It tells you how likely you are to declare something “significant” by chance even when every null is true. With \(\alpha=0.05\), one test has FWER \(0.05\), but \(K=10\) tests give \(1-0.95^{10}\approx 0.40\) and \(K=100\) tests give \(1-0.95^{100}\approx 0.99\). Almost any dataset will yield a “significant” result if you test it enough ways.

## MAIN EXAMPLE
# Family-wise error rate grows with the number of tests
alpha <- 0.05
K <- c(1, 5, 10, 20, 50, 100)
fwer <- 1 - (1-alpha)^K
cbind(K, fwer=round(fwer, 3))
##        K  fwer
## [1,]   1 0.050
## [2,]   5 0.226
## [3,]  10 0.401
## [4,]  20 0.642
## [5,]  50 0.923
## [6,] 100 0.994

Bonferroni and Holm Corrections

To hold the family-wise error rate down to \(\alpha\), we have to make each individual test stricter. Exactly how strict matters for whether real effects survive.

ImportantKey Definition

The Bonferroni correction (Bonferroni 1936; Dunn 1961) multiplies each \(p\)-value by the number of tests \(K\) (equivalently, lowers the threshold to \(\alpha/K\)). The Holm correction applies the same idea stepwise: sort the \(p\)-values from smallest to largest and compare the \(j\)th smallest to \(\alpha/(K-j+1)\).

The Bonferroni correction is useful as the simplest valid adjustment, since it guarantees \(\text{FWER} \leq \alpha\) under any dependence structure between tests. It is conservative, though, and with many tests it can miss real effects. The Holm correction (Holm 1979) is uniformly more powerful than Bonferroni while still controlling the FWER, so it is generally preferred when only family-wise control is wanted.

## MAIN EXAMPLE
# Five p-values from five separate tests
pvals <- c(0.001, 0.013, 0.021, 0.04, 0.21)

# Adjust for multiple testing
cbind(
    raw        = pvals,
    bonferroni = p.adjust(pvals, method='bonferroni'),
    holm       = p.adjust(pvals, method='holm'))
##        raw bonferroni  holm
## [1,] 0.001      0.005 0.005
## [2,] 0.013      0.065 0.052
## [3,] 0.021      0.105 0.063
## [4,] 0.040      0.200 0.080
## [5,] 0.210      1.000 0.210

NoteMust Know

A researcher runs \(K=5\) tests and reports a \(p\)-value of \(0.013\). On its own, \(0.013 < 0.05\) looks significant. The Bonferroni-adjusted value is \(0.013 \times 5 = 0.065 > 0.05\), so it is not significant once we account for the five tests. The threshold view is equivalent: compare the raw \(p\)-value to \(\alpha/K = 0.05/5 = 0.01\).

Post-Hoc Pairwise Tests

A significant global test says some group differs, but the question we usually care about is which one. Asking each pair separately reintroduces the multiple-testing problem we just controlled at the global level.

ImportantKey Definition

A post-hoc test is a follow-up comparison conducted after a significant global test (like ANOVA). Because there are \(G(G-1)/2\) pairs to test, the per-pair significance threshold must be tightened to control the family-wise error rate; Tukey’s HSD (Tukey 1953) and Holm-adjusted pairwise \(t\)-tests are common choices.

Post-hoc tests are useful after a significant ANOVA \(F\)-test (see Comparing Multiple Groups). They find which group or groups differ without inflating the false-positive rate. The number of comparisons grows quadratically in \(G\), so the per-pair threshold gets noticeably stricter as soon as there are more than a handful of groups.

# Global test: do hourly wages differ by schooling group?
summary(aov(wage ~ school_g, data=xy_hat))
##               Df Sum Sq Mean Sq F value Pr(>F)    
## school_g       2   2541  1270.3     128 <2e-16 ***
## Residuals   3291  32654     9.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# F is about 128, with a p-value far below 5%, so some group differs.

## MAIN EXAMPLE
# Post-hoc: which schooling groups differ? (pairwise t-tests, Holm-adjusted)
pairwise.t.test(xy_hat$wage, xy_hat$school_g, p.adjust.method='holm')
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  xy_hat$wage and xy_hat$school_g 
## 
##     <12     12     
## 12  6.8e-08 -      
## >12 < 2e-16 < 2e-16
## 
## P value adjustment method: holm
# All three pairs differ at the 5% level.

# Tukey's Honest Significant Differences:
# pairwise confidence intervals that already account for all comparisons
TukeyHSD(aov(wage ~ school_g, data=xy_hat))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = wage ~ school_g, data = xy_hat)
## 
## $school_g
##              diff       lwr       upr p adj
## 12-<12  0.6774177 0.3837567 0.9710787 2e-07
## >12-<12 2.2849502 1.9491477 2.6207526 0e+00
## >12-12  1.6075324 1.2634931 1.9515718 0e+00
# No interval contains 0, so all three pairs differ here too.
# The largest gap is about 2.28 dollars an hour, between >12 and <12 years.

To see what the adjustment does, the chunk below puts the raw and adjusted \(p\)-values side by side.

# Raw p-values, one per pair of schooling groups
pair_tests <- pairwise.t.test(xy_hat$wage, xy_hat$school_g, p.adjust.method='none')
p_raw <- pair_tests[['p.value']]
p_raw <- p_raw[!is.na(p_raw)] # drop the empty cell of the table
names(p_raw) <- c('12 vs <12', '>12 vs <12', '>12 vs 12')

# Adjust for the three pairs
cbind(
    raw=p_raw,
    bonferroni=p.adjust(p_raw, method='bonferroni'),
    holm=p.adjust(p_raw, method='holm'))
##                     raw   bonferroni         holm
## 12 vs <12  6.794126e-08 2.038238e-07 6.794126e-08
## >12 vs <12 2.957116e-55 8.871348e-55 8.871348e-55
## >12 vs 12  1.843717e-27 5.531152e-27 3.687435e-27
# Bonferroni multiplies each raw p-value by 3.
# Holm multiplies the smallest by 3, the next by 2, and the largest by 1.
# All three stay far below 0.05, so the adjustment changes no conclusion here.

An adjustment never lowers a \(p\)-value. A raw \(p\)-value not far below the threshold, such as \(0.013\) in the five tests above, can lose its significance, while a much smaller one keeps it.

TipTest Yourself

Repeat the post-hoc tests with murder arrest rates across the four census regions.

# Murder arrests per 100,000 in 50 US states in 1973, with Census region
xy3_hat <- USArrests
xy3_hat$Region <- state.region

# Global test: do murder rates differ by region?
summary(aov(Murder ~ Region, data=xy3_hat))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Region       3  391.2   130.4   11.14 1.28e-05 ***
## Residuals   46  538.3    11.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Post-hoc: which regions differ? (pairwise t-tests, Holm-adjusted)
pairwise.t.test(xy3_hat$Murder, xy3_hat$Region, p.adjust.method='holm')
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  xy3_hat$Murder and xy3_hat$Region 
## 
##               Northeast South   North Central
## South         7e-05     -       -            
## North Central 0.67251   0.00017 -            
## West          0.36895   0.00259 0.67251      
## 
## P value adjustment method: holm

# Tukey's Honest Significant Differences:
# pairwise confidence intervals that already account for all comparisons
TukeyHSD(aov(Murder ~ Region, data=xy3_hat))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Murder ~ Region, data = xy3_hat)
## 
## $Region
##                              diff       lwr       upr     p adj
## South-Northeast          7.006250  3.206920 10.805580 0.0000673
## North Central-Northeast  1.000000 -3.020833  5.020833 0.9104265
## West-Northeast           2.330769 -1.623231  6.284769 0.4047789
## North Central-South     -6.006250 -9.488393 -2.524107 0.0001909
## West-South              -4.675481 -8.080233 -1.270728 0.0035034
## West-North Central       1.330769 -2.319509  4.981047 0.7660715

Rerun pairwise.t.test on the region data with p.adjust.method='none' and again with 'bonferroni'. How many pairs are significant at the \(5\%\) level under each? Why can the uncorrected version flag more differences than the corrected ones, and which version should you report?

29.5 Exercises

1. Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.

2. Explain the distinction between description, explanation, and prediction using a concrete example other than minimum wages. For each objective, state what statistical tools you would use and what kind of evidence would be needed.

3. Using the USArrests dataset, regress Murder on Assault and UrbanPop. Compute \(\hat{\mathbf{b}}\) both with lm() and manually via \((\hat{\mathbf{X}}'\hat{\mathbf{X}})^{-1}\hat{\mathbf{X}}'\hat{y}\) using matrix operations in R. Verify that the two approaches give the same coefficients.

4. Using the wage2 dataset from the wooldridge package, fit a loess regression of wage on educ for three different span values (e.g., 0.4, 0.6, 0.8). Use leave-one-out cross-validation to compute the mean squared prediction error for each span. Which span minimizes prediction error, and does a more flexible fit always predict better?

5. A researcher runs \(K=20\) independent tests, each at the \(5\%\) level. Compute the family-wise error rate \(1-(1-0.05)^{20}\), and find the per-test Bonferroni threshold \(\alpha/K\) that holds the family-wise rate at \(5\%\). Then create a vector of 20 \(p\)-values, apply p.adjust() with the 'bonferroni' and 'holm' methods, and report how many results remain significant under each.

Further Reading

• https://otexts.com/fpp3/ – Hyndman and Athanasopoulos, Forecasting: Principles and Practice (3rd ed.), covering time-series smoothing, cross-validation, and prediction in R.
• https://cran.r-project.org/web/packages/np/vignettes/np.pdf – the np package vignette for nonparametric kernel estimation and bandwidth selection.
• Shmueli (2010) – a clear discussion of the describe/explain/predict distinction in statistical modeling.
• https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf – worked examples of multiple-comparison procedures and post-hoc tests in R.

Recall

This chapter wrapped up the multivariate part with the describe/explain/predict taxonomy, prediction intervals and cross-validation, time-series filtering, the matrix form of OLS, and multiple-testing corrections. The matrix-OLS sense-check on Wages1 made the algebra concrete: computing \(\hat{\mathbf{b}}=(\hat{\textbf{X}}'\hat{\textbf{X}})^{-1}\hat{\textbf{X}}'\hat{y}\) step by step with XtX <- t(X) %*% X, XtXi <- solve(XtX), and b_hat <- XtXi %*% XtY recovers exactly the coefficient vector that lm(wage ~ school + exper, data=xy_hat) returns. These tools connect back to earlier chapters. Cross-validation chooses models, and multiple-testing corrections address the search problem from Data Scientism. The matrix form writes out the OLS algebra of Multiple Regression.

Bergmeir, Christoph, Rob J. Hyndman, and Bonsoo Koo. 2018. “A Note on the Validity of Cross-Validation for Evaluating Autoregressive Time Series Prediction.” Computational Statistics & Data Analysis 120: 70–83. https://doi.org/10.1016/j.csda.2017.11.003.

Bonferroni, Carlo E. 1936. “Teoria statistica delle classi e calcolo delle probabilità.” Pubblicazioni del R. Istituto Superiore di Scienze Economiche e Commerciali di Firenze 8: 3–62.

Dunn, Olive Jean. 1961. “Multiple Comparisons Among Means.” Journal of the American Statistical Association 56 (293): 52–64. https://doi.org/10.1080/01621459.1961.10482090.

Haerdle, Wolfgang, and Philippe Vieu. 1992. “Kernel Regression Smoothing of Time Series.” Journal of Time Series Analysis 13 (3): 209–32. https://doi.org/10.1111/j.1467-9892.1992.tb00103.x.

Hart, Jeffrey D. 1991. “Kernel Regression Estimation with Time Series Errors.” Journal of the Royal Statistical Society. Series B (Methodological) 53 (1): 173–87.

Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer New York. https://doi.org/10.1007/978-0-387-84858-7_6.

Holm, Sture. 1979. “A Simple Sequentially Rejective Multiple Test Procedure.” Scandinavian Journal of Statistics 6: 65–70. https://www.jstor.org/stable/4615733.

Mullainathan, Sendhil, and Jann Spiess. 2017. “Machine Learning: An Applied Econometric Approach.” Journal of Economic Perspectives 31 (2): 87–106. https://doi.org/10.1257/jep.31.2.87.

Opsomer, Jean, Yuedong Wang, and Yuhong Yang. 2001. “Nonparametric Regression with Correlated Errors.” Statistical Science 16 (2): 134–53. http://www.jstor.org/stable/2676791.

Shmueli, Galit. 2010. “To Explain or to Predict?” Statistical Science 25 (3): 289–310. https://doi.org/10.1214/10-STS330.

Stone, M. 1974. “Cross-Validatory Choice and Assessment of Statistical Predictions.” Journal of the Royal Statistical Society: Series B (Methodological) 36 (2): 111–33. https://doi.org/10.1111/j.2517-6161.1974.tb00994.x.

Tukey, John W. 1953. “The Problem of Multiple Comparisons.” Unpublished manuscript.