BCB744 Biostatistics — Practical Assessment (Version 2)

Chapters 1–14 | Total: 170 marks

Author

Affiliation

A. J. Smit

University of the Western Cape

ImportantInstructions

• Answer all eight questions. The paper totals 170 marks and you have 30 hr to complete it.

• Submit a thoroughly annotated HTML report rendered from Quarto. Show the R code used and interpret the evidence in complete sentences.

• This is a practical assessment. Marks are awarded for exploratory data analysis, statistical reasoning, correct analysis, diagnostics, interpretation, and communication. No question asks for a stand-alone theory essay.

• Before each inferential analysis, state the biological question and the null and alternative hypotheses. Identify the observational or experimental unit.

• Use plots and numerical summaries to examine the data before choosing or fitting a test. Check assumptions at the level appropriate to the analysis.

• Report estimates or effect sizes, uncertainty where available, test statistics, degrees of freedom, and exact p-values.

• Unless a question explicitly requests otherwise, use a two-sided significance level of \(\alpha = 0.05\).

• Long, unedited console dumps are not acceptable and will be penalised. Present only output that contributes to the analysis.

• Required YAML Structure: Your Quarto document must begin with the following YAML header (fill in your own details where indicated):

  ---
  title: "BCB744 Biostatistics Assessment"
  author: "Your Name"
  date: "2026"
  format:
    html:
      embed-resources: true
      toc: false
      number-sections: false
  ---

  Your document should follow this hierarchical structure:

  # Question 1
  
  ## Preamble
  ## Introduction
  ## Methods
  ## Results
  ## Discussion
  
  # Question 2
  
  ## Preamble
  ## Introduction
  ## Methods
  ## Results
  ## Discussion
  
  [... continue for the remaining questions ...]
  
  # References

Question 1 — Mammalian Body and Brain Mass (/20)

The mammals dataset in MASS contains body mass (kg) and brain mass (g) for 62 terrestrial mammal species. Each row represents one species.

1. Inspect the structure, missingness, and numerical range of both variables. Produce an informative two-panel EDA figure showing the distribution of raw body mass and raw brain mass. Summarise what the figure and appropriate measures of centre and spread reveal. (/5)

2. Apply a base-10 logarithm to both variables. Reproduce the two distribution plots on the transformed scale and explain, using evidence from your plots and summaries, whether the transformation improves their suitability for model-based analysis. (/5)

3. Treat the 62 observed species as a finite population. Set the seed to 744, then draw 2,000 random samples without replacement at each of two sample sizes, \(n=10\) and \(n=40\). For every sample, calculate the mean log₁₀ brain mass. Plot the two sampling distributions and report their empirical means, standard errors, and central 95% ranges. Explain what the comparison demonstrates about sample size and sampling uncertainty. (/6)

4. Make a scatter plot of log₁₀ brain mass against log₁₀ body mass. Describe the form, direction, strength, and any unusual observations. State one limitation of treating the 62 species as independent replicates. Do not yet fit a regression model. (/4)

TipWorked answer — Question 1

mammals_df <- MASS::mammals |>
  tibble::rownames_to_column("species")

str(mammals_df)

'data.frame':   62 obs. of  3 variables:
 $ species: chr  "Arctic fox" "Owl monkey" "Mountain beaver" "Cow" ...
 $ body   : num  3.38 0.48 1.35 465 36.33 ...
 $ brain  : num  44.5 15.5 8.1 423 119.5 ...

colSums(is.na(mammals_df))

species    body   brain 
      0       0       0 

mammals_df |>
  summarise(
    across(
      c(body, brain),
      list(
        min = min,
        q1 = ~ quantile(.x, 0.25),
        median = median,
        mean = mean,
        q3 = ~ quantile(.x, 0.75),
        max = max,
        sd = sd,
        iqr = IQR
      )
    )
  )

  body_min body_q1 body_median body_mean body_q3 body_max body_sd body_iqr
1    0.005     0.6      3.3425    198.79 48.2025     6654 899.158  47.6025
  brain_min brain_q1 brain_median brain_mean brain_q3 brain_max brain_sd
1      0.14     4.25        17.25   283.1342      166      5712 930.2789
  brain_iqr
1    161.75

mammals_df |>
  pivot_longer(c(body, brain), names_to = "measurement", values_to = "value") |>
  mutate(
    measurement = recode(
      measurement,
      body = "Body mass (kg)",
      brain = "Brain mass (g)"
    )
  ) |>
  ggplot(aes(value)) +
  geom_histogram(bins = 18, colour = "white", fill = "#2C7FB8") +
  geom_vline(
    data = mammals_df |>
      pivot_longer(
        c(body, brain),
        names_to = "measurement",
        values_to = "value"
      ) |>
      mutate(
        measurement = recode(
          measurement,
          body = "Body mass (kg)",
          brain = "Brain mass (g)"
        )
      ) |>
      summarise(median = median(value), .by = measurement),
    aes(xintercept = median),
    linetype = 2
  ) +
  facet_wrap(~measurement, scales = "free") +
  labs(x = NULL, y = "Number of species")

Raw distributions of body mass and brain mass among 62 mammal species. The dashed lines show the medians.

There are 62 species and no missing values. Both variables are extremely right-skewed. For body mass, the mean is far larger than the median because a few very large mammals dominate the upper tail. Brain mass shows the same pattern. The medians and IQRs therefore describe a typical species more faithfully than the means and standard deviations on the raw scale.

mammals_log <- mammals_df |>
  mutate(
    log_body = log10(body),
    log_brain = log10(brain)
  )

mammals_log |>
  summarise(
    across(
      c(log_body, log_brain),
      list(mean = mean, median = median, sd = sd, iqr = IQR)
    )
  )

  log_body_mean log_body_median log_body_sd log_body_iqr log_brain_mean
1     0.5808858       0.5240363    1.356357     1.904032       1.363771
  log_brain_median log_brain_sd log_brain_iqr
1         1.236743     1.062507      1.593602

mammals_log |>
  pivot_longer(
    c(log_body, log_brain),
    names_to = "measurement",
    values_to = "log_value"
  ) |>
  mutate(
    measurement = recode(
      measurement,
      log_body = "log10 body mass",
      log_brain = "log10 brain mass"
    )
  ) |>
  ggplot(aes(log_value)) +
  geom_histogram(bins = 14, colour = "white", fill = "#41AB5D") +
  facet_wrap(~measurement, scales = "free") +
  labs(x = "Transformed value", y = "Number of species")

The log transformation compresses the extreme upper tails and makes both distributions much more nearly symmetric. It also brings each mean closer to its median. This does not itself prove that a regression will be valid, but it makes a roughly linear relationship and stable residual variance more plausible. Those properties must be checked after fitting the model.

set.seed(744)

sampling_means <- bind_rows(
  tibble(
    sample_size = 10,
    sample_mean = replicate(
      2000,
      mean(sample(mammals_log$log_brain, size = 10, replace = FALSE))
    )
  ),
  tibble(
    sample_size = 40,
    sample_mean = replicate(
      2000,
      mean(sample(mammals_log$log_brain, size = 40, replace = FALSE))
    )
  )
) |>
  mutate(sample_size = factor(sample_size, levels = c(10, 40)))

