The R boxes in this course assume a working installation of two free programs. R is the statistics language that does the calculations, and RStudio is the program in which R code is written and run. Both are available from the Posit download page; R is installed first and RStudio second.

1. Create an RStudio project for the course by choosing File, then New Project, then New Directory, then New Project, and naming the folder (for example, hsci410). A project is a folder that RStudio remembers; opening its .Rproj file later reopens the same folder with the same settings, and file paths in the code are read relative to it.
2. Install the packages used in this lesson by typing the line below into the Console pane and pressing Enter. A package is a free add-on that supplies extra functions. Installation needs an internet connection, takes a few minutes, and is done once per computer.
3. Load the packages at the top of each script with library(). Loading has to be repeated in every new R session.
4. Open a new script with File, New File, R Script, and save it inside the project. Code typed into a script is saved with the project; code typed only into the Console is lost when RStudio closes.

# Run once per computer, in the Console
install.packages(c("tidyverse", "here", "mediation", "dagitty"))

# Run at the top of every script that uses them
library(tidyverse)   # data handling, graphics, reading and writing files
library(here)        # builds file paths from the project folder

Two ways of running code appear in every lesson. Running a line sends only the current line, or the highlighted lines, to the Console; in RStudio this is the Run button or Ctrl+Enter (Cmd+Return on a Mac). Sourcing a script runs the whole file from top to bottom; this is the Source button or Ctrl+Shift+S (Cmd+Shift+S on a Mac). Running line by line suits exploration, and sourcing confirms that the complete script works from a fresh start.

Any text after a # symbol on a line is a comment. R ignores comments, so they are used to explain what each step does and why. The code boxes in this course use comments in this way, and the answer-key scripts follow the same practice.

About tidyverse. The line library(tidyverse) loads several packages at once, including dplyr (data manipulation), readr (reading and writing files), tidyr (tidying data) and ggplot2 (graphics). If the tidyverse installation fails, the four packages can be installed separately with install.packages(c("dplyr", "readr", "tidyr", "ggplot2")) and loaded with four library() lines in place of library(tidyverse). Together they supply every tidyverse function used in this lesson.

The dagitty package reads a DAG written as text and reports which variables must be adjusted for to estimate a chosen effect. The same tool is available as a free web page, dagitty.net, where the diagram is drawn with the mouse; the R version keeps the diagram in the script beside the analysis. The code below encodes the diagram from the worked example, with cvd standing for the CVD event and bp for blood pressure.

library(dagitty)

# Each "A -> B" is one arrow: A causes B
g <- dagitty("dag {
  age -> smoker ; age -> cvd
  sex -> smoker ; sex -> cvd
  smoker -> bp ; bp -> cvd
  smoker -> cvd
  smoker -> hospital ; cvd -> hospital
}")

# Variables to adjust for to estimate the TOTAL effect of smoking on cvd
adjustmentSets(g, exposure = "smoker", outcome = "cvd", effect = "total")

# Variables to adjust for to estimate the DIRECT effect (not through bp)
adjustmentSets(g, exposure = "smoker", outcome = "cvd", effect = "direct")

Reading the output. The first line answers the first call (the total effect), and the second line answers the second call (the direct effect). Each is a sufficient adjustment set: if the DAG is correct, including these variables in the model removes confounding of that effect. For the total effect of smoking, the set is age and sex. Blood pressure is left out because it is a mediator, and hospital admission is left out because it is a collider. For the direct effect, blood pressure joins the set, because the direct effect is defined as the effect that remains when the mediator is held fixed. The answer is only as good as the arrows typed in: dagitty checks the logic of the diagram, and only subject-matter knowledge can check whether the diagram matches reality.

The simulation below creates a population of 10,000 people in which smoking and diabetes are unrelated by construction. Each condition raises the chance of a hospital admission, which makes admission a collider (smoking → admission ← diabetes). The code compares the percentage with diabetes among smokers and non-smokers, first in the whole population and then among hospital patients only. The function rbinom(n, 1, p) draws n values of 0 or 1 with probability p of a 1; table() counts each combination of values; prop.table(..., margin = 1) turns the counts into proportions within each row; and subset() keeps the rows that meet a condition.

set.seed(2026)
n        <- 10000
smoker   <- rbinom(n, 1, 0.30)   # 1 = smokes (30% of people)
diabetes <- rbinom(n, 1, 0.10)   # 1 = diabetes (10%), unrelated to smoking
# Each condition raises the chance of a hospital admission (the collider)
p_admit  <- 0.05 + 0.30 * smoker + 0.30 * diabetes
hospital <- rbinom(n, 1, p_admit)
pop      <- data.frame(smoker, diabetes, hospital)

# Whole population: % with diabetes among non-smokers (0) and smokers (1)
round(100 * prop.table(table(smoker = pop$smoker, diabetes = pop$diabetes),
                       margin = 1), 1)

