# ============================================================================= # HSCI 410 Public Health Assessment and Analysis - Lesson 1: A Structured # Approach to Data Analysis # Answer key for the in-lesson R activities # Data file(s): cohort.csv, bp01.csv (download from the lesson page; save in # your R working directory. The script copies each file into the # data/raw/ folder that the lesson's project skeleton creates.) # Packages: mediation, dagitty, tidyverse (dplyr, readr, tidyr, ggplot2), here # (install once with install.packages(c("mediation", "dagitty", "tidyverse", "here"))) # Reproduces every code block in the lesson, then answers each activity question. # ============================================================================= # ==== Section 1: Introduction & Data Collection / Fitting a mediation model in R (Baron & Kenny + the mediation package) ==== # 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) # ---- Activity questions, Section 1 ------------------------------------------ c_total <- unname(coef(lm(health ~ education, data = dat))["education"]) c_direct <- unname(coef(fit_Y)["education"]) cat("\nQ1: total effect c =", round(c_total, 3), "; direct effect c' =", round(c_direct, 3), "; difference c - c' =", round(c_total - c_direct, 3), "\n") cat(" The total effect is the larger one. The gap between c and c' is the\n", " indirect effect that travels through income (about half of the total).\n") cat("\nQ2: a =", round(unname(a), 3), "; b =", round(unname(b), 3), "; a * b =", round(unname(a * b), 3), "\n") cat(" Bootstrapped ACME =", round(med$d0, 3), " 95% CI (", round(med$d0.ci[1], 3), ",", round(med$d0.ci[2], 3), ")\n") cat(" The product a*b equals the ACME to three decimals (the point estimate is\n", " the same quantity; the bootstrap only adds the CI), and the CI excludes 0.\n") cat("\nQ3: Prop. Mediated =", round(med$n0, 3), " (95% CI", round(med$n0.ci[1], 3), ",", round(med$n0.ci[2], 3), ")\n") cat(" About", round(100 * med$n0), "% of education's total effect on health is\n", " carried by income. If the ACME CI crossed zero we could not rule out that\n", " income contributes nothing to the education-health link in these data.\n") # ==== Section 1: Introduction & Data Collection / A reproducible project skeleton in RStudio ==== # NOTE (answer key): the lesson reads data/raw/cohort.csv. Download cohort.csv # from the lesson page into your working directory; the two lines after the # dir.create() calls copy it into the data/raw/ folder that the skeleton # creates, so the pipeline below runs exactly as printed on the page. # 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) # NOTE (answer key): place the downloaded raw file where the skeleton expects it file.copy("cohort.csv", "data/raw/cohort.csv", overwrite = 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 }") # ---- Activity questions, Section 2 ------------------------------------------ cat("\nQ1: folders created by the four dir.create() calls:\n") print(list.dirs(c("data", "R", "output"), recursive = TRUE)) cat(" data/raw/ holds the irreplaceable original; data/processed/ holds files\n", " that can be rebuilt from raw + code, so a cleaning mistake is never fatal.\n") cat("\nQ2: raw has", nrow(raw), "rows and", ncol(raw), "columns; clean has", nrow(clean), "rows and", ncol(clean), "columns\n") cat(" filter() dropped", nrow(raw) - nrow(clean), "rows with a missing outcome\n") cat(" brand-new column(s):", setdiff(names(clean), names(raw)), "\n") print(table(clean$age_grp)) cat(" re-encoded in place: smoker (character -> factor with levels", paste(levels(clean$smoker), collapse = " < "), ")\n") cat("\nQ3: here() resolved the raw file to:", here("data/raw/cohort.csv"), "\n") cat(" The path is built from the project root, so the same script runs on any\n", " machine, user account or operating system once the project is unzipped.\n") # ==== Section 2: Data Coding, Entry & File Management / A reproducible recoding pipeline (no overwrites, ever) ==== # NOTE (answer key): download bp01.csv from the lesson page into your working # directory; this line copies it into data/raw/ (created above). file.copy("bp01.csv", "data/raw/bp01.csv", overwrite = TRUE) 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. # ---- Activity questions, Section 3 ------------------------------------------ cat("\nQ1: bp01 had", nrow(bp_raw), "obs and", ncol(bp_raw), "vars; bp02 has", nrow(bp_clean), "obs and", ncol(bp_clean), "vars (", nrow(bp_raw) - nrow(bp_clean), "incomplete rows dropped )\n") cat(" The four new variables and their classes:\n") str(bp_clean[, c("age_ct", "age_ctsq", "age_c3", "htn")]) cat(" age_ct: continuous (age minus the mean age of", round(mean(bp_clean$age), 1), ")\n", " age_ctsq: derived continuous (age_ct squared, for a curved age effect)\n", " age_c3: categorical factor with three age bands\n", " htn: derived binary factor (140/90 rule)\n") print(table(bp_clean$age_c3)) print(table(bp_clean$htn)) prev_140 <- mean(bp_clean$systolic >= 140 | bp_clean$diastolic >= 90) prev_150 <- mean(bp_clean$systolic >= 150 | bp_clean$diastolic >= 90) cat("\nQ2: prevalence of 'hypertensive' with systolic cutoff 140:", round(100 * prev_140, 1), "%\n") cat(" prevalence with systolic cutoff 150:", round(100 * prev_150, 1), "%\n") cat(" Raising the cutoff lowers the prevalence by", round(100 * (prev_140 - prev_150), 1), "percentage points; a categorical\n", " recode is only as defensible as its threshold.\n") cat("\nQ3: files now on disk:\n") cat(" data/raw/bp01.csv exists:", file.exists("data/raw/bp01.csv"), "; data/processed/bp02.csv exists:", file.exists("data/processed/bp02.csv"), "\n") cat(" file log contents:\n") cat(readLines("data/file_log.txt"), sep = "\n") cat(" Because bp01.csv is untouched, a recoding mistake in bp02 (for example\n", " reversed htn labels) is repaired by fixing the script and re-running it;\n", " an overwritten SPSS working file cannot be re-derived.\n") # ==== Supplementary code shown in the reading ==== # The blocks below reproduce the additional R code boxes in the lesson, in the # order in which they appear. They assume the code above has been run first # (bp_raw, bp_clean, raw, clean, fit_M and the project folders must exist). # ---- Section 1: Finding the adjustment set with dagitty ---- 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") # ---- Section 1: Worked Example: a collider created by studying only hospital patients ---- 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) # ---- Section 1: Mediation, printing path a ---- coef(fit_M) # the education coefficient is path a # ---- Section 1: Project skeleton, checking the data before and after cleaning ---- glimpse(raw) # one line per variable: name, type, first values nrow(raw) # rows before filtering nrow(clean) # rows after dropping participants with no outcome count(clean, age_grp) # participants in each age band # 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")) # ---- Section 2: Worked Example: how a -999 code distorts a mean, and the fix ---- 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 # 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 # 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 # ---- Section 2: Numeric codes and R factors ---- # 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) # ---- Section 2: Worked Example: comparing two independent entries in R ---- # 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")) # ---- Section 2: Recoding pipeline, checking the result ---- 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)) # ---- Section 3: Verifying bp01.csv in R ---- 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 = "") # 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) # ---- Section 3: Worked Example: finding and fixing planted errors ---- 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 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") # ---- Section 3: Keeping a codebook in R ---- 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 # ---- Section 4: Worked Example: comparing the variance with the mean for counts ---- # 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) # ---- Section 4: Worked Example: before and after a log transformation ---- # 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 # ---- Section 4: Tabulating missing values ---- colSums(is.na(bp_raw)) # number missing per variable round(100 * colMeans(is.na(bp_raw)), 1) # percentage missing per variable # ---- Section 4: Worked Example: combining a sparse category ---- # 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) # ---- Section 4: A long-format multilevel table in R ---- # 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 ) 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 # ---- Section 4: Worked Example: two continuous variables ---- 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 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 # ---- Section 4: Worked Example: one continuous and one categorical variable ---- 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) # 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)) # ---- Section 4: Worked Example: two categorical variables ---- tab <- table(smoker = bp_clean$smoker, htn = bp_clean$htn) tab # counts round(100 * prop.table(tab, margin = 1), 1) # row percentages chisq.test(tab) # chi-squared test chisq.test(tab)$expected # expected counts if unrelated # ---- Section 4: Checking predictors for collinearity with a correlation matrix ---- round(cor(bp_clean[, c("age", "bmi", "systolic", "diastolic")]), 2)