# =============================================================================
# HSCI 410 Public Health Assessment and Analysis  -  Lesson 2: Data Cleaning
# and Descriptive Analyses
# Answer key for the in-lesson R activities
# Data file(s): phaa_survey.csv  (download from the lesson page; save in your
#               R working directory). The last block writes phaa_survey_clean.csv,
#               the analytic file that Lessons 3 and 4 read.
# Packages: psych, mice  (install once with install.packages(c("psych", "mice")))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================


# ==== Section 2: Data Cleaning Strategies / Activity: read, fix classes, and clean impossible values ====

# 1. Read the raw file. Treat blank cells as missing (na.strings).
phaa <- read.csv(file = "phaa_survey.csv",
                 stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))
summary(phaa$age)                 # before the conversion in step 2

# (answer key) keep the state before conversion so Q1 can compare before/after
age_before <- phaa$age
class(age_before)

# 2. Make sure age and bmi are numeric (in a real export they often arrive
#    as character because a cell holds text such as "unknown"; here every
#    cell is a number, so no NAs are created at this step).
phaa$age <- as.numeric(phaa$age)
phaa$bmi <- as.numeric(phaa$bmi)
summary(phaa$age)                 # after the conversion

# (answer key) look at the impossible values before the range checks remove them
phaa[!is.na(phaa$age) & (phaa$age < 18 | phaa$age > 100), c("id", "age")]
phaa[!is.na(phaa$systolic_bp) & (phaa$systolic_bp < 80 | phaa$systolic_bp > 220),
     c("id", "systolic_bp")]
phaa[!is.na(phaa$bmi) & (phaa$bmi < 13 | phaa$bmi > 60), c("id", "bmi")]
n_sbp_out <- sum(phaa$systolic_bp < 80 | phaa$systolic_bp > 220, na.rm = TRUE)

# 3a. Count and list the rows that break each hard limit BEFORE changing them
sum(phaa$age < 18 | phaa$age > 100, na.rm = TRUE)
phaa$id[which(phaa$age < 18 | phaa$age > 100)]
sum(phaa$systolic_bp < 80 | phaa$systolic_bp > 220, na.rm = TRUE)
phaa$id[which(phaa$systolic_bp < 80 | phaa$systolic_bp > 220)]
sum(phaa$bmi < 13 | phaa$bmi > 60, na.rm = TRUE)
phaa$id[which(phaa$bmi < 13 | phaa$bmi > 60)]

# 3b. Range checks: replace the impossible values with NA
phaa$age[phaa$age < 18 | phaa$age > 100] <- NA
phaa$systolic_bp[phaa$systolic_bp < 80 | phaa$systolic_bp > 220] <- NA
phaa$bmi[phaa$bmi < 13 | phaa$bmi > 60] <- NA
summary(phaa$age)                 # after the range check

# 4. Put the education levels in order (lowest to highest) so that
#    comparisons such as <= are meaningful
phaa$education <- factor(phaa$education,
  levels = c("Less than high school", "High school", "Some college",
             "Bachelor's", "Graduate degree"),
  ordered = TRUE)
levels(phaa$education)

# ---- Activity questions, Section 1 ------------------------------------------
cat("\nQ1: class of age as read from the file:", class(age_before), "\n")
cat("    NAs in age before as.numeric():", sum(is.na(age_before)),
    "; after as.numeric():", sum(is.na(as.numeric(age_before))), "\n")
cat("    NAs in age after the range check (18-100):", sum(is.na(phaa$age)), "\n")
cat("    Every raw age is a number (the two impossible entries are -3 and 220),\n",
    "   so the column is already integer and as.numeric() creates no NAs; the\n",
    "   two NAs come from the range check. A raw file with text such as\n",
    "   'unknown' would arrive as character and every such cell would become NA.\n")

cat("\nQ2:", n_sbp_out, "rows of systolic_bp were set to NA by the 80-220 range check",
    "(the values 300 and 60 in rows P0017 and P0245)\n")
