Playground: Chapter 8

ANOVA and Regression

This page practises the ideas of Chapter 8 in four parts. Practise the chapter works through the book’s exercises, with starter code, hints, and solutions. Go further goes beyond the book, with new questions and some of R’s built-in datasets. Check your understanding asks short questions whose answers open with a click. Do it yourself gives open tasks with an empty code space and no answers: you build and check the models yourself, as you will in your own thesis.

In the exercises, replace each ______ with your own code and press Run Code. The data frame study is prepared as in the chapter: one row per student, with the background variables, the stress and support scores, and the first-semester records.

Practise the chapter

These are the exercises at the end of Chapter 8, with the same numbers.

Exercise 1: Wellbeing across faculties

Test whether first-semester wellbeing differs between faculties with a one-way ANOVA, calculate \(\eta^2\), and, if the ANOVA is significant, run Tukey’s test.

Solution

wellbeing_anova <- aov(wellbeing ~ faculty, data = study)
summary(wellbeing_anova)
ss <- summary(wellbeing_anova)[[1]][["Sum Sq"]]
ss[1] / sum(ss)

F is 2.02 with a p-value of 0.09, above 0.05, and \(\eta^2\) is 0.013: faculty accounts for about 1% of the differences in wellbeing. The ANOVA is not significant, so there is no need for Tukey’s test; post-hoc tests look for which groups differ only after the ANOVA has shown that some do. Stress differed between faculties in the chapter; wellbeing, on this evidence, does not.

Exercise 2: Wellbeing and sleep

Fit a simple regression of first-semester wellbeing on sleep hours, and interpret the slope and R-squared.

Solution

wb_simple <- lm(wellbeing ~ sleep_hours, data = study)
summary(wb_simple)

The slope is 5.3: each extra hour of sleep goes with wellbeing about 5.3 points higher on the 0-to-100 scale. R-squared is 0.21, so sleep alone accounts for about 21% of the differences in wellbeing between students, a lot for a single variable, but far from everything. The slope describes an association: it does not show that sleeping more would raise a student’s wellbeing.

Exercise 3: Adding stress and support

Add stress and support to the model from Exercise 2, and describe what happens to the slope of sleep and to R-squared.

Solution

coef(wb_simple)["sleep_hours"]
coef(wb_multiple)["sleep_hours"]

The slope of sleep falls from 5.3 to 3.9, and R-squared doubles, from 0.21 to 0.43. Part of what looked like the effect of sleep belonged to stress: stressed students also sleep less, so in the simple model sleep took some of the credit. Among students with the same stress and support, each extra hour of sleep goes with wellbeing about 3.9 points higher, and each point of stress with wellbeing about 6 points lower.

Exercise 4: A categorical predictor

Add study_mode to the model from Exercise 3. Name the reference category and explain what its coefficient means. Then make part-time the reference category and fit the model again.

Solution

study_ref <- study |> mutate(study_mode = relevel(factor(study_mode), ref = "Part-time"))
coef(lm(wellbeing ~ sleep_hours + stress + support + study_mode, data = study_ref))

The reference category is Full-time, the first in alphabetical order, so the coefficient study_modePart-time, −3.0, means that part-time students’ wellbeing is about 3 points lower than full-time students’ with the same sleep, stress, and support (p < 0.001). With part-time as the reference, the coefficient becomes study_modeFull-time, +3.0: the same comparison seen from the other side. Nothing else in the model changes except the intercept. The choice of reference category changes how the result is written, never what it says.

Exercise 5: Burnout and dropout

Fit a logistic regression of considering dropout on the burnout score alone, calculating the score from the questionnaire first. Report its odds ratio, explain what it means, and find the predicted probability of considering dropout for students with burnout scores of 2 and 4.

Solution

burnout_model <- glm(dropout ~ burnout, data = study_b, family = binomial)
exp(cbind(odds_ratio = coef(burnout_model), confint.default(burnout_model)))
predict(burnout_model, newdata = data.frame(burnout = c(2, 4)), type = "response")

The odds ratio for burnout is 3.7 (95% CI 2.5 to 5.4): each extra point of burnout multiplies the odds of considering dropout by 3.7. An odds ratio is not a probability, which is why the predictions help: a student with a low burnout score of 2 has a 3% chance of considering dropout, and one with a high score of 4 a 31% chance. On the probability scale, the same one-point step matters more at high burnout than at low burnout, which is exactly what a constant odds ratio implies.

Exercise 6: A confounder in your own field

Draw a diagram, like the chapter’s diagram of caffeine, sleep, and GPA, of a confounder that might explain an association in your own field, and describe how a regression could separate the two explanations. Write your answer first, then open the model answer.

Model answer

An example from public health: towns with more fast-food outlets have higher rates of heart disease. Explanation 1: fast food causes heart disease (outlets → disease). Explanation 2: poverty is a common cause (poverty → outlets, poverty → disease), and the outlets themselves have no effect. The diagram has three boxes with arrows from poverty to both outlets and disease, and a question-marked arrow from outlets to disease. A regression of disease rate on the number of outlets and a measure of poverty (such as average income) separates the two: if the coefficient of outlets stays clearly above zero with poverty in the model, explanation 1 survives; if it shrinks towards zero, poverty was confounding the association. The regression can only allow for confounders that were measured; an unmeasured one, such as smoking rates, could still explain what remains.

Go further

These exercises go beyond the book.

Exercise 7: F by hand

Build the F statistic for the chapter’s three small groups of stress scores step by step, and compare it with aov().

Solution

F_value <- (ss_between / 2) / (ss_within / 9)
F_value

