Diagnostics and Transformations

Problem 1

Given the following small dataset:

\[\mathbf{y} = \begin{bmatrix}5\\8\\12\\20\\30\end{bmatrix}\]

\[\mathbf{X} = \begin{bmatrix}1&1\\1&2\\1&3\\1&4\\1&5\end{bmatrix}\]

a) Fit the linear regression model \(y = \beta_0 + \beta_1 x + \varepsilon\) and calculate the residuals \(\hat{\varepsilon}_i\) for each observation.

b) Calculate the leverage \(h_i\) for each observation using the formula \(h_i = H_{ii}\) where \(\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\). You may use R for the matrix calculations.

c) Calculate the standardized residuals \(r_i = \frac{\hat{\varepsilon}_i}{\hat{\sigma}\sqrt{1-h_i}}\) for each observation.

d) Calculate the studentized residuals \(t_i = \frac{\hat{\varepsilon}_i}{\hat{\sigma}_{(i)}\sqrt{1-h_i}}\) for each observation. For this, you’ll need to fit the model 5 times, each time leaving out one observation, to get \(\hat{\sigma}_{(i)}\).

e) Which observation(s), if any, would you flag as potential outliers based on the studentized residuals?

Problem 2

Using the same data from Problem 1:

a) Calculate Cook’s distance for each observation using the formula: \[D_i = \frac{r_i^2}{p} \cdot \frac{h_i}{1-h_i}\] where \(p = 2\) (number of parameters) and \(r_i\) are the standardized residuals from Problem 1.

b) Would any observations be considered influential using the rule of thumb \(D_i > 4/n\)?

c) Refit the model excluding the most influential point. How much do \(\hat{\beta}_0\) and \(\hat{\beta}_1\) change? Does this confirm the Cook’s distance calculation?

Problem 3

Consider a dataset where the relationship between y and x is quadratic: \(y = 2 + 3x - x^2 + \varepsilon\) where \(\varepsilon \sim N(0, 2)\).

Generate data with n = 100 and x from a standard normal distribution (\(N(0,1)\)).

a) Fit a linear model lm(y ~ x) and create a partial regression plot. Since there’s only one predictor, this is just a scatterplot of y vs. x with the fitted line. Does the relationship look linear?

b) Create a partial regression plot for x (here there are not any other predictors so you can just plot x on the x-axis and the residuals for the model on the y-axis). What pattern do you see? What does this suggest about the model specification?

c) Now fit a quadratic model lm(y ~ x + I(x^2)). Create the partial regression plot for this model. How does it compare to part (b)?

d) Calculate the IQR of x. Using the quadratic model, find:

The predicted y at the 25th percentile of x
The predicted y at the 75th percentile of x
The change in predicted y for an IQR increase in x

e) Interpret this IQR change. Why is this more useful than trying to interpret individual coefficients in a polynomial model?

Problem 4

Using the same data generation as Problem 3 (quadratic relationship):

a) Fit four models:

Linear: lm(y ~ x)
Quadratic: lm(y ~ x + I(x^2))
Natural spline with 3 df: lm(y ~ ns(x, df = 3))
10th order Polynomial: lm(y ~ poly(x, 10))

Create a plot showing:

The data points
All four fitted curves
The true relationship (generate a grid of x values and calculate the true y without error)

Which model comes closest to the true relationship?

b) Calculate the \(R^2\) for all three models. Which has the highest \(R^2\)? Does this necessarily mean it’s the best model?

c) For the natural spline model, what are the knot locations? (Hint: use attr(ns(x, df = 3), "knots"))

Problem 5

Suppose we fit a linear regression on the log scale: log(y) ~ x, obtaining estimates \(\hat{\beta}_0\) and \(\hat{\beta}_1\).

a) Show that the naive back-transformation \(\hat{y} = \exp(\hat{\beta}_0 + \hat{\beta}_1 x)\) is a biased estimator of \(E[y \mid x]\).

Hint: Start by writing \(E[y \mid x]\) in terms of \(\beta_0\), \(\beta_1\), and \(E[\exp(\varepsilon)]\). What is \(E[\exp(\varepsilon)]\) when \(\varepsilon \sim N(0, \sigma^2)\)?

b) Derive the correct formula for \(E[y \mid x]\) that includes the necessary correction factor.

c) Write a simulation to verify your answer from part (b). We are going to generate data from the model \(y = \exp(2 + 1.5x + \varepsilon)\) where \(x \sim N(0,1)\), \(\varepsilon \sim N(0, 0.25)\). First, let’s get the true \(E[y \mid x = 0]\) by running:

n <- 10000
e <- rnorm(n, 0, 0.25)
y <- exp(2 + e)
true_y <- mean(y)

Now let’s simulate multiple data sets. Let’s have \(n = 200\). For each simulated dataset:

Fit the model log(y) ~ x
Compute a prediction at \(x = 0\) using the naive back-transformation (i.e. just exponentiate the predicted values)
Compute a prediction at \(x = 0\) using the corrected back-transformation from part (b)