cat("    Replacing with NA keeps the other", ncol(phaa) - 1,
    "variables of those rows in the analysis.\n")

cat("\nQ3: levels(phaa$education):\n")
print(levels(phaa$education))
cat("    is.ordered:", is.ordered(phaa$education), "\n")
cat("    alphabetical order R would use by default:\n")
print(sort(unique(as.character(phaa$education))))
print(table(phaa$education))


# ==== Section 3: Descriptive Analyses / Activity: descriptives and standard plots ====

# Categorical descriptives ---------------------------------------------------
table(phaa$gender)                            # frequencies
prop.table(table(phaa$gender))                # proportions
round(prop.table(table(phaa$gender)) * 100, 1)  # %

# Cross-tab: gender by smoker, row percentages
round(prop.table(table(phaa$gender, phaa$smoker), margin = 1) * 100, 1)

# Numeric descriptives --------------------------------------------------------
summary(phaa$age)
sd(phaa$age, na.rm = TRUE)
IQR(phaa$age, na.rm = TRUE)

# psych::describe() summarises many variables at once
library(psych)
describe(phaa[, c("age", "bmi", "systolic_bp", "phys_act_min")])

# Plots: histogram, boxplot stratified by exposure, bar chart, scatter plot ---
hist(phaa$systolic_bp, main = "Systolic BP", xlab = "mmHg")

boxplot(systolic_bp ~ smoker, data = phaa,
        main = "Systolic BP by smoking status", ylab = "mmHg")

barplot(table(phaa$gender), main = "Gender distribution")

plot(phaa$age, phaa$systolic_bp, main = "Systolic BP by age",
     xlab = "Age (years)", ylab = "Systolic BP (mmHg)")

# ---- Activity questions, Section 2 ------------------------------------------
smk_pct <- round(prop.table(table(phaa$gender, phaa$smoker), margin = 1) * 100, 1)
cat("\nQ1: percentage of current smokers within each gender:\n")
print(smk_pct[, "Yes"])
cat("    highest:", names(which.max(smk_pct[, "Yes"])), "(",
    max(smk_pct[, "Yes"]), "% ); lowest:", names(which.min(smk_pct[, "Yes"])), "(",
    min(smk_pct[, "Yes"]), "% ); gap =",
    round(max(smk_pct[, "Yes"]) - min(smk_pct[, "Yes"]), 1), "percentage points\n")
print(table(phaa$gender, phaa$smoker))   # the Non-binary row rests on few people

desc <- describe(phaa[, c("age", "bmi", "systolic_bp", "phys_act_min")])
cat("\nQ2: skew by variable:\n")
print(round(desc[, c("n", "mean", "sd", "median", "skew", "kurtosis")], 2))
cat("    largest |skew|:", rownames(desc)[which.max(abs(desc$skew))],
    "; closest to zero (most nearly normal):", rownames(desc)[which.min(abs(desc$skew))], "\n")

cat("\nQ3: median systolic BP by smoking status:\n")
print(tapply(phaa$systolic_bp, phaa$smoker, median, na.rm = TRUE))
cat("    mean systolic BP by smoking status:\n")
print(round(tapply(phaa$systolic_bp, phaa$smoker, mean, na.rm = TRUE), 1))
out_no  <- boxplot.stats(phaa$systolic_bp[phaa$smoker == "No"])$out
out_yes <- boxplot.stats(phaa$systolic_bp[phaa$smoker == "Yes"])$out
cat("    boxplot outliers beyond the whiskers: non-smokers", length(out_no),
    "( range", paste(range(out_no), collapse = "-"), ") ; smokers", length(out_yes),
    if (length(out_yes) > 0) paste("( range", paste(range(out_yes), collapse = "-"), ")") else "", "\n")
print(t.test(systolic_bp ~ smoker, data = phaa))
print(wilcox.test(systolic_bp ~ smoker, data = phaa))