# Hospital patients only: this restriction conditions on the collider
inpatients <- subset(pop, hospital == 1)
nrow(inpatients)
round(100 * prop.table(table(smoker = inpatients$smoker,
                             diabetes = inpatients$diabetes), margin = 1), 1)

Reading the output. Rows are smoking status (0 = does not smoke, 1 = smokes) and columns are diabetes status, so the right-hand column is the percentage with diabetes. In the whole population, about 10% of both groups have diabetes (10.0% and 9.8%), which matches the way the data were built. Among the 1,688 hospital patients, 44% of non-smokers and only 15% of smokers have diabetes, so smoking now appears to protect against diabetes. The association is produced entirely by the sampling: a non-smoker in hospital is likely to be there because of diabetes, whereas a smoker in hospital is often there because of smoking. A study that recruited only inpatients would report this false association even with flawless data entry and modelling.

The DAG: education → income → health, with education → health directly. We will (i) run Baron & Kenny’s three regressions by hand, then (ii) replicate the result with the mediation package, which gives proper bootstrap confidence intervals for the indirect effect.

# install.packages(c("mediation", "dagitty"))
library(mediation)

# 1. Simulate data that match the DAG: education -> income -> health,
#    with a smaller direct path education -> health.
set.seed(410)
n         <- 800
education <- rnorm(n)
income    <- 0.6 * education + rnorm(n)              # path a
health    <- 0.3 * education + 0.5 * income + rnorm(n)  # direct + path b
dat       <- data.frame(education, income, health)

# 2. Baron & Kenny by hand --------------------------------------------------
#    Step 1: total effect c  (health on education)
coef(lm(health ~ education, data = dat))["education"]

#    Step 2: a (income on education)
fit_M  <- lm(income ~ education, data = dat)

#    Step 3: direct c' (education) and b (income), from health on both
fit_Y  <- lm(health ~ education + income, data = dat)
coef(fit_Y)                     # c' on education, b on income

# Indirect effect = a * b  (or equivalently c - c')
a <- coef(fit_M)["education"]
b <- coef(fit_Y)["income"]
a * b

# 3. Same answer, with bootstrap CIs, via the mediation package -------------
med <- mediate(fit_M, fit_Y,
                treat    = "education",
                mediator = "income",
                boot     = TRUE, sims = 1000)
summary(med)

What each line does.

Console output, with the four paths marked. The first three printed results come from the by-hand part of the code. The highlighted numbers are c (total effect, first result), c′ and b (the education and income coefficients in fit_Y), and a × b (the indirect effect, printed under the name education because a carries that label). The (Intercept) is the predicted health score when every predictor equals zero and is not needed here.

The value of a is computed by the code without being printed. Printing the coefficients of the mediator model shows it:

coef(fit_M)    # the education coefficient is path a

Reading the summary(med) table. Each row is one effect. ACME stands for average causal mediation effect, the indirect effect through income (a × b). ADE stands for average direct effect, the effect of education that does not pass through income (c′). Total Effect is c, and Prop. Mediated is the share of the total effect that travels through income: 0.297 ÷ 0.647 = 0.459. The Estimate column repeats the by-hand values, because the bootstrap adds uncertainty intervals without changing the estimates.

The 95% confidence interval. For the ACME, the interval runs from 0.243 to 0.350. A confidence interval shows how precisely an effect has been estimated: the values inside it are those most compatible with the data, and the method that produces it captures the true value in 95% of repeated samples. In this simulation the true indirect effect, 0.30, lies inside the interval. A narrow interval indicates a precise estimate, and an interval that includes 0 indicates that the data cannot rule out the absence of an indirect effect.

The p-value. The p-value answers the question of how often an estimate at least this far from zero would arise by chance if the true effect were zero. R prints < 2.2e-16 for any p-value below 0.00000000000000022. Here the p-value comes from the 1,000 bootstrap resamples, none of which fell on the other side of zero, so the package records it as zero. With 1,000 resamples the smallest p-value distinguishable from zero is 0.002, so the result is reported as p < 0.002. A p-value this small means the data are very hard to reconcile with an effect of zero. The stars repeat the same information using the thresholds listed on the Signif. codes line. The p-value carries no information about the size or importance of an effect; the estimate and its confidence interval carry that information.

A model reporting sentence. “In simulated data (n = 800), the estimated indirect effect of education on health through income was 0.30 units of the simulated health score per standard deviation of education (95% CI 0.24 to 0.35), the direct effect was 0.35 (95% CI 0.28 to 0.43), and the total effect was 0.65 (95% CI 0.58 to 0.72); income carried an estimated 46% of the total effect (95% CI 38% to 54%).”

