• How much reaction-time variation does our model account for?
• Does a small fit improvement justify an extra predictor?
• How uncertain are the coefficients?
• What mean does the model predict, and how uncertain is it?

Check the online lecture before the workshop. Complete the Week 6 problem set after this workshop.

Work from inside your RStudio Project folder.

Open your .Rproj file for this module and check that workshop_data is inside the project folder.

1. Download week_06_workshop.R from NOW.
2. Move it into your RStudio Project folder using your file explorer.
3. Open it in RStudio and wait for each exercise checkpoint.

Run one line at a time and check the console after each line.

Complete Part A1 with me. Use hand reaction time (rt) as the outcome and age and medicine as the predictors.

blomkvist <- read_csv("workshop_data/analysis_blomkvist.csv")
blomkvist <- drop_na(blomkvist, rt, age, medicine)
m_rt <- lm(rt ~ age + medicine, data = blomkvist)
coef(m_rt)
# (Intercept)         age    medicine
#   323.12852     5.42164     6.28484

Complete the lm() formula yourself in Part A1. As in Week 5, m_rt stores the fitted model and coef(m_rt) shows its estimated intercept and slopes.

The data are adapted from a reaction-time and aging study using the Nintendo Wii Balance Board (Blomkvist et al., 2017).

Continue Part A1: save a summary of the existing model in a new object, m_rt_summary. This does not fit the model again.

m_rt_summary <- summary(m_rt)
m_rt_summary
# Coefficients section of the printed summary:
#              Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 323.12852 27.4021096 11.792104 4.784958e-26
# age           5.42164  0.5124765 10.579292 5.027243e-22
# medicine      6.28484  4.4508131  1.412065 1.591128e-01

• Estimate: the same coefficients we just saw with coef(m_rt).
• Other columns: standard errors, t-statistics, and p-values for those estimates.
• Below this table: residual standard deviation, R-squared, adjusted R-squared, and the overall F-test.

Complete Part A1 of week_06_workshop.R with me.

m_rt_summary$r.squared
# 0.3924922
m_rt_summary$adj.r.squared
# 0.3878723
m_rt_summary$fstatistic

Read $ as: extract this named part from the saved summary.

• R-squared: about 39% of hand-RT variation is explained jointly by age and medicine.
• Adjusted R-squared accounts for the number of predictor coefficients.

R-squared is the proportion of the total outcome variation explained by the model. To calculate it, divide the explained part by the total:

\[ R^2 = \frac{ESS}{TSS}. \]

• ESS — explained sum of squares: the part explained by the fitted model.
• TSS — total sum of squares: all observed variation around the outcome mean.
• RSS — residual sum of squares: the part left unexplained by the model.

The total consists of the explained and unexplained parts: TSS = ESS + RSS.

The whole circle is TSS, the total hand-RT variation in our data. Age and medicine explain the blue part, ESS; the orange part, RSS, remains unexplained.

• R-squared = ESS/TSS = 0.3925: about 39.2% explained.
• 1 − R-squared = RSS/TSS: about 60.8% unexplained.

Discuss together, without R. These simpler, made-up values illustrate the same relationship as our fitted model.

1. If ESS (blue) = 600 and TSS = 1000, how much is RSS (orange)?
2. If R-squared = .60 is explained, what proportion and percentage are unexplained?

A different example uses 48 observations and one predictor. Its explained variation (ESS) is 300 and unexplained variation (RSS) is 700.

Now apply the discussion in R: complete the first calculations in Part A1b yourself.

1. Before calculating, predict which part is larger.
2. Work out TSS, then the explained and unexplained proportions.
3. Check that the two proportions add to 1 and convert them to percentages.
4. Explain why ESS/RSS gives a different quantity.

Adding predictors cannot lower raw R-squared on the same data, even when the extra predictor adds little useful information. Adjusted R-squared accounts for model complexity and can fall.

Here n is the number of observations and K is the number of predictor coefficients, excluding the intercept. With numeric age alone, K = 1; with age and medicine, K = 2.

\[ R^2_{\mathrm{adj}} = 1-(1-R^2)\frac{n-1}{n-K-1}. \]

Continue Part A1c with your R-squared = .30, n = 48, K = 1:

rsq <- .30
adj_rsq <- 1 - (1 - rsq) * 47 / 46
adj_rsq
# 0.2847826

Return to the 600/400 discussion example, separate from the R exercise. For the same 48 participants and outcome, adding medicine to age reduces RSS only from 400 to 398.

Discuss the supplied values; do not recalculate every entry.

• Why does raw R-squared rise while adjusted R-squared falls?
• Is the tiny improvement enough to prefer the larger model automatically?

R-squared describes the explained proportion; adjusted R-squared accounts for complexity. The overall F-test asks whether the predictors explain variation beyond an intercept-only model.

Under the null hypothesis, all population predictor slopes are zero. The predictors provide no improvement over predicting the outcome mean for everyone.

\[ F=\frac{ESS/K}{RSS/(n-K-1)}. \]

Think of F as asking: how much explanatory work do the predictors do relative to the variation left unexplained?

A large F provides evidence against the null. It tests the predictors together, rather than identifying which one helps.