sampling_summary <- sampling_means |>
  summarise(
    empirical_mean = mean(sample_mean),
    empirical_se = sd(sample_mean),
    lower_2.5 = quantile(sample_mean, 0.025),
    upper_97.5 = quantile(sample_mean, 0.975),
    .by = sample_size
  )

sampling_summary

# A tibble: 2 × 5
  sample_size empirical_mean empirical_se lower_2.5 upper_97.5
  <fct>                <dbl>        <dbl>     <dbl>      <dbl>
1 10                    1.37       0.303      0.780       1.97
2 40                    1.37       0.0982     1.18        1.56

mean(mammals_log$log_brain)

[1] 1.363771

ggplot(sampling_means, aes(sample_mean)) +
  geom_histogram(bins = 35, colour = "white", fill = "#756BB1") +
  geom_vline(
    xintercept = mean(mammals_log$log_brain),
    linewidth = 0.8
  ) +
  facet_wrap(
    ~sample_size,
    labeller = labeller(sample_size = \(x) paste0("n = ", x))
  ) +
  labs(
    x = expression("Sample mean " * log[10] * " brain mass"),
    y = "Number of simulated samples"
  )

Both sampling distributions are centred close to the finite-population mean. The distribution for \(n=40\) has a smaller empirical standard error and a narrower central 95% range than the distribution for \(n=10\). Larger samples therefore give more precise estimates of the population mean because sample-to-sample variation is reduced. Sampling without replacement also introduces a finite-population correction, which is substantial when 40 of the 62 species are selected. This simulation describes repeated samples from these 62 species; it does not by itself justify generalisation to every mammal species.

q1_provisional <- lm(log_brain ~ log_body, data = mammals_log)

mammals_log |>
  mutate(abs_residual = abs(residuals(q1_provisional))) |>
  ggplot(aes(log_body, log_brain)) +
  geom_point(colour = "#2B8CBE", alpha = 0.8) +
  geom_smooth(method = "lm", se = TRUE, colour = "#CB181D") +
  geom_text(
    data = \(x) slice_max(x, abs_residual, n = 4),
    aes(label = species),
    nudge_y = 0.12,
    check_overlap = TRUE,
    size = 3
  ) +
  labs(
    x = expression(log[10] * " body mass (kg)"),
    y = expression(log[10] * " brain mass (g)")
  )

Relationship between log-transformed body and brain mass among mammal species. Labels identify the four observations with the largest absolute residuals from a provisional straight line.

The transformed variables show a strong, positive, approximately linear association, although several species depart noticeably from the overall pattern. Species are not necessarily independent because shared evolutionary history can make closely related species more similar than expected under independent sampling. A phylogenetically informed analysis would be needed for strong comparative biological inference.

Question 2 — A Thermal-Acclimation Experiment (/16)

The following code creates a new dataset from an experiment on 24 individually tagged intertidal snails. Oxygen consumption was measured for each snail before and after a 14-day warm-acclimation treatment. Rates are in µmol O₂ h^-1.

snails <- tibble(
  snail_id = sprintf("S%02d", 1:24),
  before = c(
    5.8,
    6.1,
    4.9,
    5.4,
    6.7,
    5.2,
    5.9,
    6.4,
    5.1,
    5.6,
    6.3,
    5.0,
    5.7,
    6.2,
    4.8,
    5.5,
    6.0,
    5.3,
    6.5,
    5.4,
    5.9,
    6.1,
    5.2,
    5.8
  ),
  after = c(
    6.4,
    6.8,
    5.1,
    6.0,
    7.5,
    5.7,
    6.6,
    7.0,
    5.6,
    6.1,
    7.1,
    5.4,
    6.3,
    6.9,
    5.2,
    6.2,
    6.7,
    5.9,
    7.2,
    5.8,
    6.5,
    6.8,
    5.8,
    6.4
  )
)

1. Reshape the data into a tidy long form and produce a graphic that preserves the pairing between measurements. Add a concise numerical summary of the before, after, and within-snail change values. (/5)

2. State the experimental unit and the hypotheses. Check the distribution of the within-snail differences, explain why this is the relevant diagnostic, and conduct an appropriate inferential test. (/6)

3. Report the estimated mean change, its 95% confidence interval, the test statistic, degrees of freedom, and p-value. Interpret both the statistical and biological result. (/5)

TipWorked answer — Question 2

snails_long <- snails |>
  pivot_longer(c(before, after), names_to = "occasion", values_to = "oxygen") |>
  mutate(occasion = factor(occasion, levels = c("before", "after")))

snails_change <- snails |>
  mutate(change = after - before)

snails_change |>
  summarise(
    before_mean = mean(before),
    before_sd = sd(before),
    after_mean = mean(after),
    after_sd = sd(after),
    mean_change = mean(change),
    sd_change = sd(change),
    median_change = median(change),
    iqr_change = IQR(change)
  )

# A tibble: 1 × 8
  before_mean before_sd after_mean after_sd mean_change sd_change median_change
        <dbl>     <dbl>      <dbl>    <dbl>       <dbl>     <dbl>         <dbl>
1         5.7     0.529       6.29    0.648       0.592     0.141           0.6
# ℹ 1 more variable: iqr_change <dbl>

ggplot(snails_long, aes(occasion, oxygen, group = snail_id)) +
  geom_line(colour = "grey65", alpha = 0.8) +
  geom_point(aes(colour = occasion), size = 2.2, show.legend = FALSE) +
  stat_summary(
    aes(group = 1),
    fun = mean,
    geom = "line",
    colour = "black",
    linewidth = 1.1
  ) +
  stat_summary(
    aes(group = 1),
    fun = mean,
    geom = "point",
    colour = "black",
    size = 3
  ) +
  labs(
    x = NULL,
    y = expression("Oxygen consumption (" * mu * "mol O"[2] * " " * h^-1 * ")")
  )

The snail is the experimental unit because the treatment comparison is made within each tagged individual. Let \(d_i = \text{after}_i-\text{before}_i\). The hypotheses are \(H_0:\mu_d=0\) and \(H_A:\mu_d\ne0\).

ggplot(snails_change, aes(sample = change)) +
  stat_qq() +
  stat_qq_line() +
  labs(x = "Theoretical quantiles", y = "Observed change quantiles")

shapiro.test(snails_change$change)

    Shapiro-Wilk normality test

data:  snails_change$change
W = 0.9031, p-value = 0.02506

The paired t-test assumes that the within-pair differences are approximately normally distributed. It does not require the separate before and after samples to be normally distributed. The Q-Q plot and Shapiro-Wilk result show no important contradiction of this assumption.

q2_test <- t.test(snails$after, snails$before, paired = TRUE)
q2_test

    Paired t-test

data:  snails$after and snails$before
t = 20.533, df = 23, p-value = 2.723e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.5320579 0.6512754
sample estimates:
mean difference 
      0.5916667 

[1] 0.5916667

[1] 0.5320579 0.6512754
attr(,"conf.level")
[1] 0.95

The fitted estimate is the mean after minus before change. Its confidence interval excludes zero, and the paired test provides strong evidence against \(H_0\). Warm acclimation increased oxygen consumption in this sample by about the reported mean change in µmol O₂ h^-1. The result estimates an average response over 14 days; it does not establish whether the response persists longer or generalises beyond the sampled population.

Question 3 — Bee Repellency of Orchard Sprays (/20)

