# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 3: Introduction to Observational Studies
# Answer key for the in-lesson R activities
# Data file(s): none  (the 200-person cross-sectional dataset is simulated with set.seed(230))
# Packages: none beyond base R
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 4: Limitations, Incidence Estimation & Reporting / Build a 2x2 contingency table from raw cross-sectional data ====

# Simulate 200 people: half exposed, with higher disease prevalence among exposed.
set.seed(230)
n <- 200
exposure <- rep(c("exposed", "unexposed"), each = n/2)
disease  <- c(rbinom(100, 1, 0.30), rbinom(100, 1, 0.10))

dat <- data.frame(exposure = exposure, disease = disease)

# 2x2 cross-tabulation -- the workhorse of descriptive epi.
tab <- table(dat$exposure, dat$disease,
             dnn = c("Exposure", "Disease"))
tab

# Prevalence (proportion with disease=1) within each exposure group.
prop.table(tab, margin = 1)   # margin=1 -> divide each row by its row total

# ---- Activity questions (rActivity-230-3-1) ----
a <- tab["exposed", "1"];   b <- tab["exposed", "0"]
c <- tab["unexposed", "1"]; d <- tab["unexposed", "0"]
cat("Q1: a =", a, "(exposed, diseased), b =", b, "(exposed, not diseased),",
    "c =", c, "(unexposed, diseased), d =", d, "(unexposed, not diseased);",
    "totals:", a + b, "exposed,", c + d, "unexposed,", a + c, "diseased,", b + d, "not, n =", n, "\n")
addmargins(tab)

prev <- prop.table(tab, margin = 1)[, "1"]
cat("Q2: prevalence =", round(prev["exposed"], 2), "in the exposed and",
    round(prev["unexposed"], 2), "in the unexposed. Each is existing cases / people",
    "at one snapshot with no follow-up time, so it is prevalence, not incidence.\n")

# Q3: other seeds and a larger sample (400 per group) shrink the sampling error;
#     SE of a proportion = sqrt(p(1-p)/n), so 4x the sample halves it.
se100 <- sqrt(0.30 * 0.70 / 100); se400 <- sqrt(0.30 * 0.70 / 400)
cat("Q3: with this seed the exposed prevalence is", round(prev["exposed"], 2),
    "against a generating value of 0.30 (SE with 100 per group =", round(se100, 3),
    "). With 400 per group the SE falls to", round(se400, 3), "(half), so estimates",
    "settle closer to 0.30 and 0.10 and the contrast is more reliable; bias is unaffected.\n")
seeds <- c(1, 2, 3, 4, 5)
cat("    Exposed prevalence under five other seeds:",
    sapply(seeds, function(s) { set.seed(s); mean(rbinom(100, 1, 0.30)) }), "\n")
set.seed(230)
big <- c(mean(rbinom(400, 1, 0.30)), mean(rbinom(400, 1, 0.10)))
cat("    With 400 per group (seed 230): exposed", big[1], ", unexposed", big[2], "\n")

# ==== Section 4: Limitations, Incidence Estimation & Reporting / Activity, Building a 2x2 table and computing prevalence ====

# Simulate 200 people: half exposed, with higher disease prevalence among exposed
set.seed(230)
n        <- 200
exposure <- rep(c("exposed", "unexposed"), each = n/2)
disease  <- c(rbinom(100, 1, 0.30), rbinom(100, 1, 0.10))

dat <- data.frame(exposure = exposure, disease = disease)

# 2x2 cross-tabulation -- the workhorse of descriptive epi
tab <- table(dat$exposure, dat$disease,
             dnn = c("Exposure", "Disease"))
tab

# Prevalence (proportion with disease=1) within each exposure group
prop.table(tab, margin = 1)   # margin=1 -> divide each row by its row total

## -----------------------------------------------------------------------------
## Stretch: prevalence ratio (a quick preview of measures-of-association)
## -----------------------------------------------------------------------------
prev_exposed   <- prop.table(tab, margin = 1)["exposed",   "1"]
prev_unexposed <- prop.table(tab, margin = 1)["unexposed", "1"]
prev_ratio     <- prev_exposed / prev_unexposed
cat("Prevalence ratio:", round(prev_ratio, 2), "\n")
# NOTE (answer key): the generating prevalences are 0.30 and 0.10 (true PR = 3.0);
# with seed 230 the sample gives 0.22 / 0.05 = 4.4, a reminder that a 200-person
# sample estimates a ratio with considerable sampling error.