The degrees of freedom determine the shape of the F-distribution. They depend on:

• Model complexity: the number of predictor coefficients, K.
• Information available: the number of observations, n.

n counts observations. K counts predictor coefficients, excluding the intercept.

• df1 = K: explained variation is divided by the number of predictor coefficients.
• df2 = n − K − 1: unexplained variation is divided by the remaining degrees of freedom. We allow for K slopes and one intercept.

Start Part A1d: our exercise has 48 observations and one predictor.

exercise_n <- 48
exercise_K <- 1
exercise_df1 <- exercise_K
exercise_df2 <- exercise_n - exercise_K - 1

Work out df1 and df2, then use them to calculate (ESS / df1) / (RSS / df2). Keep using ESS = 300 and RSS = 700.

Diagram only: F = 3, df1 = 3, df2 = 46. The orange upper-tail area is the p-value.

Part A1d: continue with your calculated F, df1 = 1, and df2 = 46.

f_stat <- (300 / exercise_df1) / (700 / exercise_df2)
pf(f_stat, df1 = exercise_df1, df2 = exercise_df2, lower.tail = FALSE)

The default is the left area; lower.tail = FALSE gives the right area. Do not double it.

After the break, we will look at uncertainty in coefficients and predicted means.

Spend 20 minutes on Parts A2 and A3. These apply familiar ideas: confidence intervals, coefficient tests, and predictions from a fitted model.

• A2: obtain medicine’s 95% and 99% intervals, read its t-statistic and p-value, and rebuild its 95% interval from the estimate and t-statistic.
• A3: predict mean RT for age = 60 and medicine = 3, with confidence limits.

The script explains the newer commands. In A3, adding interval = “confidence” to predict() gives limits for the predicted population mean.

Same idea as Week 4: the interval contains population coefficient values we cannot rule out at the chosen confidence level, under the model assumptions. Here it concerns medicine’s slope after accounting for age.

confint(m_rt, parm = "medicine")
#              2.5 %   97.5 %
# medicine -2.478923 15.048602
confint(m_rt, parm = "medicine", level = .99)

parm = "medicine" selects the medicine coefficient, not a medicine score. level sets the confidence level: .95 by default, or .99 for 99%.

Review Part A2: why is the 99% interval wider than the 95% interval?

m_rt_summary
# Coefficients section of the printed summary:
#              Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 323.12852 27.4021096 11.792104 4.784958e-26
# age           5.42164  0.5124765 10.579292 5.027243e-22
# medicine      6.28484  4.4508131  1.412065 1.591128e-01

In Part A2, read the medicine row in the printed summary.

• Under t value, find medicine’s t-statistic.
• Under Pr(>|t|), find its p-value.

What do these tell us about medicine’s slope after accounting for age? How does this hypothesis differ from the overall F-test?

For the null hypothesis that a population coefficient is 0:

\[ t = \frac{\widehat{\beta}}{SE(\widehat{\beta})}, \qquad df = n - K - 1. \]

• K counts fitted predictor coefficients, excluding the intercept.
• The coefficient test uses a t-distribution, not an F-distribution.
• A coefficient interval is estimate ± t-cutoff × SE, as in Week 4.

Check Part A2 with medicine’s actual values from m_rt: estimate = 6.28484, t = 1.412065, n = 266, and K = 2 (age and medicine).

beta_hat <- 6.28484
t_value <- 1.412065
se_beta <- beta_hat / t_value
df_beta <- 266 - 2 - 1
tau_beta <- qt(.975, df = df_beta)
beta_hat - tau_beta * se_beta
beta_hat + tau_beta * se_beta

Compare with confint(m_rt, parm = "medicine"): about −2.48 to 15.05, allowing for rounding. For 95% use .975; for 99% use .995.

new_person <- tibble(age = 60, medicine = 3)
new_person

This table specifies the person whose population mean we want to predict. We will get the mean and its confidence limits in one call.

In Part A3, the table must supply both predictor values.

predict(m_rt, newdata = new_person, interval = "confidence")
#        fit      lwr      upr
# 1 667.2814 646.9741 687.5887

This interval describes uncertainty around the predicted mean.

fit is the predicted mean; lwr and upr are its confidence limits. This is not an interval for an individual person’s reaction time, and it is not sigma.

Complete Part B: predict hand reaction time from age only, using the same complete Blomkvist data.

For this actual model, n = 266, K = 1, TSS = 9,608,662 (total), and RSS = 5,881,593 (unexplained). These are summaries calculated from the fitted model/data, rounded to whole numbers.

1. Calculate ESS: how much variation is explained?
2. Calculate the explained and unexplained proportions, F, and its upper-tail p-value.
3. Check your R-squared and F against summary(), allowing for rounding.

Recommended reading: Andrews (2021), Chapter 9, Normal Linear Models.

Complete the Week 6 problem set and week_06_consolidation.R after the workshop. The consolidation practises methods already introduced in class, using new values. It is not optional stretch work.

week_06_stretch.R optionally shows how to calculate RSS from model residuals and TSS from the outcome data, then reconstruct the model-fit statistics. These extra calculations are not required exam preparation.