Two cautions. (1) The estimate of the indirect effect is only as credible as the DAG. If an unmeasured variable confounds income and health, the indirect estimate is biased even though every regression runs cleanly. (2) Both the Baron & Kenny arithmetic and mediate() as called here assume that the effect of income on health is the same at every level of education, which is called the assumption of no exposure-mediator interaction. To relax it, the outcome model is fitted as lm(health ~ education + income + education:income, data = dat), where education:income is the interaction term; mediate() then reports separate indirect effects, labelled ACME (control) and ACME (treated), for education values of 0 and 1.

The structured approach starts with structured files. The convention below, with one RStudio project, one folder per stage, and one numbered script per task, scales from a homework assignment to a journal-ready paper. The steps are as follows.

1. Open the RStudio project created in the setup box (File, Open Project), so that R works inside the project folder.
2. Run the four dir.create() lines below once, in the Console. They create the folders data/raw, data/processed, R, and output/figures inside the project.
3. Download cohort.csv from the link above and move it into the new data/raw folder. The file describes a small smoking cohort: 600 participants, 36 of whom have no recorded outcome because they were lost to follow-up.
4. Create a new script, save it as R/01_load_clean.R, paste in the code from library(tidyverse) onward, and run it line by line.
5. Compare the console with the output shown below the code.

# Create directories from R (or by hand). Run once at project start.
dir.create("data/raw",        recursive = TRUE)
dir.create("data/processed",  recursive = TRUE)
dir.create("R");  dir.create("output/figures", recursive = TRUE)

# tidyverse: dplyr (manipulation), ggplot2 (graphics), readr (file IO),
# tidyr (reshape), stringr (text). Install once.
# install.packages(c("tidyverse", "here"))
library(tidyverse)
library(here)                                # builds file paths from the project folder

# A canonical pipeline: read -> clean -> save -> analyse
raw <- read_csv(here("data/raw/cohort.csv"))
clean <- raw |>
  filter(!is.na(outcome)) |>
  mutate(age_grp = cut(age, c(0, 30, 50, 70, Inf)),
         smoker  = factor(smoker, levels = c("No", "Yes")))
write_csv(clean, here("data/processed/cohort_clean.csv"))

# Sketch a DAG to anchor the analysis (see the earlier DAG course)
# library(dagitty)
# g <- dagitty("dag { smoker -> outcome ; age -> smoker ; age -> outcome }")

Checking the data before and after cleaning. When read_csv() runs, it prints a short report of what it found. chr marks text (character) columns and dbl marks numeric columns (“double” is R’s name for a number that can have decimals). The message is informational, and adding show_col_types = FALSE to read_csv() switches it off.

Before any cleaning, glimpse() shows every variable on one line, with its type and its first few values:

glimpse(raw)    # one line per variable: name, type, first values

The output can be read against a short codebook for the file:

Counting rows before and after the filter confirms what the pipeline did, and count() shows how many participants fall in each new age band:

nrow(raw)       # rows before filtering
nrow(clean)     # rows after dropping participants with no outcome
count(clean, age_grp)   # participants in each age band

Reading the output. The raw file has 600 rows and the cleaned file 564, so 36 participants were removed. In filter(!is.na(outcome)), is.na() asks whether each outcome is missing and ! means “not”, so the filter keeps the rows whose outcome is known. In the age bands, (30,50] means older than 30 and up to and including 50; a round bracket excludes the boundary and a square bracket includes it.

Dropping missing outcomes is a decision to record. Keeping only participants with a known outcome is a complete-case decision: the analysis will describe only the people who were followed to the end. If those lost to follow-up differ from the rest (for example, if sicker participants were more likely to drop out), the remaining 564 no longer represent the cohort, and the results may be biased. The decision therefore belongs in the file log, with the number of rows removed and the reason, as in the two lines below. The date in the log will be the date on which the code is run.

# Record the complete-case decision in the file log
cat(format(Sys.Date()), "cohort_clean.csv:", nrow(raw) - nrow(clean),
    "participants with no recorded outcome excluded;", nrow(clean), "remain\n",
    file = here("data/file_log.txt"), append = TRUE)
readLines(here("data/file_log.txt"))

The pipe operator |> (or %>%) is the workhorse of the tidyverse: it chains verb -> verb -> verb so analysis reads top to bottom. Combined with here() for paths, your project is movable, shareable, and version-control-friendly out of the box.

The rule about missing-value codes matters because software treats a code such as −999 as a real number until it is told otherwise. The five-person example below shows what happens to the mean age when one person’s unknown age is stored as −999. The function tibble() builds a small table, and c() combines values into a vector.

library(tidyverse)

# Five participants; the third did not report age, coded -999
ages <- tibble(id = 1:5, age = c(34, 51, -999, 47, 62))
mean(ages$age)                       # -999 is treated as a real age

R has averaged the code with the real ages: (34 + 51 − 999 + 47 + 62) ÷ 5 = −805 ÷ 5 = −161 years, an impossible value. A code of 999 would have produced an equally wrong mean of 238.6 years, and a code hidden among many real values distorts the mean by less, which makes it harder to notice. Converting the code to NA fixes the problem:

# Fix 1: convert the code to NA, R's missing-value code
ages_fixed <- ages |> mutate(age = na_if(age, -999))
ages_fixed$age
mean(ages_fixed$age)                 # NA: R will not average an unknown value
mean(ages_fixed$age, na.rm = TRUE)   # drop the NA, then average the rest

na_if(age, -999) replaces every −999 in age with NA. A plain mean() then returns NA, because the average of a set that includes an unknown value is itself unknown; the argument na.rm = TRUE (“NA remove”) tells R to drop the missing value and average the rest: (34 + 51 + 47 + 62) ÷ 4 = 194 ÷ 4 = 48.5 years. The same conversion can be made when the file is first read, by listing every value that should count as missing:

# Fix 2: declare the code as missing when the file is read
demo_file <- tempfile(fileext = ".csv")   # a throwaway file for the demo
write_csv(ages, demo_file)
ages_read <- read_csv(demo_file, na = c("", "NA", "-999"),
                      show_col_types = FALSE)
ages_read$age

In practice. For a real file, the second fix is written as read_csv("data/raw/file.csv", na = c("", "NA", "-999")). The empty string "" covers blank cells and "NA" covers cells that already contain the letters NA. Every missing-value code used in the study belongs in this list and in the codebook.

The rule to use numbers for every variable comes from software in which text values were awkward to analyse. In R, the usual practice is to keep the numeric codes in the raw file, documented in the codebook, and to convert each coded variable into a factor in the cleaning script. A factor is R’s type for a categorical variable: it stores the codes, attaches a text label to each one, and fixes the order of the categories, so that tables and models display Female and Male in place of 1 and 2. A raw file that already stores text, such as the Yes and No values in cohort.csv, is also acceptable provided every value is spelled consistently, and the same factor() step then fixes the order of the categories.

# The raw file stores sex as 1 = Female, 2 = Male (see the codebook)
codes <- tibble(id = 1:4, sex = c(1, 2, 2, 1))
codes <- codes |>
  mutate(sex_f = factor(sex, levels = c(1, 2), labels = c("Female", "Male")))
codes
table(codes$sex_f)

Reading the code. In factor(sex, levels = c(1, 2), labels = c("Female", "Male")), levels lists the codes in the order wanted and labels gives the text for each code in the same order. In the printed table, <dbl> marks the original numeric column and <fct> marks the new factor. Keeping both columns side by side makes the conversion easy to check. The first level (here Female) becomes the reference category when the factor is used in a regression model.

Double data entry produces two files that should be identical, and the comparison step takes only a few lines of code. In the example below, two people entered the same five forms (sbp is systolic blood pressure in mmHg). The function anti_join(x, y, by = ...) returns the rows of x that have no exact match in y on the listed columns, so running it in both directions lists every record on which the two entries disagree.

# The same five forms, entered independently by two people
entry_a <- tibble(id = 1:5, age = c(34, 51, 29, 47, 62),
                  sbp = c(128, 141, 117, 135, 150))
entry_b <- tibble(id = 1:5, age = c(34, 15, 29, 47, 62),
                  sbp = c(128, 141, 117, 153, 150))

# Rows of entry_a with no exact match in entry_b, then the reverse
anti_join(entry_a, entry_b, by = c("id", "age", "sbp"))
anti_join(entry_b, entry_a, by = c("id", "age", "sbp"))

Reading the output. Records 2 and 4 disagree. In record 2 the first person typed an age of 51 and the second typed 15, a pair of swapped digits; in record 4 the systolic pressure is 135 in one entry and 153 in the other. The analyst checks each disagreement against the original paper form and corrects the wrong entry. An error made identically by both people cannot be detected this way, so double entry reduces errors substantially without removing every one.

The "save a new version, don't overwrite" rule is automatic if your transformations live in a script. The script, together with the log line it writes, becomes the file log. Download bp01.csv (the blood pressure study from the file log above: 1092 records, one per participant, 8 variables, 45 records incomplete) and save it as data/raw/bp01.csv in the project skeleton built earlier. For brevity, the script carries out the second and third steps of the log in one pass, so the file it saves as bp02.csv already contains the four derived variables (1047 observations and 12 variables); in a larger project each step would get its own numbered file, as in the log.

library(tidyverse)

# Read raw, never modify in place
bp_raw <- read_csv("data/raw/bp01.csv")

# Tidy: drop incomplete rows, build derived variables, lock factors
bp_clean <- bp_raw |>
  drop_na(systolic, diastolic, age) |>
  mutate(
    age_ct    = age - mean(age),                              # centred
    age_ctsq  = age_ct^2,                                       # quadratic term
    age_c3    = cut(age, c(0, 35, 55, Inf),
                    labels = c("young", "middle", "older")),
    htn       = factor(systolic >= 140 | diastolic >= 90,
                       levels = c(FALSE, TRUE),
                       labels = c("normotensive", "hypertensive"))
  )