# ==== Section 3: Descriptive Analyses / Activity: Cronbach's alpha and exploratory factor analysis ====

# 1. Internal consistency for each candidate scale ----------------------------
library(psych)
dep_items <- phaa[, c("dep1","dep2","dep3","dep4","dep5","dep6","dep7")]
anx_items <- phaa[, c("anx1","anx2","anx3","anx4","anx5")]

alpha(dep_items)             # raw_alpha >= 0.70 is acceptable, >= 0.80 is good
alpha(anx_items)
dep_alpha <- alpha(dep_items)            # store the result to read it at 3 decimals
round(dep_alpha$total$raw_alpha, 3)
round(dep_alpha$alpha.drop[, "raw_alpha", drop = FALSE], 3)  # alpha if an item is dropped

# 2. How many underlying factors? ------------------------------------------
fa.parallel(dep_items, fa = "fa")   # parallel analysis for factors

# 3. Inspect a one-factor solution for the depression items ----------------
factanal(x = na.omit(dep_items), factors = 1)
# na.omit(): factanal() needs complete rows and dep4 has 15 NAs, so these 15
# people are left out of this check (a complete-case analysis).
# With a single factor there is nothing to rotate, so no rotation is requested.
# Loadings >= 0.40 contribute meaningfully to the factor.

# 4. Two-factor solution combining dep + anx items -------------------------
combined <- na.omit(cbind(dep_items, anx_items))
factanal(x = combined, factors = 2, rotation = "varimax")
# dep1-7 should load on one factor; anx1-5 on the other.

# 5. Build derived scale variables for use in later lessons ------------------
# Score = mean of the answered items x number of items, calculated only when
# at most one item is missing (otherwise NA). A missing item is never
# counted as 0, which would understate the person's score.
n_dep <- rowSums(!is.na(dep_items))     # items answered by each person
n_anx <- rowSums(!is.na(anx_items))
phaa$dep_score <- ifelse(n_dep >= 6, rowMeans(dep_items, na.rm = TRUE) * 7, NA)
phaa$anx_score <- ifelse(n_anx >= 4, rowMeans(anx_items, na.rm = TRUE) * 5, NA)
summary(phaa$dep_score)

# 6. Save the cleaned, scale-augmented file under a new name -----------------
write.csv(phaa, "phaa_survey_clean.csv", row.names = FALSE)

# ---- Activity questions, Section 3 ------------------------------------------
a_dep <- alpha(dep_items)
a_anx <- alpha(anx_items)
cat("\nQ1: raw alpha, depression items:", round(a_dep$total$raw_alpha, 3),
    "; anxiety items:", round(a_anx$total$raw_alpha, 3), "\n")
cat("    Reliability if an item is dropped (depression):\n")
print(round(a_dep$alpha.drop[, c("raw_alpha", "std.alpha")], 3))
cat("    item-total correlations (r.drop):\n")
print(round(a_dep$item.stats$r.drop, 3))
cat("    item contributing least:", rownames(a_dep$alpha.drop)[which.max(a_dep$alpha.drop$raw_alpha)],
    "(highest alpha if dropped, lowest r.drop)\n")

fp <- fa.parallel(dep_items, fa = "fa", plot = FALSE)
cat("\nQ2: fa.parallel suggests", fp$nfact, "factor(s) and", fp$ncomp, "component(s)\n")
cat("    observed eigenvalues (FA):", round(fp$fa.values, 2), "\n")
cat("    simulated (parallel) eigenvalues (FA):", round(fp$fa.sim, 2), "\n")
f1 <- factanal(x = na.omit(dep_items), factors = 1)
# (answer key) with one factor, rotation = "varimax" changes nothing:
f1_v <- factanal(x = na.omit(dep_items), factors = 1, rotation = "varimax")
cat("    one-factor loadings identical with and without varimax:",
    isTRUE(all.equal(unclass(f1$loadings), unclass(f1_v$loadings))), "\n")
