# =============================================================================
# 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  (install once with install.packages("psych"))
# 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"))

# (answer key) keep the state before conversion so Q1 can compare before/after
age_before <- phaa$age
class(age_before)
summary(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)

# (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)

# 3. Range checks: replace 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

# 4. Order education and income so that comparisons (<= ...) are meaningful
phaa$education <- factor(phaa$education,
  levels = c("Less than high school", "High school", "Some college",
             "Bachelor's", "Graduate degree"),
  ordered = TRUE)

# ---- 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, plots, and a stratified Table 1 ====

# 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 -----------------
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")

# ---- 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)

# 2. How many underlying factors? ------------------------------------------
fa.parallel(dep_items)        # scree plot suggests one factor

# 3. Inspect a one-factor solution for the depression items ----------------
# NOTE (answer key): factanal() needs complete rows and dep4 has 15 missing
# values, so the page's original call factanal(x = dep_items, ...) stops with
# "'x' must contain finite values only". na.omit() drops those rows first, as
# step 4 already does for the combined items.
factanal(x = na.omit(dep_items), factors = 1, rotation = "varimax")
# 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 ------------------
phaa$dep_score <- rowSums(dep_items, na.rm = TRUE)
phaa$anx_score <- rowSums(anx_items, na.rm = TRUE)
summary(phaa$dep_score)

# 6. Save the cleaned, scale-augmented file --------------------------------
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, 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, rotation = "varimax")
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("\nSaved phaa_survey_clean.csv with", nrow(phaa), "rows and", ncol(phaa), "columns\n")