# Persist as a new versioned file - and a small log line
write_csv(bp_clean, "data/processed/bp02.csv")
cat("bp02.csv", format(Sys.Date()), nrow(bp_clean), "obs",
    "\n", file = "data/file_log.txt", append = TRUE)

## At any time you can rebuild bp02 from bp01 by re-running this script.

What each function does.

What htn means. The variable htn (short for hypertension, or high blood pressure) classifies each participant as hypertensive if the systolic pressure is 140 mmHg or higher or the diastolic pressure is 90 mmHg or higher, and as normotensive (normal blood pressure) otherwise. Systolic pressure is the higher number in a reading such as 140/90 and diastolic the lower. The 140/90 threshold is a widely used clinical definition of hypertension; some guidelines use lower thresholds, which is one reason the activity below asks what happens when the cut-off moves. This simple definition ignores on_treatment, so a participant whose blood pressure is controlled by medication is classified as normotensive. A real analysis would decide how to classify treated participants and record that decision in the codebook.

Checking the result. The lines below confirm the row and column counts, show the derived variables for the first three participants, and produce the percentages quoted in the activity answers.

nrow(bp_raw)           # rows before cleaning
nrow(bp_clean)         # rows after drop_na()
ncol(bp_clean)         # 8 original variables + 4 derived
mean(bp_clean$age)     # the value subtracted to make age_ct

head(bp_clean[, c("id", "age", "age_ct", "age_ctsq", "age_c3")], 3)

# Prevalence of hypertension under the 140/90 rule
count(bp_clean, htn) |> mutate(percent = round(100 * n / sum(n), 1))

# The same count with the systolic cut-off raised to 150
bp_clean |>
  mutate(htn150 = systolic >= 150 | diastolic >= 90) |>
  count(htn150) |>
  mutate(percent = round(100 * n / sum(n), 1))

Reading the output. The 45 incomplete records have gone (1092 to 1047 rows) and four variables have been added (12 columns). count() counts the rows in each category, and the added percent column divides each count by the total: 313 of 1047 participants (29.9%) are hypertensive under the 140/90 rule, and 251 (24.0%, printed as 24) under a 150/90 rule. Centring is easiest to see in the first rows. The mean age is 48.05 years, so participant BP0001, aged 49, has age_ct = 49 − 48.05 = 0.95, and BP0002, aged 43, has age_ct = 43 − 48.05 = −5.05. A centred value of 0 identifies a participant of exactly average age, so in a regression model that uses age_ct, the intercept is the predicted outcome for a participant of average age in this sample. The squared term grows quickly away from the mean: BP0003, aged 82, has age_ct = 33.95 and age_ctsq = 33.95 × 33.95 ≈ 1153.

SPSS is a statistics package that is often used through menus. A workflow built only by clicking through menus, in SPSS or any similar program, is difficult to re-derive unless the underlying commands are also saved (SPSS calls them syntax). A script can be re-run six months from now, by a colleague, on a different computer.

The verification steps listed above translate into a handful of R functions. The code below runs them on the raw blood pressure file from the previous section, one check at a time, and the output under each piece of code is what it produces. The table after the code explains how to read each check and what a problem would look like.

library(tidyverse)
bp_raw <- read_csv("data/raw/bp01.csv", show_col_types = FALSE)

# 1. Range, quartiles, mean and number missing for each continuous variable
summary(bp_raw[, c("age", "bmi", "systolic", "diastolic")])

# 2. Number of missing values (NA) in every variable
colSums(is.na(bp_raw))

# 3. The five smallest and five largest values
head(sort(bp_raw$systolic), 5)
tail(sort(bp_raw$systolic), 5)
head(sort(bp_raw$age), 5)
tail(sort(bp_raw$age), 5)

# 4. Histogram: does the shape look reasonable?
hist(bp_raw$systolic, breaks = 30, xlab = "Systolic BP (mmHg)", main = "")

The histogram appears in the Plots pane; panel A of the figure below shows it.

# 5. Frequency table for each categorical variable, showing NA if any
table(bp_raw$sex, useNA = "ifany")
table(bp_raw$smoker, useNA = "ifany")
table(bp_raw$on_treatment, useNA = "ifany")

# 6. Consistency check: diastolic pressure should never exceed systolic
bp_raw |> filter(diastolic > systolic)

# 7. Every identifier should appear exactly once (compare with nrow)
n_distinct(bp_raw$id)

How to read each check.

The raw file passes every check above, so the example below plants three typical errors on a copy of the data, leaving the raw file untouched, and shows how the same checks reveal them and how the cleaning script fixes them.