The built-in OrchardSprays dataset records the decrease in the volume of a sugar solution after honeybees were exposed to eight orchard-spray treatments. The decrease is used as an index of bee repellency. There are eight observations per treatment. For this assessment, treat each orchard position as an independent experimental unit and treatment as the sole explanatory factor; do not model rowpos or colpos.

1. Inspect the data and produce a figure that displays all observations as well as a suitable distributional summary for each treatment. Calculate a suite of suitable summary statistics and report them neatly in a table using the gt package. (/6)

2. State the experimental unit and hypotheses for an overall treatment comparison. Select an appropriate inferential approach and justify it using the study design, the distribution of model errors, the variability among treatments, and the influence of unusual observations. Show the diagnostic evidence supporting your decision. (/7)

3. Conduct the analysis selected in (b). If the overall result justifies further investigation, compare treatment pairs using a procedure that controls the family-wise error rate. Report the overall result and the adjusted pairwise evidence, then explain which treatments differ and in which direction without claiming more than the analysis supports. (/7)

TipWorked answer — Question 3

sprays <- as_tibble(OrchardSprays) |>
  mutate(treatment = factor(treatment))

spray_summary <- sprays |>
  summarise(
    n = n(),
    mean = mean(decrease),
    median = median(decrease),
    sd = sd(decrease),
    iqr = IQR(decrease),
    .by = treatment
  ) |>
  arrange(treatment)

spray_summary |>
  gt::gt(rowname_col = "treatment") |>
  gt::tab_header(
    title = "Bee repellency by orchard-spray treatment",
    subtitle = "Decrease in sugar-solution volume"
  ) |>
  gt::cols_label(
    n = "n",
    mean = "Mean",
    median = "Median",
    sd = "SD",
    iqr = "IQR"
  ) |>
  gt::fmt_number(
    columns = c(mean, median, sd, iqr),
    decimals = 2
  )

Bee repellency by orchard-spray treatment
Decrease in sugar-solution volume
	n	Mean	Median	SD	IQR
A	8	4.62	4.00	3.20	2.25
B	8	7.62	7.50	3.29	3.00
C	8	25.25	16.50	24.43	7.00
D	8	35.00	32.00	13.44	16.25
E	8	63.12	53.00	26.91	23.50
F	8	69.00	70.00	29.19	34.00
G	8	68.50	72.00	20.14	9.50
H	8	90.25	81.00	24.22	21.25

ggplot(sprays, aes(treatment, decrease, colour = treatment)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.18, show.legend = FALSE) +
  geom_jitter(width = 0.08, height = 0, size = 2.2, show.legend = FALSE) +
  labs(x = "Spray treatment", y = "Decrease in sugar-solution volume")

The experimental unit is the individual spray observation represented by a row, assuming that the eight values within each treatment arose from independent replicates. The hypotheses are \(H_0:\mu_A=\mu_B=\cdots=\mu_H\) and \(H_A:\) at least one treatment mean differs.

q3_lm <- lm(decrease ~ treatment, data = sprays)

par(mfrow = c(2, 2))
plot(q3_lm)

par(mfrow = c(1, 1))

shapiro.test(residuals(q3_lm))

    Shapiro-Wilk normality test

data:  residuals(q3_lm)
W = 0.92159, p-value = 0.0005796

bartlett.test(decrease ~ treatment, data = sprays)

    Bartlett test of homogeneity of variances

data:  decrease by treatment
Bartlett's K-squared = 42.031, df = 7, p-value = 5.128e-07

sort(cooks.distance(q3_lm), decreasing = TRUE)[1:5]

         9         27         29         58         30 
0.16736129 0.12550127 0.11642103 0.09818933 0.09601947 

4 / nrow(sprays)

[1] 0.0625

The residual plots should be read together rather than turning each diagnostic into an automatic pass/fail rule. Here the distributions differ in spread and include unusual observations. Bartlett’s test is sensitive to non-normality, but its result and the residual-versus-fitted plot both warn against assuming equal variance. I therefore use Welch’s one-way ANOVA, which relaxes the equal-variance assumption while retaining a comparison of group means.

oneway.test(decrease ~ treatment, data = sprays, var.equal = FALSE)

    One-way analysis of means (not assuming equal variances)

data:  decrease and treatment
F = 33.515, num df = 7.000, denom df = 22.843, p-value = 1.476e-10

The Welch test rejects equality of all treatment means. A Games-Howell procedure would be a natural unequal-variance post-hoc method, but it requires an additional package. A transparent base-R alternative is pairwise Welch t-tests with Holm correction:

pairwise.t.test(
  x = sprays$decrease,
  g = sprays$treatment,
  p.adjust.method = "holm",
  pool.sd = FALSE
)

    Pairwise comparisons using t tests with non-pooled SD 

data:  sprays$decrease and sprays$treatment 

  A       B       C       D       E       F       G      
B 0.57461 -       -       -       -       -       -      
C 0.48607 0.57461 -       -       -       -       -      
D 0.00623 0.01027 1.00000 -       -       -       -      
E 0.00875 0.01027 0.13889 0.26397 -       -       -      
F 0.00850 0.01027 0.08408 0.16494 1.00000 -       -      
G 0.00091 0.00120 0.02928 0.02998 1.00000 1.00000 -      
H 0.00050 0.00061 0.00249 0.00354 0.48607 0.68070 0.57461

P value adjustment method: holm 

Only pairs with Holm-adjusted p < 0.05 should be described as supported differences. Treatments D, E, F, G, and H each had a larger mean decrease than A and B. Treatment G also exceeded C and D, while H exceeded C and D. No other pair had an adjusted p-value below 0.05. The omnibus result does not mean every treatment differs from every other treatment, and non-significant pairs should not be declared equivalent.

Question 4 — Allometric Scaling of Mammalian Brains (/22)

Return to the log-transformed mammal data from Question 1. Treat this as a model-building question rather than a repeat of the EDA.

1. Quantify the association between log₁₀ body mass and log₁₀ brain mass. Select and justify an appropriate inferential measure from the form and distribution of the data, state the hypotheses, and report its estimate, 95% confidence interval, test statistic, degrees of freedom, and p-value. (/4)

2. Fit a suitable model to predict log₁₀ brain mass from log₁₀ body mass. Present the fitted equation, interpret both coefficients and the slope in allometric terms, and report the slope estimate with its 95% confidence interval, test statistic, degrees of freedom, p-value, and model \(R^2\). (/6)

3. Assess whether the fitted model satisfies its distributional and structural assumptions. Quantify case influence, identify the most influential species, and explain whether the evidence warrants deleting any observation. (/4)

4. Obtain the mean-response 95% confidence interval and individual-species 95% prediction interval for a mammal weighing 10 kg. Interpret the difference between the two intervals on the original brain-mass scale. (/4)

5. Write one concise, publication-style results paragraph integrating the principal numerical findings from (a)–(d). Distinguish the fitted mean response from prediction for a new species and avoid causal language. (/4)

TipWorked answer — Question 4

\(H_0:\rho=0\) and \(H_A:\rho\ne0\), where \(\rho\) is the population Pearson correlation between the two transformed measurements.

q4_cor <- cor.test(
  mammals_log$log_body,
  mammals_log$log_brain,
  method = "pearson"
)
q4_cor

    Pearson's product-moment correlation

data:  mammals_log$log_body and mammals_log$log_brain
t = 26.409, df = 60, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9335522 0.9755352
sample estimates:
      cor 
0.9595748 

