# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 7: Conceptualization, Measurement, and Causal Specification
# Answer key for the in-lesson R activities
# Data file(s): none  (the 2,000-person 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 1: Construct Validity & Measurement / Watch attenuation bias shrink a true effect toward zero ====

set.seed(230)
n <- 2000

# True exposure X, outcome Y with true slope = 1
X <- rnorm(n, mean = 10, sd = 2)
Y <- 2 + 1*X + rnorm(n, sd = 1)

# Noisy version of X (e.g., FFQ-measured dietary intake)
X_noisy <- X + rnorm(n, sd = 2)

# Compare slopes from "perfect" vs "noisy" exposure
coef(lm(Y ~ X))["X"]               # expect ~1.00 (truth)
coef(lm(Y ~ X_noisy))["X_noisy"]   # expect < 1.00 (attenuated)

# Stretch: how does noise size change the attenuation?
sds <- c(0, 0.5, 1, 2, 4)
sapply(sds, function(s) {
  Xn <- X + rnorm(n, sd = s)
  unname(coef(lm(Y ~ Xn))[2])
})

# ---- Activity questions (rActivity-230-7-1) ----
# NOTE (answer key): the stretch sapply() above consumes random numbers, so the
# slopes below are recomputed from the same seed to match the lesson's printed output.
set.seed(230)
X <- rnorm(n, mean = 10, sd = 2)
Y <- 2 + 1*X + rnorm(n, sd = 1)
X_noisy <- X + rnorm(n, sd = 2)
b_clean <- unname(coef(lm(Y ~ X))["X"])
b_noisy <- unname(coef(lm(Y ~ X_noisy))["X_noisy"])
slopes  <- sapply(sds, function(s) { Xn <- X + rnorm(n, sd = s); unname(coef(lm(Y ~ Xn))[2]) })
names(slopes) <- paste0("sd=", sds)

reliab <- var(X) / (var(X) + 2^2)
cat("Q1: clean slope =", round(b_clean, 3), ", noisy slope =", round(b_noisy, 3),
    "-> attenuated by", round(100 * (1 - b_noisy / b_clean)), "%.",
    "Reliability = var(X)/(var(X)+4) =", round(reliab, 3), "(var(X) =", round(var(X), 2), ").\n")

cat("Q2: slopes by noise SD:\n"); print(round(slopes, 3))
cat("    Expected attenuation factor var(X)/(var(X)+sd^2):",
    paste(round(var(X) / (var(X) + sds^2), 3), collapse = " "), "\n")
cat("    Cutting noise SD from 2 to 1 moves the slope from", round(slopes["sd=2"], 2), "to",
    round(slopes["sd=1"], 2), ": from about half of the truth to about four-fifths of it.\n")

cat("Q3: with FFQ-level noise (reliability ~0.2-0.4) a true slope of 1.0 is reported as",
    round(slopes["sd=4"], 2), "-", round(slopes["sd=2"], 2),
    "; apparent null findings can reflect noisy measurement rather than no effect.",
    "Larger n does not help: attenuation depends on the noise-to-signal ratio.\n")

# ==== Section 4: Final Assessment / Activity: Attenuation bias from noisy exposure measurement ====

set.seed(230)
n <- 2000

# True exposure X, outcome Y with true slope = 1
X <- rnorm(n, mean = 10, sd = 2)
Y <- 2 + 1*X + rnorm(n, sd = 1)

# Noisy version of X (e.g., FFQ-measured dietary intake)
X_noisy <- X + rnorm(n, sd = 2)

# Compare slopes from "perfect" vs "noisy" exposure
coef(lm(Y ~ X))["X"]             # expect ~1.00 (truth)
coef(lm(Y ~ X_noisy))["X_noisy"] # expect < 1.00 (attenuated)

## -----------------------------------------------------------------------------
## Stretch: how does noise size change the attenuation?
## -----------------------------------------------------------------------------
sds <- c(0, 0.5, 1, 2, 4)
sapply(sds, function(s) {
  Xn <- X + rnorm(n, sd = s)
  unname(coef(lm(Y ~ Xn))[2])
})
# Larger SD of measurement error -> larger attenuation toward zero.