21  Exercises: A fishy counting trick

TipConditional variables

Prerequisite for these exercises was covered within the session and is therefore not repeated here. The exercises build on this material.

21.1 Learning objectives

The learning objectives for this session are:

• Udvise forståelse for usikkerheds- og sandsynlighedsbegrebet samt grundliggende begreber indenfor biostatistik
• Udvise kendskab til de grundliggende overordnede studietyper samt at skelne mellem forklarende, eksplorative og prædiktive studier
• Redegøre for forskellige typer af tilfældig og ikke tilfældig variation
• Udvise forståelse for statistiske værktøjers begrænsninger og muligheder
• Forstå statistiske problemstillinger, der er centrale for medicin med industrial specialisering og forstå, hvordan de biostatistiske værktøjer kan appliceres på disse problemstillinger

21.2 Exercises: Revisiting smoking regulations

You are working in the Medical Regulations Department, where a historical review of evidence used to inform tobacco regulation is underway.

As part of the review, your team has been given mortality data from the British Doctors Study, the landmark study initiated by Richard Doll and Austin Bradford Hill in the 1950s. The study followed British doctors over many years and recorded their smoking habits and subsequent mortality.

The dataset available to you is an aggregated version of the historical study. It contains the number of deaths and the total person-years of observation for doctors grouped by age and smoking status.

At an initial briefing, one analyst presents a simple summary of the data:

“The mortality rate appears to be higher among smokers. The difference is quite striking.”

This catches the attention of the regulatory team. However, another epidemiologist raises an important concern:

“We should be careful about interpreting that comparison. Mortality changes substantially with age, and the smoking groups don’t necessarily contribute the same amount of observation at every age.”

The team therefore decides that simply comparing the number of deaths between smokers and non-smokers is not sufficient. Doctors contribute different amounts of observation time, and age may need to be accounted for when comparing mortality rates.

You are asked to investigate the historical data more systematically.

Your task is to examine mortality rates in relation to smoking status while accounting for differences in person-time and age. The department is particularly interested in whether a statistical model can describe the observed mortality patterns and quantify the association between smoking and mortality.

The data are aggregated, so the analysis will require you to think carefully about what the response variable represents and how the amount of observation differs between groups.

21.2.1 A: Let’s start off by considering when in general to use Poisson, logistic or normal regression.

Exercise 1A: Poisson vs. Binomial vs. Gaussian

You are designing a study to analyze the number of hospital readmissions for heart failure patients within 30 days of discharge.

• 1: Should you use a Poisson, Binomial, or Gaussian model? Justify your choice.
• 2: What are the key assumptions of the distribution you chose? How would you check them?
• 3: If the outcome were the probability of readmission (yes/no), which model would you use? Why?
• 4: If the outcome were the length of stay (in days), which model would you use? Why?

Click for the solution

1. Poisson.
2. Positive real numbers. Exponential. Mean and variance is the same.
3. Logistic regression.
4. Poisson.

Exercise 1B: Reviewing a colleagues paper.

You’re assigned as an internal peer reviewer for a study titled:

Association Between Corticosteroid Therapy and Hospital-Acquired Infections in an Intensive Care Unit

Here is the paper extract:

Background

Researchers conducted a retrospective cohort study of 612 patients admitted to an intensive care unit (ICU). The outcome was the number of hospital-acquired infections each patient experienced during their ICU stay. Patients contributed different lengths of follow-up because some were discharged early while others remained hospitalized for several weeks.

The investigators collected the following variables:

• Number of hospital-acquired infections (instances from electronic journal)
• Corticosteroid treatment (Yes/No)
• Age (years)
• Diabetes (Yes/No)
• Length of ICU stay (days)

A Poisson regression model was fitted.

Statistical Methods

Poisson regression was used to estimate the association between corticosteroid therapy and infection occurrence. Age was included as the offset variable to account for differences in patient risk over time. The model included corticosteroid treatment as main treatment effect and adjusting for diabetes.

The fitted model was

\[ \text{Infections} = Poisson(λ) \\ log(λ_i​) = β0​ + β1​(Steroid) + β2​(Diabetes) + \text{offset}(Age) \]

results

Variable	Coefficient	Standard Error	p-value
Intercept	-2.31	0.28	<0.001
Corticosteroid therapy	0.47	0.11	<0.001
Diabetes	0.29	0.10	0.004

Raw coefficients can be seen from the table above.

Patients receiving corticosteroids had 0.47 times higher odds of developing hospital-acquired infections compared with patients who did not receive corticosteroids (p < 0.001).

These findings indicate that corticosteroid therapy substantially increases the odds of hospital-acquired infection after adjustment for age and diabetes.

Discussion

After adjusting for age using an offset variable, corticosteroid therapy remained an independent predictor of hospital-acquired infection. The regression coefficient of 0.47 suggests a clinically meaningful increase in infection risk.

These results support careful monitoring of patients receiving corticosteroids in the ICU.

Your task:

• Identify all errors you can — you may ask for tips.

Click for tips

There are several mistakes:

• What was age doing in the model?
• What is the table outputting and how did they make inference in the result and discussion section about that?
• Is this Poisson or Logistic regression?

Click for the solution

1. Wrong offset

• The model uses Age as the offset.
• In Poisson regression, the offset should represent the logarithm of exposure or person-time (e.g., log(length of ICU stay)), not a patient characteristic.
• Age should be included as a regular covariate if adjustment is desired.

2. Wrong interpretation