bp_bad <- bp_raw                       # a copy; bp_raw is untouched
bp_bad$age[c(5, 120, 640)] <- 999      # planted: unknown age typed as 999
bp_bad$sex[c(10, 302)] <- "F"          # planted: abbreviation
bp_bad$sex[450] <- "male"              # planted: lower-case spelling

summary(bp_bad$age)
tail(sort(bp_bad$age), 5)
table(bp_bad$sex, useNA = "ifany")

The maximum age of 999 and a mean of 50.69 (up from 48.05) reveal the planted code, and the frequency table shows four categories where the codebook allows two. Both problems are fixed in the cleaning script:

# The fix belongs in the cleaning script, never in the raw file
bp_fixed <- bp_bad |>
  mutate(age = na_if(age, 999),
         sex = case_when(sex %in% c("Female", "F") ~ "Female",
                         sex %in% c("Male", "male") ~ "Male"))
summary(bp_fixed$age)
table(bp_fixed$sex, useNA = "ifany")

Reading the fix. na_if(age, 999) turns the code into NA. case_when() checks each condition in turn: %in% asks whether the value is one of the listed spellings, and the value after ~ is assigned when it is. Any value that matches no condition becomes NA, so a spelling that the analyst did not anticipate appears as missing in the next frequency table and is never silently relabelled. The fixed counts, 554 women and 538 men, match the raw file. After the fix, age has 14 missing values (11 original and 3 converted); where possible, the three ages are looked up on the original forms and entered through the script.

Variable labels can be stored in R in several ways. The simplest, and the one that works with every package, is a codebook table saved beside the data. The function tribble() builds a small table row by row: the first line names the columns, each name starting with ~, and every following line describes one variable. Category labels are attached separately with factor(), as shown in an earlier section.

codebook <- tribble(
  ~variable,      ~label,                             ~units_or_codes,
  "id",           "Participant identifier",           "BP0001 to BP1092",
  "age",          "Age",                              "years",
  "sex",          "Sex",                              "Female, Male",
  "smoker",       "Smoking status",                   "No, Yes",
  "bmi",          "Body mass index",                  "kg/m2",
  "systolic",     "Systolic blood pressure",          "mmHg",
  "diastolic",    "Diastolic blood pressure",         "mmHg",
  "on_treatment", "Takes blood-pressure medication",  "No, Yes"
)
codebook
write_csv(codebook, "data/codebook_bp01.csv")   # saved beside the data

Saving the codebook as a .csv file in the data folder keeps it with the data, and any collaborator can open it in R or in a spreadsheet. When a variable is added or redefined in the cleaning script, the codebook is updated in the same session.

The check for a Poisson model is to compare the mean of the counts with their variance. The invented example below records clinic visits in one year for ten patients at each of two clinics; var() computes the variance.

# Clinic visits in one year for ten patients at each of two clinics
visits_a <- c(0, 2, 4, 3, 2, 6, 1, 3, 2, 3)
visits_b <- c(0, 1, 1, 0, 9, 2, 0, 12, 1, 0)
mean(visits_a)
var(visits_a)
mean(visits_b)
var(visits_b)

Reading the output. Both clinics average 2.6 visits per patient. For clinic A, the squared distances from 2.6 are 6.76, 0.36, 1.96, 0.16, 0.36, 11.56, 2.56, 0.16, 0.36, and 0.16; they sum to 24.4, and 24.4 ÷ 9 = 2.71, because var() divides by n − 1 = 9. A variance this close to the mean is what a Poisson model expects. In clinic B, most patients have no visits or one visit and two have 9 and 12, so the variance (18.27) is about seven times the mean. That pattern is overdispersion, and a negative binomial model would be the safer starting point for these counts.

When an outcome has a long right tail, taking the natural logarithm of each value, log() in R, often makes the distribution close to symmetric, because the logarithm shrinks large values much more than small ones. The simulated example below uses hospital length of stay in days. The function rlnorm() draws right-skewed values for the simulation, and round(..., 1) keeps one decimal place.

# Simulated hospital length of stay (days) for 1000 patients
set.seed(4104)
los <- round(rlnorm(1000, meanlog = 1.2, sdlog = 0.8), 1)
summary(los)        # mean well above median: a sign of right skew
summary(log(los))   # after the log transformation
hist(los, breaks = 40)          # before: long right tail
hist(log(los), breaks = 30)     # after: roughly symmetric

Reading the output. Before the transformation the mean (4.57 days) is well above the median (3.40 days), because the few very long stays pull the mean upwards; this gap is a quick numerical sign of right skew. After the transformation the mean (1.19) and median (1.22) are close, matching the symmetric histogram. In a real analysis the new variable is created in the cleaning script, for example mutate(log_los = log(los)), and the transformation is recorded in the codebook. The logarithm is defined only for values above zero, so an outcome that contains zeros needs adjustment first; the next lesson returns to transformations.