The between-group sum of squares is 1.31 (every score counts its group’s distance from the grand mean), the within-group sum of squares is 0.55, and F = (1.31 / 2) / (0.55 / 9) = 10.7, exactly what aov() reports, with p = 0.004. Dividing by the degrees of freedom turns each sum into an average, a mean square, so that three group averages can be compared fairly with twelve individual scores.

Exercise 8: Regression as conditional means

A regression line estimates the average outcome at each value of the predictor. Group students by their sleep, rounded to the nearest hour, calculate the average wellbeing in each group, and set the averages beside the line’s predictions.

Solution

bands |>
  mutate(line = predict(wb_line, newdata = data.frame(sleep_hours = sleep_band)))

Where most students are, the line and the group averages agree within a point: at 6 hours, 58.2 against 57.9, and at 7 hours, 63.7 against 63.2 (202 students each). At the edges, where groups are tiny (7 students at 9 hours, 3 at 10), the averages bounce around the line. This is what a regression does: it replaces a list of noisy group averages with one straight line that borrows strength from all the students, which works well as long as the averages really do lie close to a straight line.

Exercise 9: Simpson’s paradox

R’s built-in UCBAdmissions data records applications to the six largest departments of the University of California, Berkeley, in 1973. Calculate the admission rate of men and women overall, then in each department.

Solution

admissions |>
  summarise(rate = sum(Freq[Admit == "Admitted"]) / sum(Freq),
            applicants = sum(Freq), .by = Gender)

Overall, 45% of men and 30% of women were admitted. Within departments, the picture reverses: women’s admission rate is higher in four of the six departments, and far higher in department A (82% against 62%). The explanation is a confounder, the department: women applied mostly to departments C to F, which admitted few applicants of either gender, while men applied mostly to A and B, which admitted most. Comparing the overall rates mixes the effect of gender with the effect of department, exactly as the chapter’s caffeine example mixed caffeine with sleep. Holding the department constant, as a regression would, removes the apparent disadvantage.

Exercise 10: A missed curve in the residuals

The chapter found that GPA rises with study hours and then levels off. Fit a straight line of GPA on study hours, and plot its residuals against study hours with a smooth trend. Then add the squared term and plot again.

Solution

curved <- lm(gpa ~ study_hours + I(study_hours^2), data = study)

For the straight line, the smooth trend of the residuals forms an arch: below zero for students who study little (the line predicts too high), above zero in the middle, and below zero again beyond about 50 hours. A residual pattern means the model has missed something systematic, here the levelling off. With the squared term, the residuals scatter evenly around zero, and R-squared rises from 0.004 to 0.036. Both values are small, since study hours explain little of GPA either way; what matters is that only the curve describes the shape correctly. The curve peaks at about 32 hours a week.

Check your understanding

Answer each question in your own words first, then click to see a model answer.

1. What does the F statistic of an ANOVA compare?

Answer

The variation between the group averages with the variation within the groups, each divided by its degrees of freedom. An F near 1 means that the averages differ no more than chance would produce given how much individuals differ; a large F means that they differ more. The same differences between averages can give a large or a small F, depending on the spread within the groups.

2. In a multiple regression, what does “holding the other predictors constant” mean?

Answer

That the coefficient describes the difference in the outcome between students who differ by one unit in that predictor but have the same values of all the other predictors in the model. It is a statistical comparison, not an experiment: nobody was held constant, and the comparison allows only for the predictors that were measured and included.

3. What is a confounder?

Answer

A variable that affects both the predictor and the outcome, and so creates, hides, or distorts an association between them. In the chapter, sleep affects both caffeine intake and GPA, which made caffeine look harmful to grades; with sleep held constant, caffeine’s effect disappeared.

4. A residual plot shows a funnel: the residuals spread more and more as the fitted values grow. What does it mean?

Answer

The spread of the outcome is not the same everywhere: predictions for high values are less precise than for low values. The coefficients are still unbiased, but their standard errors, and so the p-values and confidence intervals, are unreliable. Remedies include transforming the outcome (for example, taking its logarithm) or using robust standard errors.

5. An odds ratio of 2 does not mean that the probability doubles. Why not?

Answer

Odds are the probability of an event divided by the probability that it does not happen. Doubling the odds nearly doubles a small probability (from 5% to about 10%), but it cannot double a large one: odds of 1 (a 50% probability) doubled become odds of 2, a probability of 67%. For that reason, logistic regression results are best reported with predicted probabilities for typical cases as well as odds ratios.

Do it yourself

These tasks have no starter code and no answers. Each code space below is empty and ready to run: write your own code, as you would for a thesis. The data (students, semesters, questionnaire, scores, and study) and the dplyr package are already loaded, and R’s built-in datasets are always available.

1. Build a regression model for first-semester wellbeing with predictors of your choice, justified by a research question. Check the residual plot, and write the paragraph reporting the model for a results section, with every number calculated by code.

2. Test whether the effect of supervisor support on wellbeing differs between full-time and part-time students, using an interaction term. Centre support first, interpret the interaction coefficient, and draw an interaction plot of the two slopes.

3. Fit a logistic regression for another yes-or-no variable, such as lives_away, with two or three predictors you can justify. Report the odds ratios with confidence intervals, and the predicted probabilities for two contrasting students.

4. R’s built-in trees data records the girth, height, and timber volume of 31 black cherry trees. Build a model that predicts volume from girth and height, check whether a straight line is adequate, and explain what you would change if it is not.

Work on your own computer

Download the Chapter 8 project

The project contains the data and all four parts of this page as an R script: the exercises with blanks, the questions to check your understanding (answers in solutions.R), and the open tasks, each with space to write your code.

Download chapter08.zip, unzip it, and double-click chapter08.Rproj.
Or type this one line in RStudio’s Console:

usethis::use_course("https://polla-fattah.github.io/data2thesis_r/playground/chapter08.zip")