Repeat this process 5000 times and compare the average of each prediction method to the true value of \(E[y \mid x = 0]\). Which estimator is closer to the truth?

d) Based on your simulation results, approximately what percentage bias does the naive approach have? How does this relate to \(\sigma^2\)?

e) Suppose \(\sigma^2\) were much larger, say \(\sigma^2 = 1\). Without running a new simulation, predict approximately what the percentage bias would be for the naive approach in this case.

Mathematical fact you may use: If \(Z \sim N(\mu, \sigma^2)\), then \(E[\exp(Z)] = \exp(\mu + \sigma^2/2)\).

Problem 6

This problem explores Type I error inflation when we test for non-linearity and add terms based on significance tests, compared to using flexible modeling from the start.

Simulation setup:

True model: \(y = \varepsilon\) where \(\varepsilon \sim N(0, 1)\) (there is no relationship between \(x\) and \(y\))
Sample size: n = 50
Generate x from \(N(0,1)\)
Number of simulations: 5000

Strategy 1: Test-then-add approach

Fit lm(y ~ x)
Test if adding \(x^2\) is significant at \(\alpha = 0.05\) by fitting lm(y ~ x + I(x^2)) and checking if the p-value for the \(x^2\) coefficient is < 0.05
If significant, use the quadratic model; otherwise use the linear model
Perform an F-test to test whether the model with the predictor(s) is better than the intercept only model. Use \(\alpha = 0.05\).
Record if you reject \(H_0\)

Strategy 2: Natural spline by default

Fit lm(y ~ ns(x, df = 3)) (always use this flexible model)
Perform an F-test to test whether the model with the predictor(s) is better than the intercept only model. Use \(\alpha = 0.05\).
Record whether you reject \(H_0\) at \(\alpha = 0.05\)

Please try to code up this simulation on your own! Stuck on getting the p-value out? See if you can look at the help files or use Google to figure it out; learning how to code like this is a great exercise! If you are still feeling stuck, here is a blog post I wrote a few years ago that goes through this same exercise – I strongly encourage you to try it yourself first, though!

a) Implement Strategy 1. What proportion of the time do you:

Find the quadratic term “significant” (even though the true model is linear)?
Reject \(H_0\) of your F-test (this is called a Type I error)?

b) The nominal Type I error rate should be 5%. Is the Type I error rate from Strategy 1 close to 5%? If not, explain why testing for significance and then using that result to build your model inflates Type I error.

c) Implement Strategy 2. What is the Type I error rate when testing \(H_0\) using this approach?

d) Reflect on the trade-off: Strategy 1 (test-then-add) inflates Type I error when the relationship is truly linear, but Strategy 2 (always use splines) might be overly complex. What would you recommend to a researcher in practice?

Problem 7

A cubic spline with knots at \(k_1, k_2, \ldots, k_K\) can be represented using the following basis functions:

\(b_1(x) = 1\) (intercept)
\(b_2(x) = x\)
\(b_3(x) = x^2\)
\(b_4(x) = x^3\)
\(b_{4+j}(x) = (x - k_j)^3_+ = \begin{cases} (x - k_j)^3 & \text{if } x > k_j \\ 0 & \text{otherwise} \end{cases}\) for \(j = 1, \ldots, K\)

The fitted model is: \(\hat{y} = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \sum_{j=1}^K \beta_{4+j}(x - k_j)^3_+\)

Consider the following dataset with n = 10:

x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y <- c(2, 2.5, 4, 5, 5.5, 6, 8, 11, 15, 20)

We’ll place knots at \(k_1 = 3\) and \(k_2 = 6\).

a) Manually construct the design matrix \(\mathbf{X}\) for the cubic spline with these two knots. Your matrix should have 10 rows (one per observation) and 6 columns: \([1, x, x^2, x^3, (x-3)^3_+, (x-6)^3_+]\).

Show the full \(10 \times 6\) matrix with all values filled in.

b) Fit the cubic spline model in R using your basis functions.

Report the estimated coefficients \(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2, \hat{\beta}_3, \hat{\beta}_4, \hat{\beta}_5\).

c) Using your fitted coefficients, manually calculate \(\hat{y}\) for three observations: - At \(x = 2\) (before the first knot) - At \(x = 4.5\) (between the two knots)
- At \(x = 8\) (after the second knot)

For each, show which basis functions are “active” (non-zero) and write out the prediction formula explicitly.

d) Compare your manual cubic spline to R’s built-in bs() function:

library(splines)
model_bs <- lm(y ~ bs(x, knots = c(3, 6), degree = 3))

Plot both fitted curves along with the data. They should be identical (or very close). Create a single plot showing:

The original data points
Your manual cubic spline fit
The bs() function fit (use a different color/style)
Vertical lines at the knot locations

e) Natural cubic splines add the additional constraint that the function is linear (not cubic) beyond the boundary knots. This means setting the second derivative to zero at the boundaries.

Fit a natural cubic spline using ns(x, knots = c(3, 6)) and compare it to your unconstrained cubic spline. Plot both curves. Do you notice any differences?