cat("    one-factor loadings:\n")
print(round(as.vector(f1$loadings), 2))
cat("    all seven >= 0.40:", all(f1$loadings >= 0.40), "; smallest =",
    round(min(f1$loadings), 2), "; proportion of variance explained =",
    round(sum(f1$loadings^2) / 7, 2), "\n")

f2 <- factanal(x = combined, factors = 2, rotation = "varimax")
cat("\nQ3: two-factor loadings (combined items):\n")
print(round(unclass(f2$loadings), 2))
L <- unclass(f2$loadings)
cat("    dep3 loads", round(L["dep3", 1], 2), "on Factor 1 and", round(L["dep3", 2], 2),
    "on Factor 2; anx2 loads", round(L["anx2", 1], 2), "on Factor 1 and",
    round(L["anx2", 2], 2), "on Factor 2\n")
cat("    largest cross-loading of any item:", round(max(apply(L, 1, min)), 2), "\n")
cat("    smallest loading of any item on its own factor:",
    round(min(apply(L, 1, max)), 2), "\n")
# (answer key) the scores are prorated (mean of answered items x number of items),
# so a skipped item is not counted as 0; compare with the simple row sum:
inc_dep <- !complete.cases(dep_items)
cat("    people missing one depression item:", sum(inc_dep),
    "; prorated minus simple row sum, range:",
    paste(round(range(phaa$dep_score[inc_dep] - rowSums(dep_items[inc_dep, ], na.rm = TRUE)), 2), collapse = " to "),
    "; mean:", round(mean(phaa$dep_score[inc_dep] - rowSums(dep_items[inc_dep, ], na.rm = TRUE)), 2), "\n")
cat("\nSaved phaa_survey_clean.csv with", nrow(phaa), "rows and", ncol(phaa), "columns\n")


# ==== Supplementary code shown in the reading ====
# These are the worked-code boxes from the lesson page, in page order. They do not
# change phaa_survey_clean.csv, which the scale activity above has already written.

# ---- Section 1: a first look at the raw file ----
# Read the raw file. Blank cells and the text "NA" become missing values (NA).
phaa <- read.csv("phaa_survey.csv", stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))

dim(phaa)                 # number of rows (people) and columns (variables)
str(phaa)                 # each variable's type and its first few values
summary(phaa[, c("age", "systolic_bp", "bmi", "dep1")])
table(phaa$smoker, useNA = "ifany")   # counts for a text (character) variable
colSums(is.na(phaa))      # how many values are missing in each column

# ---- Section 2: outlier rules and soft limits (raw data) ----
# Worked example: the eleven readings from the text
x <- c(96, 102, 106, 110, 112, 115, 118, 121, 125, 131, 184)
quantile(x, probs = c(0.25, 0.75))     # Q1 and Q3
IQR(x)                                  # Q3 minus Q1
boxplot.stats(x)$out                    # values beyond the fences

# The same steps on BMI in the course file (before any cleaning)
q     <- quantile(phaa$bmi, probs = c(0.25, 0.75), na.rm = TRUE)
iqr   <- IQR(phaa$bmi, na.rm = TRUE)
lower <- q[1] - 1.5 * iqr
upper <- q[2] + 1.5 * iqr
c(lower, upper)
flag <- which(phaa$bmi < lower | phaa$bmi > upper)
length(flag)                            # how many rows are flagged
phaa[flag, c("id", "bmi")]              # who they are, and their values

# z-scores: scale() computes (x - mean) / SD for every value
z <- as.numeric(scale(phaa$bmi))
phaa$id[which(abs(z) > 3)]              # IDs with |z| greater than 3
round(z[which(phaa$bmi == 31.6)], 2)   # z for the BMI of 31.6

