• Connect the Week 6 lecture to the R output from lm().
• Extract R-squared, adjusted R-squared, and the model F-statistic.
• Use confint() for coefficient confidence intervals.
• Extract specific values from the coefficient table.
• Use predict() with interval = "confidence".
• Connect given sums of squares to R-squared and the F-statistic.
• Calculate an F-test p-value and reconstruct a coefficient interval.

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

The lecture moves from fitting a regression model to interpreting model output.

• Coefficients describe the fitted regression equation.
• Confidence intervals describe uncertainty around coefficients.
• R-squared describes how much variation the model explains.
• The F-test asks whether the model explains more than an intercept-only model.
• predict() can give a confidence interval for the predicted mean.

Workshop focus: fit one model and extract the quantities the problem set asks for.

For Week 6, use this workflow:

1. fit the model with lm().
2. save summary(model) as an object.
3. extract model-fit values from the summary object.
4. use confint() for coefficient intervals.
5. use predict() for new values and confidence intervals.

Do not stop at running the model. The important skill is understanding what the extracted value means.

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.

Part A of week_06_workshop.R follows the slides. We will complete it together with the Blomkvist data.

blomkvist <- read_csv("workshop_data/analysis_blomkvist.csv")
blomkvist <- drop_na(blomkvist, rt, age, medicine)

How well do age and medicine predict hand reaction time?

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

m_rt <- lm(rt ~ age + medicine, data = blomkvist)
m_rt
summary(m_rt)

lm() fits the model. summary() shows the coefficient table and model-fit information.

The formula says: predict hand reaction time from age and medicine score.

m_rt_summary <- summary(m_rt)
m_rt_summary

Saving the summary makes it easier to extract specific values.

m_rt_summary$r.squared
# 0.3924922

Read $ as: extract this named part from the object.

Complete Part A1 of week_06_workshop.R with me.

m_rt_summary$r.squared
m_rt_summary$adj.r.squared
m_rt_summary$fstatistic

m_rt_summary$r.squared
# 0.3924922

m_rt_summary$adj.r.squared
# 0.3878723

• R-squared: proportion of outcome variation explained by the model.
• Adjusted R-squared: similar, but takes the number of predictors into account.

R-squared is useful, but higher is not automatically better if we keep adding predictors.

• TSS is the total variation in the outcome.
• ESS is what the model explains about the data.
• RSS is what is left unexplained by the model.
• R-squared is the explained part divided by the total part.

For an ordinary least-squares model with an intercept, TSS = ESS + RSS. Explained variation does not establish causation.

m_rt_summary$fstatistic
#    value    numdf    dendf
# 84.95811  2.00000 263.00000

• value is the F-statistic.
• numdf is the numerator degrees of freedom.
• dendf is the denominator degrees of freedom.

The F-test asks whether the model explains more variation than an intercept-only model.

For an overall regression F-test with n observations and K predictor coefficients, the degrees of freedom are K and n - K - 1.

f_value <- 2.8
n <- 76
K <- 3
pf(f_value, df1 = K, df2 = n - K - 1,
   lower.tail = FALSE)

pf() is the F-distribution counterpart of pt(). The p-value is the upper-tail probability under the null; do not double it.

f_parts <- m_rt_summary$fstatistic
f_parts
pf(f_parts[1], df1 = f_parts[2], df2 = f_parts[3],
   lower.tail = FALSE)

The three entries are the F-statistic, numerator df, and denominator df, in that order.

For a question with a supplied F-statistic, n, and K, use those supplied values. You do not need to fit a model first.

Sometimes a question gives RSS, TSS, sample size n, and the number of predictor coefficients K (excluding the intercept).

rss <- 5837337.647099
tss <- 9608662.329574
n <- 266
K <- 2
ess <- tss - rss
rsq <- ess / tss
rsq

These are the sums of squares for our fitted model. R-squared is about 0.3925.

If a question instead gives ESS and RSS, obtain TSS with tss <- ess + rss, then use the same ratio ess / tss.

