Practical #9

Advanced Statistical Programming using R — Modelling & Analysis

Author

Leonhard Kestel, Lisa Bondo Andersen, Cynthia Huang

Published

June 19, 2026

Quiz

Before starting, work through this QUIZ to check your understanding of this week’s lecture on analysis & inference (with a short recap of last week’s open data & reproducibility material).

Overview

In yesterday’s lecture you saw how a vague research question (“we will use time series analysis”) becomes a specific, estimable model — one that names the response variable, the predictors, and the method. Today you walk that whole path: from a loose question, to a classified analysis type, to a written model specification, to a hypothesis test in R with the {infer} package.

  • Part 1 - sandbox:
    • starwars teaching dataset
    • answers are known and we give you solutions.
  • Part 2 - the real thing:
    • Your own project
    • Open a GitHub Issue for your next analysis step

Work in your project group. At least one person should have the group project open in RStudio, with proposal.qmd and your dataset to hand.

Note

You will need the infer package. Install it inside your project (so renv records it) and snapshot afterwards:

install.packages("infer")
renv::snapshot()   # run in the R console

Part 1: Sandbox — an analysis plan for starwars

~40 min

{dplyr} ships with starwars, a dataset of 87 Star Wars characters. A few of its columns:

  • height — height in cm (<int>, some missing)
  • mass — weight in kg (<dbl>, some missing)
  • gender — "masculine", "feminine", or NA
  • species, homeworld, eye_color — categorical descriptors

Take this deliberately vague research question:

“Are masculine and feminine Star Wars characters different heights?”

In this part you sharpen it into something estimable and test it. Keep your answers short — the point is to rehearse the four moves you will repeat on your own data in Part 2.

Exercise 1.1: Rewrite the question as a specific claim

Turn the vague question above into one sentence of this form:

We will estimate outcome ~ predictor using [method], and report [specific output].

Name actual starwars columns and a method you have studied.

We will compare mean height between masculine and feminine characters (height ~ gender) using a difference in means test, and report the estimated difference, a 95% confidence interval, and a p-value.

Note how much the vague version left open: which “size” variable? which groups compared? compared how? The specific version answers all three — and a group member could now start coding from it.

Exercise 1.2: Classify the analysis

As in the lecture, classify the analysis along two independent axes:

Which goal?

  • Descriptive — what is the state of things? (a summary, a mean, a rate)
  • Explanatory — is there a real relationship, and how strong is it?
  • Predictive — how well can we predict a new, unseen case?

Exploratory or confirmatory? — was the question fixed before you looked at the data (confirmatory), or did the data suggest it (exploratory)? This is separate from the goal: any goal can be pursued either way.

Goal: this is explanatory — we are asking whether gender is related to height, not just describing one group.

Exploratory or confirmatory? We did not pre-register this question before seeing the data, so it is exploratory: the result is a candidate finding, not a confirmed one. And because starwars is observational (no one assigned characters a gender), any difference is associational, not causal. Stating these limits is a sign of good judgement, not weakness.

Exercise 1.3: Write the model specification

Write out the full specification for the claim from Ex 1.1:

  • Response variable — which column, and what type?
  • Predictor — which column, and what type?
  • Method — the test you will run
  • Expected output — the exact objects you will produce
  • Interpretation sentence — one sentence connecting the result to the question
  • Response: height — height in cm (<int>)
  • Predictor: gender — masculine / feminine, a two-level grouping (<chr> → treated as a factor)
  • Method: difference in means between the two groups, using a theory-based t null distribution
  • Expected output: the observed difference in means, a null-distribution plot with the observed statistic shaded, a p-value, and a 95% CI
  • Interpretation: “This tells us whether masculine characters are systematically taller or shorter than feminine characters in this dataset.”

Exercise 1.4: Run the test with {infer}

As in the lecture, the theory-based {infer} pipeline has two short branches that share a specify() step:

specify() → assume(distribution)                     →┐
specify() → calculate(stat) → obs_stat               →┴→ get_p_value() / get_confidence_interval() / visualize()

First prepare the data — assume("t") needs a numeric response and a two-group predictor with no missing values:

library(dplyr)
library(infer)

sw <- starwars |>
  select(gender, height) |>
  filter(gender %in% c("masculine", "feminine")) |>
  tidyr::drop_na()

sw |> count(gender)   # how many characters in each group?

Now run the difference-in-means test for height ~ gender:

# 1. The theoretical (t) null distribution
null_dist <- sw |>
  specify(height ~ gender) |>
  assume("t")

# 2. The observed difference in means (masculine − feminine)
obs_diff <- sw |>
  specify(height ~ gender) |>
  calculate(stat = "diff in means", order = c("masculine", "feminine"))

# 3. The p-value
null_dist |>
  get_p_value(obs_stat = obs_diff, direction = "two-sided")

# 4. A 95% confidence interval around the observed difference
null_dist |>
  get_confidence_interval(point_estimate = obs_diff, level = 0.95)

# 5. Visualise where the observed difference falls in the null distribution
null_dist |>
  visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "two-sided")

What does each step do, and what does the result tell you?

  • specify(height ~ gender) declares the response and the explanatory variable.
  • assume("t") builds the null distribution from the theoretical t-distribution — no simulation needed, because mathematicians already derived the shape for a difference in means.
  • calculate(stat = "diff in means", order = ...) reduces the data to one number: the masculine − feminine difference. order fixes the sign so the result is interpretable.
  • get_p_value() reports how often a t-null this extreme would occur; get_confidence_interval() gives a plausible range for the true difference.
  • visualize() + shade_p_value() plots the t-null with the observed statistic shaded.