• The text interprets the results as odds (“0.47 times higher odds”).
• Standard Poisson regression estimates log-ratios, not (log) odds ratios. Saying odds ratios or incidence rate ratio (IRRs) is double wrong (see next point).

3. Coefficients not exponentiated

• The table reports coefficients (0.47 and 0.29), but the narrative treats them as effect sizes.
• The correct interpretation requires exponentiation:
  • exp(0.47) = 1.60, meaning a 60% higher infection rate.
  • exp(0.29) = 1.34, meaning a 34% higher infection rate.

21.2.2 B: Analyzing the historic “Doctors” dataset

Let’s start by loading the dataset:

doctors <- medisr::doctors

You should also read the documentation ?medisr::doctors.

Exercise 1B

Approach the dataset naïvely by simply counting up the number of deaths for smokers and non-smokers.

• How many died?

;

Click for the solution

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

doctors |>
    group_by(smoking) |>
    summarise(n_deaths = sum(deaths))

# A tibble: 2 × 2
  smoking    n_deaths
  <fct>         <dbl>
1 non-smoker      101
2 smoker          630

• non-smoker: 101
• smoker: 630

Exercise 2B

Check the balance in the groups.

• How many years in total did the researchers follow the smoking group and how many the non-smoking group?

;

Click for the solution

doctors |>
    group_by(smoking) |>
    summarise(n_years = sum(personyears))

# A tibble: 2 × 2
  smoking    n_years
  <fct>        <dbl>
1 non-smoker   39220
2 smoker      142247

• non-smoker: 39220
• smoker: 142247

Exercise 3B

Finally, check whether the age groups are composed evenly and whether they are comparable between the smoking groups.

1. Which group is nearly half of the participants between 35 to 44 years old?
2. In the other group, how many of the age categories has higher percent distribution?

Tip: For each group, compute the percent of all participants belonging to each of the five age groups and visually compare them one by one.

;

Click for the solution

doctors |>
  group_by(smoking) |>
  mutate(percent = 100 * personyears / sum(personyears)) |>
  ungroup() |>
  select(age, smoking, percent)

# A tibble: 10 × 3
     age smoking    percent
   <dbl> <fct>        <dbl>
 1     1 smoker       36.8 
 2     2 smoker       30.4 
 3     3 smoker       20.1 
 4     4 smoker        8.90
 5     5 smoker        3.74
 6     1 non-smoker   47.9 
 7     2 non-smoker   27.2 
 8     3 non-smoker   14.6 
 9     4 non-smoker    6.59
10     5 non-smoker    3.73

age smoker `non-smoker`

1 1 36.8 47.9 2 2 30.4 27.2 3 3 20.1 14.6 4 4 8.90 6.59 5 5 3.74 3.73

1. non-smokers
2. 2—5, so the 4 older groups.

Exercise 4B

Now, run the model naïvely, without accounting for these differences.

The model will be:

\[ deaths = Poisson(λ) \\ log(λ) = β0 + β1_{smoking} \] Report the (incorrect) estimand for incidence rate ratio for smoking vs. non-smoking.

;

Click for the solution

m4b <- glm(deaths ~ smoking,
                     family = poisson(link = "log"),
                     data = doctors)

m4b |> parameters::parameters(exponentiate=TRUE)

# A tibble: 2 × 9
  Parameter     Coefficient    SE    CI CI_low CI_high     z df_error
  <chr>               <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>
1 (Intercept)         20.2  2.01   0.95  16.5    24.4   30.2      Inf
2 smokingsmoker        6.24 0.669  0.95   5.08    7.74  17.1      Inf
# ℹ 1 more variable: p <dbl>

21.3 Parameter | IRR | SE | 95% CI | z | p

(Intercept) | 20.20 | 2.01 | [16.51, 24.40] | 30.21 | < .001 smoking [smoker] | 6.24 | 0.67 | [ 5.08, 7.74] | 17.08 | < .001

Exercise 5B

Let’s try to correct for these issues.

Provide the scientific regression equation we will run

Click for the solution

The model will be:

\[ deaths = Poisson(λ) \\ log(λ) = β0 + β1_{smoking} + β2_{age} + \text{offset}(personyears) \]

Note that the offset does not have a β coefficient associated with it!

Exercise 6B

The following R code contains two mistakes. Identify and fix them:

m5b <- glm(deaths ~ smoking + age,
                     family = poisson(link = "logit"),
                     offset = personyears,
                     data = doctors)

Click for the solution

# Mistake 1: the offset should be log-personyears!
# Mistake 2: the link should be log (logit is for logistic regression)!
m6b <- glm(deaths ~ smoking + age,
                     family = poisson(link = "log"),
                     offset = log(personyears),
                     data = doctors)

21.3.1 C: Writing a Clinical Report

Exercise 3A:

The final model is provided in Exercise 6B. If you didn’t make that exercise, go ahead and see the solution and copy the code. You will not make inference based on the model.

You will write a clincal report based on the inference we got from the model.

Structuring the Report

• 1: How would you explain the incidence rate ratio (IRR) to a clinician who is not familiar with statistics?
• 2: Write a 1-paragraph summary of the findings from your Poisson model for a clinical audience.
• 3: Compare these results to the naïve model from Exercise 4B
• 4: How would you visually present the effect of smoking on mortality in a figure for the report? Describe the plot.
• 5: Now implement the plot using {ggplot2}.
• 6: Make a legend to the plot so that others can read it with the intention you had.

No solution is given here, as you can be creative 😄

Please discuss your ideas with the instructors!