The scatter plot from Question 1 supports a linear measure of association on the log scales. The estimated correlation is positive and strong, with its confidence interval well above zero.

q4_lm <- lm(log_brain ~ log_body, data = mammals_log)
summary(q4_lm)

Call:
lm(formula = log_brain ~ log_body, data = mammals_log)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.74503 -0.21380 -0.02676  0.18934  0.84613 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.92713    0.04171   22.23   <2e-16 ***
log_body     0.75169    0.02846   26.41   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3015 on 60 degrees of freedom
Multiple R-squared:  0.9208,    Adjusted R-squared:  0.9195 
F-statistic: 697.4 on 1 and 60 DF,  p-value: < 2.2e-16

confint(q4_lm)

                2.5 %    97.5 %
(Intercept) 0.8436923 1.0105616
log_body    0.6947503 0.8086215

The fitted equation is

\[ \widehat{\log_{10}(\text{brain})} = \hat{\beta}_0 + \hat{\beta}_1\log_{10}(\text{body}). \]

The slope is an allometric exponent: a ten-fold increase in body mass is associated with an expected increase of \(\hat{\beta}_1\) log₁₀ units in brain mass, or a multiplication of expected brain mass by \(10^{\hat{\beta}_1}\). The slope test evaluates \(H_0:\beta_1=0\) against \(H_A:\beta_1\ne0\).

par(mfrow = c(2, 2))
plot(q4_lm)

par(mfrow = c(1, 1))

q4_influence <- mammals_log |>
  mutate(cooks_d = cooks.distance(q4_lm)) |>
  arrange(desc(cooks_d))

q4_influence |>
  dplyr::select(species, body, brain, cooks_d) |>
  slice_head(n = 6)

          species    body   brain    cooks_d
1           Human  62.000 1320.00 0.12202210
2      Musk shrew   0.048    0.33 0.05105702
3   Water opossum   3.500    3.90 0.05090140
4   Rhesus monkey   6.800  179.00 0.04651557
5 Ground squirrel   0.101    4.00 0.04079438
6             Pig 192.000  180.00 0.03787819

4 / nrow(mammals_log)

[1] 0.06451613

shapiro.test(residuals(q4_lm))

    Shapiro-Wilk normality test

data:  residuals(q4_lm)
W = 0.98268, p-value = 0.5293

The residual-versus-fitted and scale-location plots assess functional form and variance; the Q-Q plot assesses residual normality; and Cook’s distance measures how strongly each row affects the fitted model. A value above the common \(4/n\) screening threshold deserves investigation, but it is not an automatic deletion rule. A biologically valid extreme species should remain unless there is evidence of measurement or transcription error. A sensitivity fit with and without an influential point may be reported.

q4_new <- data.frame(log_body = log10(10))

q4_conf <- predict(q4_lm, newdata = q4_new, interval = "confidence")
q4_pred <- predict(q4_lm, newdata = q4_new, interval = "prediction")

10^q4_conf

       fit      lwr      upr
1 47.73236 39.68098 57.41738

10^q4_pred

       fit      lwr      upr
1 47.73236 11.75855 193.7636

Exponentiating converts the fitted values and interval limits back to grams. The confidence interval estimates uncertainty in the mean brain mass for the modelled population of 10 kg mammals. The prediction interval is wider because it also includes residual among-species variation and is intended for a single new 10 kg species.

A suitable publication-style result is:

Across the 62 mammal species, log₁₀ brain mass was strongly and positively associated with log₁₀ body mass (Pearson’s \(r\) = 0.96, 95% CI [0.934, 0.976], \(t_{60}\) = 26.41, \(p\) < 0.001). The fitted relationship was $ = 0.927 + 0.752 \(\log_{10}(\text{body mass})\), with body mass explaining 92.1% of the variation in log₁₀ brain mass. The allometric slope was 0.752 (95% CI [0.695, 0.809]). For a 10 kg mammal, the estimated mean brain mass was 47.7 g (95% CI [39.7, 57.4]), whereas the 95% interval for a single new 10 kg species was [11.8, 193.8] g. The wider latter interval reflects residual among-species variation, and the observational analysis supports association rather than causation.

Question 5 — Curvature in an Enzyme Assay (/18)

The built-in DNase dataset contains optical-density measurements from an assay of deoxyribonuclease. Use only observations from Run == 1. Treat conc as the predictor and density as the response.

1. Inspect and plot the relationship. Fit two candidate models for the expected response: the first assumes a constant change in optical density per unit increase in concentration; the second adds the square of concentration as an additional predictor, entered in the original concentration units. Overlay both fitted relationships on the observations and write their equations. (/6)

2. Formally determine whether the additional term improves fit sufficiently to retain the more complex candidate. State the hypotheses in terms of that added coefficient and report the comparison statistic, degrees of freedom, and p-value. (/5)

3. Diagnose the selected model’s assumptions and case influence. Interpret the fitted relationship only within the observed concentration range, and use fitted values to quantify the expected optical density at the minimum, median, and maximum observed concentrations. Explain what these estimates reveal about how the response changes across the concentration gradient. (/7)

TipWorked answer — Question 5

dnase_1 <- as_tibble(DNase) |>
  filter(Run == 1)

str(dnase_1)

tibble [16 × 3] (S3: tbl_df/tbl/data.frame)
 $ Run    : Ord.factor w/ 11 levels "10"<"11"<"9"<..: 4 4 4 4 4 4 4 4 4 4 ...
 $ conc   : num [1:16] 0.0488 0.0488 0.1953 0.1953 0.3906 ...
 $ density: num [1:16] 0.017 0.018 0.121 0.124 0.206 0.215 0.377 0.374 0.614 0.609 ...

range(dnase_1$conc)

[1]  0.04882812 12.50000000

q5_linear <- lm(density ~ conc, data = dnase_1)
q5_quadratic <- lm(density ~ conc + I(conc^2), data = dnase_1)

q5_grid <- tibble(
  conc = seq(min(dnase_1$conc), max(dnase_1$conc), length.out = 200)
) |>
  mutate(
    linear = predict(q5_linear, newdata = pick(everything())),
    quadratic = predict(q5_quadratic, newdata = pick(everything()))
  ) |>
  pivot_longer(c(linear, quadratic), names_to = "model", values_to = "density")

ggplot(dnase_1, aes(conc, density)) +
  geom_point(size = 2.5) +
  geom_line(
    data = q5_grid,
    aes(colour = model),
    linewidth = 1
  ) +
  labs(
    x = "DNase concentration",
    y = "Optical density",
    colour = "Model"
  )

The two fitted candidate equations are

\[ \widehat{\text{density}} = 0.2597 + 0.1344(\text{conc}) \]

and

\[ \widehat{\text{density}} = 0.0942 + 0.308(\text{conc}) - 0.0144(\text{conc}^2). \]

The nested hypotheses are \(H_0:\beta_2=0\) (the straight line is adequate) and \(H_A:\beta_2\ne0\) (the quadratic term improves fit).

anova(q5_linear, q5_quadratic)

Analysis of Variance Table

Model 1: density ~ conc
Model 2: density ~ conc + I(conc^2)
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1     14 0.67824                                  
2     13 0.07680  1   0.60144 101.81 1.622e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(q5_quadratic)