Df1 <- K
Df2 <- n - K - 1
(ess / Df1) / (rss / Df2)
# 84.95811

1 - (1 - rsq) * (n - 1) / (n - K - 1)
# 0.3878723

Complete Part A1b with me, then compare these answers with m_rt_summary.

Use the values supplied in a calculation question. Do not substitute values from a different fitted model.

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

confint(m_rt)
#                   2.5 %     97.5 %
# (Intercept) 269.173085 377.083962
# age           4.412560   6.430718
# medicine     -2.478923  15.048602

confint() gives confidence intervals for all coefficients in the model.

If a 95% coefficient interval excludes 0, that coefficient would usually be statistically significant at alpha = .05.

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.

The normal-model likelihood provides the estimation framework. For the coefficients, maximising likelihood is equivalent to minimising the sum of squared residuals.

confint(m_rt, parm = "age")
#       2.5 %   97.5 %
# age 4.41256 6.430718

parm = "age" asks R for the confidence interval for only the age coefficient.

confint(m_rt, parm = "medicine", level = .99)
#              0.5 %   99.5 %
# medicine -5.263471 17.83315

level = .99 asks for a 99% confidence interval.

m_rt_summary$coefficients
#              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

This table contains the estimate, standard error, t-statistic, and p-value for each coefficient.

Rows come first, then columns:

m_rt_summary$coefficients["medicine", "t value"]
# 1.412065

Read this as:

• from the coefficient table,
• take the row called medicine,
• and the column called t value.

For its p-value, use the same row and the column named Pr(>|t|):

m_rt_summary$coefficients["medicine", "Pr(>|t|)"]

Continue Part A2 in week_06_workshop.R with me.

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

predict(m_rt, newdata = new_person)
# 667.2814

The prediction is the fitted mean reaction time for a person with age = 60 and medicine = 3.

The same workflow works with more predictors: add each one with + in the formula and supply every predictor column in newdata.

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.

The problem set uses this logic: create the new data, then pass it to predict().

Finish Part A3 in week_06_workshop.R with me.

Part B of week_06_workshop.R uses the same data but predicts foot reaction time.

blomkvist_foot <- read_csv("workshop_data/analysis_blomkvist.csv")
blomkvist_foot <- drop_na(blomkvist_foot, rt_foot, age, medicine)
m_foot <- lm(rt_foot ~ age + medicine, data = blomkvist_foot)

How well do age and medicine predict foot reaction time?

Use the same extraction logic as Part A.

Part B is the main in-class practice. Begin Part C if you finish; otherwise complete it after the workshop. Part C consolidates core skills and is separate from the optional stretch script.

• Which quantity describes explained variation, and which describes residual spread?
• What null hypothesis does the model F-test test?
• How does a coefficient confidence interval differ from a confidence interval for a predicted mean?

Use this worked example after the workshop to practise the same interval calculation used in Part A2.

Made-up example: one predictor, 68 observations, coefficient estimate 2.6, and a nonzero t-statistic of 3.4 for a null coefficient of 0.

beta_hat <- 2.6
t_value <- 3.4
se_beta <- beta_hat / t_value
df_beta <- 68 - 1 - 1
tau_beta <- qt(.975, df = df_beta)
beta_hat - tau_beta * se_beta
beta_hat + tau_beta * se_beta

For 95% use .975; for 99% use .995. These are cumulative probabilities, not confidence levels passed to qt().

The stretch script uses residuals() to obtain RSS from a fitted model, then checks the model-fit calculations. The core script already practises using given sums of squares.

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

Problem-set reminder: complete the Week 6 problem set after this workshop. Focus on extracting and interpreting model output, not only on running the code.

Before next week, you should be comfortable with:

• fitting a multiple regression model with lm().
• saving and using summary(model).
• extracting R-squared, adjusted R-squared, and F-statistic values.
• calculating an F-test p-value from a statistic and its two df.
• using confint() for coefficient intervals.
• reconstructing a coefficient SE and interval from an estimate and t-statistic.
• using predict() for fitted means and confidence intervals.