The first question about any predictor is how much of it is missing. The two lines below give the count and the percentage of missing values for every variable in the raw blood pressure file. colMeans(is.na()) gives the proportion missing, because R counts each TRUE as 1 and each FALSE as 0.

colSums(is.na(bp_raw))                    # number missing per variable
round(100 * colMeans(is.na(bp_raw)), 1)   # percentage missing per variable

The decision that follows belongs in the analysis log in plain words. An entry for this file might read as follows.

The analyst types the log entry by hand, using the counts from the output above; 45 of 1092 participants is 4.1%.

The cross-tabulation below splits the cleaned blood pressure data into ten-year age bands and counts normotensive and hypertensive participants in each. The breaks c(17, 29, 39, ...) work because ages are whole numbers, so the band above 17 and up to 29 holds ages 18 to 29.

# bp_clean comes from the recoding pipeline in the section on data coding
bp_dec <- bp_clean |>
  mutate(age_dec = cut(age, c(17, 29, 39, 49, 59, 69, 79, Inf),
                       labels = c("18-29", "30-39", "40-49", "50-59",
                                  "60-69", "70-79", "80+")))
table(bp_dec$age_dec, bp_dec$htn)

# Merge the two oldest bands into 70+
bp_dec <- bp_dec |>
  mutate(age_dec6 = cut(age, c(17, 29, 39, 49, 59, 69, Inf),
                        labels = c("18-29", "30-39", "40-49", "50-59",
                                   "60-69", "70+")))
table(bp_dec$age_dec6, bp_dec$htn)

Reading the output. Before merging, the 80+ band holds only 17 participants, and only 3 of them are normotensive. A logistic model that compares each band with a reference band would estimate the 80+ effect from those 3 people, which gives a very wide confidence interval: with 50–59 as the reference band, the odds of hypertension in the 80+ band are estimated at about 11 times the odds in the reference band, with a 95% confidence interval from about 3 to more than 38. A band with no normotensive participants at all would prevent the model from converging (the fitting procedure cannot settle on an answer). Merging 70–79 and 80+ into 70+ gives 24 normotensive and 49 hypertensive participants, enough for a stable estimate. The 18–29 band has only 4 hypertensive participants and might likewise be merged with 30–39. Each merge is a judgement made before the outcome model is fitted, and it is recorded in the log.

Multilevel data are usually stored in long format, with one row per measurement and columns that identify the person and the centre. The small table below, typed into R for illustration with tribble(), holds ten blood pressure measurements on five people in two centres.

# Long format: one row per measurement (illustrative data typed into R)
bp_long <- tribble(
  ~centre, ~person, ~visit, ~sbp,
  "A",     "A01",   1,      132,
  "A",     "A01",   2,      128,
  "A",     "A02",   1,      141,
  "A",     "A02",   2,      139,
  "A",     "A02",   3,      137,
  "B",     "B01",   1,      118,
  "B",     "B01",   2,      121,
  "B",     "B02",   1,      150,
  "B",     "B03",   1,      126,
  "B",     "B03",   2,      124
)

Three counts characterise the hierarchy. count(bp_long, centre, person) counts the rows for each combination of centre and person, which is the number of measurements per person; distinct() keeps one row per person, so counting those rows by centre gives the number of people per centre; and n_distinct() counts the different person identifiers.

count(bp_long, centre, person)                        # measurements per person
bp_long |> distinct(centre, person) |> count(centre)  # people per centre
n_distinct(bp_long$person)                            # people overall

Reading the output. People have between one and three measurements, centre A has two people and centre B three, and there are five people in total. Here every person identifier begins with the centre letter, so each person is uniquely identified across centres. If two centres had both numbered their patients 01, 02, 03, the analyst would build a unique identifier first, for example with mutate(uid = paste(centre, person)), so that no merging step could confuse two people.

The correlation coefficient r summarises how closely two continuous variables follow a straight line. cor() computes it, and cor.test() adds a 95% confidence interval and a p-value.

cor(bp_clean$age, bp_clean$systolic)        # correlation coefficient r
cor.test(bp_clean$age, bp_clean$systolic)   # adds a 95% CI and a p-value

The scatterplot shows the shape of the relationship, and a simple linear regression puts a number on the slope. confint() gives 95% confidence intervals for the coefficients.

plot(systolic ~ age, data = bp_clean)         # scatterplot
fit1 <- lm(systolic ~ age, data = bp_clean)   # simple linear regression
coef(fit1)
confint(fit1)                                 # 95% CIs for the coefficients