Call:
lm(formula = density ~ conc + I(conc^2), data = dnase_1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.124207 -0.047433  0.006407  0.055317  0.102491 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.094221   0.029249   3.221  0.00669 ** 
conc         0.308026   0.017854  17.252 2.43e-10 ***
I(conc^2)   -0.014366   0.001424 -10.090 1.62e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07686 on 13 degrees of freedom
Multiple R-squared:  0.9858,    Adjusted R-squared:  0.9836 
F-statistic: 451.5 on 2 and 13 DF,  p-value: 9.728e-13

The nested-model F-test provides evidence that the quadratic term improves fit if its p-value is below 0.05. The comparison must be made between models fitted to exactly the same observations.

par(mfrow = c(2, 2))
plot(q5_quadratic)

par(mfrow = c(1, 1))

q5_prediction_points <- tibble(
  conc = c(
    min(dnase_1$conc),
    median(dnase_1$conc),
    max(dnase_1$conc)
  ),
  position = c("Minimum", "Median", "Maximum")
) |>
  mutate(
    fitted_density = predict(q5_quadratic, newdata = pick(conc))
  ) |>
  dplyr::select(position, conc, fitted_density)

q5_prediction_points

# A tibble: 3 × 3
  position    conc fitted_density
  <chr>      <dbl>          <dbl>
1 Minimum   0.0488          0.109
2 Median    1.17            0.435
3 Maximum  12.5             1.70 

The selected model describes a curved increase in optical density over the measured concentration range. The fitted values at the minimum, median, and maximum concentrations quantify the response on the scale of the assay. Their successive changes show that optical density increases rapidly over part of the gradient but that the increase is not constant across the full range. These predictions are interpolations supported by the observed concentrations; the fitted equation should not be used to infer behaviour beyond that range. The diagnostic plots determine whether the additional term has adequately removed systematic curvature without introducing serious variance, normality, or influence problems.

Question 6 — Sex, Body Mass, and Heart Mass in Cats (/24)

The cats dataset in MASS contains sex (Sex), body mass (Bwt, kg), and heart mass (Hwt, g) for 144 adult cats.

1. Inspect the data and use a scatter plot, group-specific fitted lines, and concise summaries to explore how heart mass varies with body mass and sex. (/5)

2. Fit a model that estimates the partial association of body mass and sex with heart mass, initially assuming that the body-mass association has the same slope in both sexes. Write the fitted equation, state term-specific hypotheses, identify the reference sex, interpret every coefficient with units, and report coefficient uncertainty and overall model fit. (/8)

3. Compare three candidate models: (i) body mass as the sole predictor; (ii) body mass and sex with a common body-mass slope; and (iii) body mass and sex with the body-mass slope allowed to differ between sexes. Use appropriate formal comparisons to determine what additional structure the data support. Explain why insufficient evidence for different slopes is not proof that the population slopes are identical. (/5)

4. For the selected model, assess distributional and structural assumptions, quantify whether the predictors contain redundant information, and identify influential observations. Conclude with a concise results paragraph containing the principal estimates and inferential evidence while distinguishing association from causation. (/6)

TipWorked answer — Question 6

cats_df <- as_tibble(MASS::cats) |>
  mutate(Sex = relevel(factor(Sex), ref = "F"))

cats_df |>
  summarise(
    n = n(),
    mean_bwt = mean(Bwt),
    sd_bwt = sd(Bwt),
    mean_hwt = mean(Hwt),
    sd_hwt = sd(Hwt),
    .by = Sex
  )

# A tibble: 2 × 6
  Sex       n mean_bwt sd_bwt mean_hwt sd_hwt
  <fct> <int>    <dbl>  <dbl>    <dbl>  <dbl>
1 F        47     2.36  0.274     9.20   1.36
2 M        97     2.9   0.467    11.3    2.54

ggplot(cats_df, aes(Bwt, Hwt, colour = Sex)) +
  geom_point(alpha = 0.75) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "Body mass (kg)",
    y = "Heart mass (g)",
    colour = "Sex"
  )

q6_body <- lm(Hwt ~ Bwt, data = cats_df)
q6_additive <- lm(Hwt ~ Bwt + Sex, data = cats_df)
q6_interaction <- lm(Hwt ~ Bwt * Sex, data = cats_df)

summary(q6_additive)