Read the output honestly: the feminine group is small (count it in the prep step — only a handful of characters), and Star Wars heights vary enormously (Yoda at 66 cm, Chewbacca at 228 cm, droids included). With a small group and high spread, expect a wide confidence interval and weak evidence of any difference. That is a perfectly good result — see the note on null findings in Part 2.

assume() works because the t-distribution has a known formula. When your statistic has no standard distribution (a ratio, a trimmed mean, a custom score), you can instead simulate the null by shuffling group labels — a permutation test. The pipeline shape barely changes:

null_sim <- sw |>
  specify(height ~ gender) |>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  calculate(stat = "diff in means", order = c("masculine", "feminine"))

null_sim |>
  get_p_value(obs_stat = obs_diff, direction = "two-sided")

For a standard test like this one, assume("t") and generate() give essentially the same answer. Great explainer: Andrew Heiss, Null Worlds — https://nullworlds.andrewheiss.com.

Question specify() assume() calculate(stat = ...)
Two groups, continuous response → difference in means y ~ group "t" "diff in means"
Two categorical variables → association y ~ x (both factors) "Chisq" "Chisq"
One mean vs a fixed value response = y "t" "mean"
Two continuous variables → correlation y ~ x "t" "correlation"

Prefer base R? t.test(), chisq.test(), cor.test(), and lm() answer the same questions — use whichever you can explain. (The summary of an lm() is inference: it reports a t-test per coefficient.)

Part 2: Your project — refine your own analysis plan

~70 min

Now do exactly what you just did in the sandbox, but on your own dataset and as a group. A plan is only useful if a group member can pick it up and start coding without asking what you meant: “We will analyse the effect of weather on cycling” is a topic, not a plan.

Open proposal.qmd alongside yesterday’s slides. At least one person should have the dataset loaded.

Exercise 2.1: Rewrite each research question as a specific claim

For each research question you are keeping, write one sentence in the form from Ex 1.1, naming the actual columns from your dataset and a method you have studied.

Vague: “We will use regression to analyse the effect of weather on cycling.”

Specific: “We will regress gesamt (total daily count) on sonnenstunden (sunshine hours) and niederschlag (precipitation), controlling for day of week, using linear regression. We will report a coefficient table and a residual plot. This addresses RQ1 because it quantifies the relationship between weather and cycling volume while accounting for weekly patterns.”

Only the second version gives a group member enough to start working — and gives you something to defend in the oral exam.

Exercise 2.2: Classify each analysis

Using the two axes from Ex 1.2, label each research question in proposal.qmd: its goal (descriptive / explanatory / predictive) and whether it is exploratory or confirmatory.

Tip

If you find yourself writing “explanatory” for an observational dataset, add a sentence on what you cannot conclude — just as in the starwars sandbox, most project datasets support associational claims, not causal ones. And if the question came from looking at the data, label it exploratory and treat the finding as provisional.

Exercise 2.3: Write and sanity-check a model specification

Pick the single research question you are most confident about — the one you will start working on first. Write its full specification (response, predictor(s), method, expected output, interpretation), then confirm the data can actually support it:

library(dplyr)

dat <- readr::read_csv("data/raw/your-data.csv")
glimpse(dat)

# Do the variables you named exist, with the right type?
# Is there enough variation and coverage to estimate what you want?
dat |> count(your_group_variable)          # are groups large enough?
dat |> summarise(across(where(is.numeric), \(x) sum(is.na(x))))  # missingness

If a variable is missing, the wrong type, or too sparse, revise the specification now rather than discovering the problem mid-analysis.

Important

The oral exam tests reasoning, not memorisation — you will be asked why you chose this method and how it relates to your research question. If your proposal was drafted with heavy AI assistance, read this specification line by line and make sure you can defend every choice. “The model said so” is not an answer you can give — you need to own the decision.

Exercise 2.4: Run an {infer} test on your own data

Apply the workflow from Ex 1.4 to one testable claim from your project. Use the lookup table in Part 1 to pick the right specify() / assume() / calculate() combination, then run the pipeline:

library(infer)

# Replace the placeholders, following the Ex 1.4 pattern
obs_stat <- dat |>
  specify(your_outcome ~ your_predictor) |>   # or  response = your_outcome
  calculate(stat = "diff in means",           # see the lookup table for your case
            order = c("group_a", "group_b"))

null_dist <- dat |>
  specify(your_outcome ~ your_predictor) |>
  assume("t")                                 # the right distribution for your test

Write down, in one sentence, what the result means for your research question.

Tip

You do not need a “significant” result. A clearly-reported null or boring finding is just as valuable — articulating why you tested it matters more than the p-value. And because you almost certainly chose this test after exploring your data, label the finding as exploratory / provisional: it is a candidate worth reporting, not a confirmed conclusion. Note what fresh data would be needed to confirm it.

Exercise 2.5: Open a GitHub Issue for the analysis

Turn the specification from Ex 2.3 into a GitHub Issue in your project repository so a group member can pick it up:

  • Goal — what you are trying to produce (the expected output)
  • Data — the dataset and the specific variables
  • Owner — who is responsible
  • Target — a rough date for a first draft
Tip

Assign the issue to a group member before the session ends. Issues that are assigned get done; issues that are not get forgotten.

Reflection log

Take a few minutes to add this week’s entry to your reflection log. Then commit and push everything — the refined proposal.qmd, any analysis script, and the reflection log entry can all go in one commit.

End-of-practical checklist

Before you leave today, make sure your group has:

Resources

Statistical inference

{infer}