# Masking: repeat without the impossible BMI of 95
bmi_ok <- phaa$bmi[phaa$bmi < 60]
round(c(mean = mean(phaa$bmi), sd = sd(phaa$bmi)), 2)
round(c(mean = mean(bmi_ok), sd = sd(bmi_ok)), 2)
round((31.6 - mean(bmi_ok)) / sd(bmi_ok), 2)

# Soft limits: count values that deserve a second look (nothing is changed)
sum(phaa$systolic_bp < 90 | phaa$systolic_bp > 180)
sum(phaa$bmi < 16 | phaa$bmi > 50)

# ---- Section 2: the boxes below follow the cleaning activity on the page, so the
#      cleaning steps are repeated here (without the printed checks) ----
phaa <- read.csv(file = "phaa_survey.csv", stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))
phaa$age[phaa$age < 18 | phaa$age > 100] <- NA
phaa$systolic_bp[phaa$systolic_bp < 80 | phaa$systolic_bp > 220] <- NA
phaa$bmi[phaa$bmi < 13 | phaa$bmi > 60] <- NA
phaa$education <- factor(phaa$education,
  levels = c("Less than high school", "High school", "Some college",
             "Bachelor's", "Graduate degree"),
  ordered = TRUE)

# ---- Section 2: counting and comparing missing values ----
# How much is missing, and where? (after the range checks above)
n_miss <- colSums(is.na(phaa))
n_miss[n_miss > 0]                       # only the columns with gaps
round(100 * n_miss[n_miss > 0] / nrow(phaa), 1)   # the same, as percentages
sum(complete.cases(phaa))                # rows with no gap in any column
nrow(na.omit(phaa))                      # na.omit() keeps only those rows

# Do people with and without a recorded income differ on what was observed?
inc_miss <- is.na(phaa$income)           # TRUE = income missing
table(inc_miss)
round(tapply(phaa$age, inc_miss, mean, na.rm = TRUE), 1)
round(tapply(phaa$systolic_bp, inc_miss, mean, na.rm = TRUE), 1)
round(prop.table(table(inc_miss, phaa$smoker), margin = 1) * 100, 1)

# ---- Section 2: a simulation of MCAR and MNAR missingness ----
set.seed(410)
income <- round(rlnorm(10000, meanlog = log(55), sdlog = 0.6))  # simulated, $1000s
mean(income)                            # the true mean in this population

# MCAR: everyone has the same 30% chance of leaving the question blank
inc_mcar <- income
inc_mcar[runif(10000) < 0.30] <- NA
mean(inc_mcar, na.rm = TRUE)

# MNAR: the top quarter of earners leave it blank far more often
p_miss <- ifelse(income > quantile(income, 0.75), 0.70, 0.15)
inc_mnar <- income
inc_mnar[runif(10000) < p_miss] <- NA
mean(inc_mnar, na.rm = TRUE)

# Standard error of a mean = SD / square root of n (smaller = more precise)
sd(income) / sqrt(10000)
sd(inc_mcar, na.rm = TRUE) / sqrt(sum(!is.na(inc_mcar)))

# ---- Section 2: converting missing-data codes to NA ----
age_raw <- c(34, -9, 51, 999, 47)    # -9 = refused, 999 = not asked
mean(age_raw)                        # wrong: the codes are averaged as if they were ages
age_raw[age_raw %in% c(-9, 999)] <- NA
age_raw
mean(age_raw, na.rm = TRUE)          # na.rm = TRUE leaves the NAs out
# The codes can also be converted while reading a file:
# mydata <- read.csv("mydata.csv", na.strings = c("", "NA", "-9", "999"))

# ---- Section 2 (going further): a minimal multiple imputation ----
library(mice)
mi_vars <- phaa[, c("age", "gender", "smoker", "bmi", "phys_act_min",
                    "social_support_score")]
mi_vars$gender <- factor(mi_vars$gender)
mi_vars$smoker <- factor(mi_vars$smoker)
imp <- mice(mi_vars, m = 5, seed = 2026, printFlag = FALSE)  # 5 completed datasets
fit <- with(imp, lm(phys_act_min ~ 1))   # lm(y ~ 1) estimates the mean of y
summary(pool(fit), conf.int = TRUE)      # pooled with Rubin's rules