Call:
lm(formula = Hwt ~ Bwt + Sex, data = cats_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5833 -0.9700 -0.0948  1.0432  5.1016 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.4149     0.7273  -0.571    0.569    
Bwt           4.0758     0.2948  13.826   <2e-16 ***
SexM         -0.0821     0.3040  -0.270    0.788    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.457 on 141 degrees of freedom
Multiple R-squared:  0.6468,    Adjusted R-squared:  0.6418 
F-statistic: 129.1 on 2 and 141 DF,  p-value: < 2.2e-16

confint(q6_additive)

                 2.5 %   97.5 %
(Intercept) -1.8528230 1.022918
Bwt          3.4929923 4.658546
SexM        -0.6831776 0.518984

Because female (F) is explicitly the reference category:

• The intercept estimates female heart mass at Bwt = 0 kg. It is required algebraically but is outside the data and has no direct biological meaning.
• The Bwt coefficient is the expected change in heart mass, in grams, per 1 kg increase in body mass, holding sex constant.
• The SexM coefficient is the adjusted mean difference in heart mass between males and females of the same body mass.

The term-specific null hypotheses are \(H_0:\beta_1=0\) and \(H_0:\beta_2=0\), each against a two-sided non-zero alternative. Coefficient confidence intervals convey the precision of the estimates, while adjusted \(R^2\) summarises the proportion of variability explained with a penalty for model size.

q6_compare_sex <- anova(q6_body, q6_additive)
q6_compare_interaction <- anova(q6_additive, q6_interaction)

q6_compare_sex

Analysis of Variance Table

Model 1: Hwt ~ Bwt
Model 2: Hwt ~ Bwt + Sex
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    142 299.53                           
2    141 299.38  1    0.1548 0.0729 0.7875

q6_compare_interaction

Analysis of Variance Table

Model 1: Hwt ~ Bwt + Sex
Model 2: Hwt ~ Bwt * Sex
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    141 299.38                              
2    140 291.05  1    8.3317 4.0077 0.04722 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(q6_interaction)

Call:
lm(formula = Hwt ~ Bwt * Sex, data = cats_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7728 -1.0118 -0.1196  0.9272  4.8646 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9813     1.8428   1.618 0.107960    
Bwt           2.6364     0.7759   3.398 0.000885 ***
SexM         -4.1654     2.0618  -2.020 0.045258 *  
Bwt:SexM      1.6763     0.8373   2.002 0.047225 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.442 on 140 degrees of freedom
Multiple R-squared:  0.6566,    Adjusted R-squared:  0.6493 
F-statistic: 89.24 on 3 and 140 DF,  p-value: < 2.2e-16

confint(q6_interaction)

                 2.5 %      97.5 %
(Intercept) -0.6620801  6.62470490
Bwt          1.1024137  4.17041438
SexM        -8.2416012 -0.08919944
Bwt:SexM     0.0208271  3.33170228

The first comparison provides insufficient evidence that adding a constant sex difference improves the body-mass-only model. The second comparison does support allowing the body-mass slope to differ by sex. The interaction model is therefore retained. In that model, the Bwt coefficient is the female slope, SexM is the estimated male–female difference at 0 kg and is not biologically interpretable in isolation, and Bwt:SexM is the amount by which the male slope differs from the female slope. The male slope is the sum of the Bwt and Bwt:SexM coefficients. A non-significant interaction would indicate insufficient evidence of different slopes, not prove exact equality; in these data, however, the interaction comparison is significant at \(\alpha=0.05\).

par(mfrow = c(2, 2))
plot(q6_interaction)

par(mfrow = c(1, 1))

# Compute a variance-inflation factor for every non-intercept model column.
q6_x <- model.matrix(q6_interaction)[, -1, drop = FALSE]
q6_vif <- vapply(
  seq_len(ncol(q6_x)),
  \(j) {
    response <- q6_x[, j]
    others <- q6_x[, -j, drop = FALSE]
    1 / (1 - summary(lm(response ~ others))$r.squared)
  },
  numeric(1)
)
names(q6_vif) <- colnames(q6_x)
q6_vif

      Bwt      SexM  Bwt:SexM 
 9.753216 64.736123 96.871572 

cats_df |>
  mutate(
    cooks_d = cooks.distance(q6_interaction),
    leverage = hatvalues(q6_interaction),
    studentised_residual = rstudent(q6_interaction)
  ) |>
  arrange(desc(cooks_d)) |>
  slice_head(n = 6)

# A tibble: 6 × 6
  Sex     Bwt   Hwt cooks_d leverage studentised_residual
  <fct> <dbl> <dbl>   <dbl>    <dbl>                <dbl>
1 M       3.9  20.5  0.186    0.0580                 3.62
2 F       3    13    0.101    0.140                  1.59
3 M       3.7  11    0.0759   0.0408                -2.73
4 M       3.5  17.2  0.0378   0.0275                 2.35
5 M       3.6  11.8  0.0280   0.0337                -1.81
6 M       2.1  10.1  0.0265   0.0408                 1.59

4 / nrow(cats_df)

[1] 0.02777778

The residual plots assess linearity, constant variance, normality, and influence. Variance-inflation factors quantify redundancy among the three model columns. Because an uncentred continuous predictor appears in an interaction, some non-essential inflation is expected; centring body mass would reduce this without changing fitted values or the interaction test. Cook’s distance flags observations for checking and sensitivity analysis, not automatic deletion.

A suitable results paragraph is:

Heart mass increased with body mass in both sexes, but the fitted association differed between females and males. Adding a constant sex difference to the body-mass-only model did not improve fit (\(F_{1,141}\) = 0.07, \(p\) = 0.788), whereas allowing the slopes to differ did (\(F_{1,140}\) = 4.01, \(p\) = 0.0472). In the retained model, the estimated female slope was 2.64 g heart mass per kg body mass and the male slope was 4.31 g kg^-1. The fitted slope difference was 1.68 g kg^-1 (95% CI [0.02, 3.33]). The model explained 65.7% of the observed variation in heart mass. These observational data support sex-specific associations between body and heart mass, not a causal effect.

Question 7 — Theophylline Concentration Profiles (/25)

The built-in Theoph dataset contains serum theophylline concentrations (conc, mg L^-1) measured repeatedly over approximately 25 hours after oral administration to 12 subjects. It also records subject body mass (Wt, kg), dose (Dose, mg kg^-1), and elapsed time (Time, h). The object has a specialised grouped-data class rather than being an ordinary tibble.

1. Convert Theoph to a plain tibble, convert Subject to an unordered factor, and arrange observations by subject and time. Audit the resulting data by reporting its dimensions, missing values, observations per subject, and whether Wt and Dose remain constant within each subject. Explain why these checks are necessary before calculating subject-level quantities. (/5)

2. Create a new dataset containing exactly one row per subject and the following variables:

• body mass and dose;
• number of concentration measurements;
• observed peak concentration;
• time at which the observed peak occurred; and
• area under the concentration-time curve from the first to last observation, calculated by joining successive observations with straight lines and summing the resulting trapezoid areas.

Present this dataset as a clearly labelled and appropriately rounded gt table. (/7)

3. Produce one figure that shows all subject-specific concentration profiles as lines and points. Superimpose a population summary obtained by dividing elapsed time into the intervals 0–0.5, >0.5–1, >1–2, >2–4, >4–8, >8–12, and >12–25 hours, then calculating the median and interquartile range of concentration within each interval. The figure must distinguish individual trajectories from the population summary and have a self-contained caption. (/6)

4. Determine whether body mass is associated with the subject-level area under the curve calculated in (b). Produce a labelled EDA figure, select and justify a suitable inferential measure, state the hypotheses, check the relevant assumptions, and report the estimate, 95% confidence interval, test statistic, degrees of freedom, and p-value. Interpret the evidence without treating a non-significant result as proof of no association. (/7)

TipWorked answer — Question 7

theoph <- as_tibble(datasets::Theoph) |>
  mutate(Subject = factor(Subject, ordered = FALSE)) |>
  arrange(Subject, Time)

dim(theoph)

[1] 132   5

colSums(is.na(theoph))

Subject      Wt    Dose    Time    conc 
      0       0       0       0       0 

q7_audit <- theoph |>
  summarise(
    observations = n(),
    distinct_weights = n_distinct(Wt),
    distinct_doses = n_distinct(Dose),
    .by = Subject
  )

q7_audit

# A tibble: 12 × 4
   Subject observations distinct_weights distinct_doses
   <fct>          <int>            <int>          <int>
 1 6                 11                1              1
 2 7                 11                1              1
 3 8                 11                1              1
 4 11                11                1              1
 5 3                 11                1              1
 6 2                 11                1              1
 7 4                 11                1              1
 8 9                 11                1              1
 9 12                11                1              1
10 10                11                1              1
11 1                 11                1              1
12 5                 11                1              1

The data contain 132 observations from 12 subjects, with 11 observations per subject and no missing values. Each subject has one recorded body mass and one dose. Confirming ordering and within-subject constancy is essential because peak timing and trapezoid areas depend on the temporal sequence, while subject-level attributes must not be accidentally averaged across inconsistent values.

theoph_subject <- theoph |>
  summarise(
    Wt = first(Wt),
    Dose = first(Dose),
    n_measurements = n(),
    peak_conc = max(conc),
    time_to_peak = Time[which.max(conc)],
    auc = sum(
      diff(Time) * (head(conc, -1) + tail(conc, -1)) / 2
    ),
    .by = Subject
  )

theoph_subject |>
  gt::gt(rowname_col = "Subject") |>
  gt::tab_header(
    title = "Subject-level theophylline summaries",
    subtitle = "Observed concentrations over the sampled 24-hour period"
  ) |>
  gt::cols_label(
    Wt = "Body mass (kg)",
    Dose = "Dose (mg kg⁻¹)",
    n_measurements = "Measurements",
    peak_conc = "Peak concentration (mg L⁻¹)",
    time_to_peak = "Time to peak (h)",
    auc = "AUC (mg h L⁻¹)"
  ) |>
  gt::fmt_number(
    columns = c(Wt, Dose, peak_conc, time_to_peak, auc),
    decimals = 2
  )

Subject-level theophylline summaries
Observed concentrations over the sampled 24-hour period
	Body mass (kg)	Dose (mg kg⁻¹)	Measurements	Peak concentration (mg L⁻¹)	Time to peak (h)	AUC (mg h L⁻¹)
6	80.00	4.00	11	6.44	1.15	73.78
7	64.60	4.95	11	7.09	3.48	90.75
8	70.50	4.53	11	7.56	2.02	88.56
11	65.00	4.92	11	8.00	0.98	80.09
3	70.50	4.53	11	8.20	1.02	99.29
2	72.40	4.40	11	8.33	1.92	91.53
4	72.70	4.40	11	8.60	1.07	106.80
9	86.40	3.10	11	9.03	0.63	86.33
12	60.50	5.30	11	9.75	3.52	119.98
10	58.20	5.50	11	10.21	3.55	138.37
1	79.60	4.02	11	10.50	1.12	148.92
5	54.60	5.86	11	11.40	1.00	121.29

For successive observations \((t_i,c_i)\) and \((t_{i+1},c_{i+1})\), the trapezoid contribution is

\[ (t_{i+1}-t_i)\frac{c_i+c_{i+1}}{2}. \]

Summing these contributions produces one observed area under the concentration-time curve for each subject.

q7_breaks <- c(0, 0.5, 1, 2, 4, 8, 12, 25)

theoph_binned <- theoph |>
  mutate(
    time_interval = cut(
      Time,
      breaks = q7_breaks,
      include.lowest = TRUE,
      right = TRUE
    )
  )

theoph_population <- theoph_binned |>
  filter(!is.na(time_interval)) |>
  summarise(
    plot_time = median(Time),
    median_conc = median(conc),
    q1_conc = quantile(conc, 0.25),
    q3_conc = quantile(conc, 0.75),
    .by = time_interval
  ) |>
  arrange(plot_time)

ggplot(theoph, aes(Time, conc, group = Subject, colour = Subject)) +
  geom_line(alpha = 0.55, linewidth = 0.6) +
  geom_point(alpha = 0.7, size = 1.5) +
  geom_ribbon(
    data = theoph_population,
    aes(
      x = plot_time,
      ymin = q1_conc,
      ymax = q3_conc,
      group = 1
    ),
    inherit.aes = FALSE,
    fill = "grey55",
    alpha = 0.25
  ) +
  geom_line(
    data = theoph_population,
    aes(x = plot_time, y = median_conc, group = 1),
    inherit.aes = FALSE,
    colour = "black",
    linewidth = 1.2
  ) +
  geom_point(
    data = theoph_population,
    aes(x = plot_time, y = median_conc),
    inherit.aes = FALSE,
    colour = "black",
    size = 2.4
  ) +
  guides(colour = guide_legend(ncol = 2)) +
  labs(
    x = "Time since administration (h)",
    y = expression("Serum theophylline concentration (mg L"^-1 * ")"),
    colour = "Subject"
  )

Serum theophylline concentration over time for 12 subjects. Thin coloured lines and points show individual profiles. The black line shows median concentration within seven elapsed-time intervals, and the grey band spans the corresponding interquartile range.

The individual trajectories show rapid uptake followed by a slower decline, with appreciable subject-to-subject differences in peak height, peak timing, and persistence. The binned median and IQR summarise the common temporal pattern without hiding individual profiles.

ggplot(theoph_subject, aes(Wt, auc, label = Subject)) +
  geom_point(size = 2.5, colour = "#2C7FB8") +
  geom_smooth(method = "lm", se = TRUE, colour = "#CB181D") +
  geom_text(nudge_y = 3, check_overlap = TRUE, size = 3) +
  labs(
    x = "Body mass (kg)",
    y = expression("Observed AUC (mg h L"^-1 * ")")
  )

Association between body mass and the observed area under the serum theophylline concentration-time curve. Labels identify subjects.

shapiro.test(theoph_subject$Wt)

    Shapiro-Wilk normality test

data:  theoph_subject$Wt
W = 0.9742, p-value = 0.9495

shapiro.test(theoph_subject$auc)

    Shapiro-Wilk normality test

data:  theoph_subject$auc
W = 0.92514, p-value = 0.3314

q7_cor <- cor.test(
  theoph_subject$Wt,
  theoph_subject$auc,
  method = "pearson"
)
q7_cor

    Pearson's product-moment correlation

data:  theoph_subject$Wt and theoph_subject$auc
t = -1.1225, df = 10, p-value = 0.2879
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.7621143  0.2962494
sample estimates:
      cor 
-0.334528 

The subject is the independent observational unit. The hypotheses are \(H_0:\rho=0\) and \(H_A:\rho\ne0\). The scatter plot is reasonably compatible with a straight-line association, neither variable shows strong evidence against normality, and no single point obviously determines the pattern, so Pearson’s correlation is defensible. The estimated association is moderate and negative, but its confidence interval spans both substantial negative and modest positive values (\(r\) = -0.335, 95% CI [-0.762, 0.296], \(t_{10}\) = -1.12, \(p\) = 0.288). The data therefore do not provide strong evidence of an association, but the small sample and wide interval mean that an association cannot be ruled out.

Question 8 — Exercise and Resting Pulse (/25)

The survey dataset in MASS contains questionnaire responses from 237 university students. Relevant variables are sex (Sex), exercise frequency (Exer), resting pulse rate (Pulse, beats min^-1), and age (Age, years). This question compares students who reported exercising frequently with those who reported no exercise; students in the intermediate "Some" category are outside the target comparison.

1. Starting from MASS::survey, retain only Sex, Exer, Pulse, and Age. Keep the "Freq" and "None" exercise categories, remove incomplete rows, and:

• recode exercise as an ordered factor labelled "None" and "Frequently";
• create age bands "16–18", "19–21", and "22+"; and
• report, in a compact audit, the original number of rows, the number excluded as "Some", the number removed for incomplete required values, and the final number analysed.

Show enough checks to demonstrate that the resulting variables and category order are correct. (/6)

2. Calculate sample size, mean, median, standard deviation, and IQR of pulse rate for each exercise-by-sex combination. Present the result as a polished gt table with units and sensible rounding. (/5)

3. Construct one figure comparing the pulse-rate distributions between exercise groups. It must show the distribution shape, a compact quartile summary, individual observations coloured by age band, separate panels for sex, and the mean with a 95% confidence interval. Use position adjustments that keep all layers readable and provide a self-contained caption. (/7)

4. Ignoring sex for this specified unadjusted comparison, determine whether mean pulse rate differs between the two exercise groups. State the observational unit and hypotheses, use graphical and formal evidence to assess the relevant assumptions, select and conduct an appropriate analysis, and report the estimated difference with its 95% confidence interval, test statistic, degrees of freedom, and p-value. Discuss one reason why this observational comparison should not be interpreted causally. (/7)

TipWorked answer — Question 8

survey_raw <- as_tibble(MASS::survey)

q8_original_n <- nrow(survey_raw)
q8_some_n <- sum(survey_raw$Exer == "Some", na.rm = TRUE)

survey_target <- survey_raw |>
  dplyr::select(Sex, Exer, Pulse, Age) |>
  filter(Exer %in% c("Freq", "None"))

q8_incomplete_n <- sum(!complete.cases(survey_target))

survey_clean <- survey_target |>
  drop_na() |>
  mutate(
    exercise = factor(
      Exer,
      levels = c("None", "Freq"),
      labels = c("None", "Frequently"),
      ordered = TRUE
    ),
    age_band = cut(
      Age,
      breaks = c(-Inf, 19, 22, Inf),
      right = FALSE,
      labels = c("16–18", "19–21", "22+")
    )
  ) |>
  dplyr::select(-Exer)

tibble(
  stage = c(
    "Original rows",
    "Excluded: Some exercise",
    "Removed: incomplete required values",
    "Final analysed rows"
  ),
  n = c(
    q8_original_n,
    q8_some_n,
    q8_incomplete_n,
    nrow(survey_clean)
  )
)

# A tibble: 4 × 2
  stage                                   n
  <chr>                               <int>
1 Original rows                         237
2 Excluded: Some exercise                98
3 Removed: incomplete required values    28
4 Final analysed rows                   111

str(survey_clean)

tibble [111 × 5] (S3: tbl_df/tbl/data.frame)
 $ Sex     : Factor w/ 2 levels "Female","Male": 2 2 2 1 1 2 1 2 2 2 ...
 $ Pulse   : int [1:111] 104 87 83 74 80 68 89 78 72 72 ...
 $ Age     : num [1:111] 17.6 16.9 18.8 35.8 28.5 ...
 $ exercise: Ord.factor w/ 2 levels "None"<"Frequently": 1 1 2 2 2 2 2 2 2 2 ...
 $ age_band: Factor w/ 3 levels "16–18","19–21",..: 1 1 1 3 3 1 2 1 1 1 ...

count(survey_clean, exercise, Sex, age_band)

# A tibble: 12 × 4
   exercise   Sex    age_band     n
   <ord>      <fct>  <fct>    <int>
 1 None       Female 16–18        3
 2 None       Female 19–21        2
 3 None       Female 22+          2
 4 None       Male   16–18        7
 5 None       Male   19–21        2
 6 None       Male   22+          1
 7 Frequently Female 16–18       25
 8 Frequently Female 19–21        9
 9 Frequently Female 22+          7
10 Frequently Male   16–18       29
11 Frequently Male   19–21       15
12 Frequently Male   22+          9

The final dataset contains one row per student and explicitly orders the exercise comparison from no exercise to frequent exercise. The audit distinguishes exclusions made by the research question from removals caused by missing required data.

survey_summary <- survey_clean |>
  summarise(
    n = n(),
    mean = mean(Pulse),
    median = median(Pulse),
    sd = sd(Pulse),
    iqr = IQR(Pulse),
    .by = c(exercise, Sex)
  ) |>
  arrange(exercise, Sex)

survey_summary |>
  gt::gt(groupname_col = "exercise", rowname_col = "Sex") |>
  gt::tab_header(
    title = "Resting pulse rate by exercise frequency and sex",
    subtitle = "Pulse rate in beats min⁻¹"
  ) |>
  gt::cols_label(
    n = "n",
    mean = "Mean",
    median = "Median",
    sd = "SD",
    iqr = "IQR"
  ) |>
  gt::fmt_number(
    columns = c(mean, median, sd, iqr),
    decimals = 1
  )

Resting pulse rate by exercise frequency and sex
Pulse rate in beats min⁻¹
	n	Mean	Median	SD	IQR
None
Female	7	71.4	70.0	11.4	9.0
Male	10	80.5	80.0	15.2	25.8
Frequently
Female	41	73.6	72.0	12.5	12.0
Male	53	70.7	70.0	9.6	12.0

ggplot(
  survey_clean,
  aes(exercise, Pulse, fill = exercise)
) +
  geom_violin(
    trim = FALSE,
    alpha = 0.22,
    colour = "grey35",
    show.legend = FALSE
  ) +
  geom_boxplot(
    width = 0.16,
    outlier.shape = NA,
    alpha = 0.45,
    show.legend = FALSE
  ) +
  geom_jitter(
    aes(colour = age_band),
    width = 0.10,
    height = 0,
    alpha = 0.75,
    size = 1.8
  ) +
  stat_summary(
    fun.data = mean_cl_normal,
    geom = "pointrange",
    colour = "black",
    linewidth = 0.7
  ) +
  facet_wrap(~Sex) +
  labs(
    x = "Exercise frequency",
    y = expression("Resting pulse rate (beats min"^-1 * ")"),
    colour = "Age band"
  )

Resting pulse-rate distributions among students reporting no exercise or frequent exercise, shown separately by sex. Violins show distribution shape, internal boxes show medians and quartiles, coloured points are individual students classified by age band, and black points and intervals show means and 95% confidence intervals.

The figure preserves the raw observations and distribution shapes while also showing quartiles and uncertainty around each mean. The unequal group sizes remain visible rather than being concealed by bars.

ggplot(survey_clean, aes(sample = Pulse)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~exercise) +
  labs(x = "Theoretical quantiles", y = "Observed pulse-rate quantiles")

survey_clean |>
  summarise(
    shapiro_w = unname(shapiro.test(Pulse)$statistic),
    p_value = shapiro.test(Pulse)$p.value,
    variance = var(Pulse),
    .by = exercise
  )

# A tibble: 2 × 4
  exercise   shapiro_w p_value variance
  <ord>          <dbl>   <dbl>    <dbl>
1 None           0.971   0.831     200.
2 Frequently     0.978   0.107     121.

q8_test <- t.test(Pulse ~ exercise, data = survey_clean, var.equal = FALSE)
q8_test

    Welch Two Sample t-test

data:  Pulse by exercise
t = 1.3306, df = 19.642, p-value = 0.1986
alternative hypothesis: true difference in means between group None and group Frequently is not equal to 0
95 percent confidence interval:
 -2.737842 12.352361
sample estimates:
      mean in group None mean in group Frequently 
                76.76471                 71.95745 

The student is the observational unit. With the factor ordered as "None" then "Frequently", the estimated contrast is mean pulse in the no-exercise group minus mean pulse in the frequent-exercise group. The hypotheses are \(H_0:\mu_{\text{None}}-\mu_{\text{Frequently}}=0\) and \(H_A:\mu_{\text{None}}-\mu_{\text{Frequently}}\ne0\). The Q-Q plots and within-group tests do not show strong departures from normality. Welch’s two-sample procedure is a cautious choice because it does not require equal population variances and accommodates the unequal sample sizes.

The estimated difference is 4.81 beats min^-1, with a 95% confidence interval from -2.74 to 12.35 beats min^-1 (\(t_{19.64}\) = 1.33, \(p\) = 0.199). The interval includes zero, so the data do not provide strong evidence of a mean difference. This does not demonstrate equal population means. Exercise was self-reported rather than experimentally assigned, and the groups may differ in sex, age, health, medication, fitness, or other factors; the result therefore cannot be interpreted as a causal effect of exercise.

Coverage and marking map

Question	Main assessed work	Chapters	Marks
1	Data structure, descriptive statistics, visualisation, distributions, transformation, sampling uncertainty, EDA	1–6	20
2	Hypotheses, assumptions, paired t-test, estimation and interpretation	5–7, 10	16
3	Group EDA, one-way ANOVA alternatives, residual diagnostics, multiplicity-controlled comparisons	5–6, 8, 10–11	20
4	Correlation, simple regression, diagnostics, confidence and prediction intervals	9, 11–12	22
5	Linear versus polynomial regression, nested comparison, diagnostics and bounded interpretation	11, 13	18
6	Multiple regression, categorical predictors, nested models, diagnostics, collinearity and influence	11, 14	24
7	Longitudinal-data wrangling, subject-level integration, layered visualisation, association	2–5, 9–10	25
8	Filtering, recoding, grouped summaries, layered faceted visualisation, two-group inference	2–3, 5–7, 10	25
Total		1–14	170