Reading the output. The correlation is r = 0.48: systolic pressure tends to be higher in older participants, with considerable scatter. As a rough guide, values of |r| below about 0.3 are described as weak, values from about 0.3 to 0.7 as moderate, and values above about 0.7 as strong; the labels are conventions, and the scatterplot should always be inspected as well. In the cor.test() output, the 95% confidence interval for r runs from 0.43 to 0.52, and the p-value (< 2.2e-16) means that a correlation this large would almost never arise by chance if age and systolic pressure were unrelated; t and df are the test statistic and its degrees of freedom. The regression slope for age is 0.585: on average, systolic pressure is 0.59 mmHg higher for each additional year of age (95% CI 0.52 to 0.65). The intercept, 99.4 mmHg, is the predicted value at age 0, which lies far outside the data and has no practical meaning; centred age would give an intercept for a participant of average age.

A model reporting sentence. “Systolic blood pressure was moderately correlated with age (r = 0.48, 95% CI 0.43 to 0.52); in a simple linear regression, each additional year of age was associated with a 0.59 mmHg higher systolic pressure (95% CI 0.52 to 0.65), before adjustment for other variables.”

For two groups, the group means and a simple linear regression answer the question together. tapply(x, group, mean) applies mean() to x separately within each group.

table(bp_clean$smoker)                             # group sizes
tapply(bp_clean$systolic, bp_clean$smoker, mean)   # mean systolic by group
fit2 <- lm(systolic ~ smoker, data = bp_clean)
coef(fit2)
confint(fit2)

For three or more groups, one-way ANOVA tests whether the group means differ. The example uses the three age bands created earlier.

# Three or more groups: one-way ANOVA across the age bands
tapply(bp_clean$systolic, bp_clean$age_c3, mean)
summary(aov(systolic ~ age_c3, data = bp_clean))

Reading the output. Mean systolic pressure is 126.3 mmHg among the 811 non-smokers and 131.5 mmHg among the 236 smokers. In the regression, the intercept (126.3) is the mean for the reference group, non-smokers, because No is the first level of smoker, and the coefficient smokerYes (5.15) is the difference between the two group means, with a 95% CI of 2.57 to 7.73 mmHg that excludes zero. In the ANOVA table, Df gives the degrees of freedom (number of groups minus 1, and number of people minus number of groups), Sum Sq and Mean Sq measure the variation between and within the groups, F value is the ratio of between-group to within-group variation, and Pr(>F) is the p-value. An F of 111.4 with p < 2e-16 indicates that mean systolic pressure differs across the age bands; ANOVA does not identify which bands differ, so the group means (115.6, 126.8, and 136.8 mmHg) are reported alongside it.

A model reporting sentence. “Mean systolic pressure was 131.5 mmHg among smokers and 126.3 mmHg among non-smokers, an unadjusted difference of 5.1 mmHg (95% CI 2.6 to 7.7). Mean systolic pressure also differed across age bands (one-way ANOVA, F = 111.4 on 2 and 1044 degrees of freedom, p < 0.001), rising from 115.6 mmHg in participants aged 35 or younger to 136.8 mmHg in those older than 55.”

A cross-tabulation counts the participants in every combination of two categorical variables, and row percentages make the groups comparable. prop.table(tab, margin = 1) divides each count by its row total.

tab <- table(smoker = bp_clean$smoker, htn = bp_clean$htn)
tab                                            # counts
round(100 * prop.table(tab, margin = 1), 1)    # row percentages

The chi-squared test asks whether the pattern of counts differs from what would be expected if the two variables were unrelated. The expected counts are worth printing, because the test is unreliable when any expected count is small.

chisq.test(tab)                                # chi-squared test
chisq.test(tab)$expected                       # expected counts if unrelated

Reading the output. Hypertension is more common among smokers (35.2%, 83 of 236) than among non-smokers (28.4%, 230 of 811). The chi-squared statistic is 3.73 on 1 degree of freedom with p = 0.054. R applies a small correction (Yates’ continuity correction) automatically for a table with two rows and two columns. A p-value of 0.054 indicates weak evidence against no association: the data are somewhat unusual under no association, and a difference at least this large would arise in about 5 of every 100 samples of this size if there were no association. The crude difference will also change once age, a likely confounder, is taken into account. All expected counts are far above 5, so the test is appropriate here; when any expected count falls below 5, Fisher’s exact test, fisher.test(tab), is used instead.

A model reporting sentence. “Hypertension was more common among smokers (35.2%) than among non-smokers (28.4%), although the evidence for an unadjusted association was weak (chi-squared test, p = 0.054).”

A correlation matrix shows the correlation between every pair of continuous variables at once. cor() applied to several columns returns the matrix, and round(..., 2) keeps two decimal places.

round(cor(bp_clean[, c("age", "bmi", "systolic", "diastolic")]), 2)

Reading the output. Each cell is the correlation between the row variable and the column variable, and the diagonal is always 1. A common rule of thumb flags pairs with |r| above about 0.8 as possible collinearity problems. Systolic and diastolic pressure (r = 0.84) cross that line, which is expected because both measure blood pressure; a model would normally include one of them, or a combined measure, as a predictor. The correlation matrix helps detect collinearity, whereas the choice of confounders to include in a model comes from the causal diagram.