# Complete-case estimate and its standard error, for comparison
mean(phaa$phys_act_min, na.rm = TRUE)
sd(phaa$phys_act_min, na.rm = TRUE) / sqrt(sum(!is.na(phaa$phys_act_min)))

# ---- Section 2: log(), log1p(), exp() and the reciprocal ----
crp <- c(0.4, 0.8, 1.1, 1.5, 2.3, 3.0, 4.8, 9.5, 21.0)  # a right-skewed biomarker (mg/L)
mean(crp)                  # arithmetic mean, pulled up by the 21.0
median(crp)
log_crp <- log(crp)        # natural log of each value
round(log_crp, 2)
mean(log_crp)              # the mean on the log scale
exp(mean(log_crp))         # back-transformed: the geometric mean, in mg/L

visits <- c(0, 0, 1, 2, 5, 12)
log(visits)                # log(0) is -Inf, which cannot be analysed
log1p(visits)              # log1p(x) = log(x + 1), and log1p(0) = 0
expm1(log1p(visits))       # expm1() undoes log1p()

1 / c(2, 5, 10)            # the reciprocal reverses the order
-1 / c(2, 5, 10)           # a minus sign restores it

# ---- Section 2: cleaning text, recoding, and checking for duplicates ----
city <- c("Vancouver", "vancouver ", "VANCOUVER", "Vancoouver", " Burnaby")
city <- trimws(tolower(city))    # remove stray spaces, then make lower case
city
city[grepl("^vanc", city)] <- "vancouver"   # regular expression: starts with "vanc"
table(city)

# Recoding with ifelse() and cut()
smoker01 <- ifelse(phaa$smoker == "Yes", 1, 0)   # 1 = smoker, 0 = non-smoker
table(phaa$smoker, smoker01)
bmi_group <- cut(phaa$bmi, breaks = c(0, 18.5, 25, 30, Inf), right = FALSE,
                 labels = c("Under 18.5", "18.5 to 24.9", "25 to 29.9", "30 or more"))
table(bmi_group, useNA = "ifany")

# Duplicates: is any ID present more than once?
sum(duplicated(phaa$id))

# ---- Section 2: the header of a cleaning script and a cleaning log ----
# ---------------------------------------------------------------
# clean_phaa.R   Cleaning script for phaa_survey.csv
# Raw file: phaa_survey.csv (never edited)
# Output:   phaa_survey_clean.csv (written at the end of the script)
# ---------------------------------------------------------------
cleaning_log <- data.frame(
  date     = "2026-09-30",
  analyst  = "AB",
  id       = c("P0033", "P0612", "P0017", "P0245", "P0078"),
  variable = c("age", "age", "systolic_bp", "systolic_bp", "bmi"),
  old      = c(-3, 220, 300, 60, 95),
  new      = NA,
  reason   = c("below 18: impossible", "above 100: impossible",
               "above 220: implausible", "below 80: implausible",
               "above 60: implausible")
)
cleaning_log
write.csv(cleaning_log, "cleaning_log.csv", row.names = FALSE)

# ---- Section 3: mean, median and standard deviation ----
income5 <- c(30, 35, 40, 45, 50)        # five incomes, $1000s per year
mean(income5)
median(income5)
income6 <- c(income5, 500)              # add one very high earner
mean(income6)
median(income6)

# Standard deviation, step by step
sbp <- c(110, 116, 120, 124, 130)
dev <- sbp - mean(sbp)                  # distance of each value from the mean
dev
dev^2                                   # squared distances
sum(dev^2) / (length(sbp) - 1)          # the variance
sqrt(sum(dev^2) / (length(sbp) - 1))    # the standard deviation
sd(sbp)                                 # R's built-in function gives the same

# ---- Section 3: annotated summary() and describe() output ----
summary(phaa$diastolic_bp)
library(psych)
describe(phaa[, c("diastolic_bp", "discrimination_score", "social_support_score")])

# ---- Section 3: odds and odds ratios from a table ----
# The outbreak table from the text
outbreak <- matrix(c(45, 15, 30, 60), nrow = 2,
                   dimnames = list(Exposure = c("Exposed", "Unexposed"),
                                   Group    = c("Case", "Control")))
outbreak
odds_cases    <- outbreak["Exposed", "Case"] / outbreak["Unexposed", "Case"]
odds_controls <- outbreak["Exposed", "Control"] / outbreak["Unexposed", "Control"]
odds_cases
odds_controls
odds_cases / odds_controls              # the odds ratio

# The same steps with the course data: hypertension by smoking status
tab <- table(Smoker = phaa$smoker, Hypertension = phaa$hypertension)
tab
odds_smokers     <- tab["Yes", "Yes"] / tab["Yes", "No"]
odds_non_smokers <- tab["No", "Yes"] / tab["No", "No"]
round(c(smokers = odds_smokers, non_smokers = odds_non_smokers,
        OR = odds_smokers / odds_non_smokers), 3)

# ---- Section 3: a Table 1 stratified by smoking status ----
# Three small helper functions. function(v) { ... } defines a reusable recipe
# that is applied to the whole sample and then to each smoking group.
grp <- phaa$smoker                       # the grouping variable ("No", "Yes")

mean_sd <- function(v) {                 # "mean (SD)" and the standardised difference
  f <- function(x) sprintf("%.1f (%.1f)", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
  m <- tapply(v, grp, mean, na.rm = TRUE)
  s <- tapply(v, grp, sd, na.rm = TRUE)
  d <- (m["Yes"] - m["No"]) / sqrt((s["Yes"]^2 + s["No"]^2) / 2)
  c(Overall = f(v), tapply(v, grp, f), SMD = sprintf("%.2f", d))
}
n_pct <- function(v, level) {            # "n (%)" and the standardised difference
  f <- function(x) sprintf("%d (%.1f%%)", sum(x == level, na.rm = TRUE),
                           100 * mean(x == level, na.rm = TRUE))
  p <- tapply(v == level, grp, mean, na.rm = TRUE)
  d <- (p["Yes"] - p["No"]) / sqrt((p["Yes"] * (1 - p["Yes"]) + p["No"] * (1 - p["No"])) / 2)
  c(Overall = f(v), tapply(v, grp, f), SMD = sprintf("%.2f", d))
}
n_miss <- function(v) {                  # number of missing values
  c(Overall = sum(is.na(v)), tapply(is.na(v), grp, sum), SMD = "")
}

table1 <- rbind(
  "N"                     = c(Overall = length(grp), table(grp), SMD = ""),
  "Age, mean (SD)"        = mean_sd(phaa$age),
  "  Age missing, n"      = n_miss(phaa$age),
  "Woman, n (%)"          = n_pct(phaa$gender, "Woman"),
  "Degree, n (%)"         = n_pct(phaa$education >= "Bachelor's", TRUE),
  "BMI, mean (SD)"        = mean_sd(phaa$bmi),
  "SBP, mean (SD)"        = mean_sd(phaa$systolic_bp),
  "Activity, mean (SD)"   = mean_sd(phaa$phys_act_min),
  "  Activity missing, n" = n_miss(phaa$phys_act_min),
  "Income missing, n (%)" = n_pct(is.na(phaa$income), TRUE)
)
noquote(table1)

# ---- Section 3: a Q-Q plot and the Shapiro-Wilk test ----
qqnorm(phaa$systolic_bp, main = "Systolic BP: normal Q-Q plot")
qqline(phaa$systolic_bp)          # the reference line
shapiro.test(phaa$systolic_bp)
shapiro.test(phaa$diastolic_bp)