# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 6: Ecological and Group-Level Studies
# Answer key for the in-lesson R activities
# Data file(s): none  (the three-group 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 3: Inferential Errors & Sources of Ecologic Bias / Build the ecological fallacy from scratch ====

set.seed(230)

# Three groups; means line up positively, individuals are flat within group.
group <- rep(c("A", "B", "C"), each = 50)
x <- c(rnorm(50, 2), rnorm(50, 5), rnorm(50, 8))
y <- c(rnorm(50, 3), rnorm(50, 6), rnorm(50, 9))

# Pooled (individual-level) correlation -- dominated by group means
cor(x, y)

# Within each group (truth at the person level)
tapply(seq_along(x), group, function(i) cor(x[i], y[i]))

# Group-level (ecologic) correlation: nearly perfect by construction
gmean <- aggregate(cbind(x, y), list(group = group), mean)
cor(gmean$x, gmean$y)

# Stretch: visualise the discrepancy
plot(x, y, col = factor(group), pch = 19,
     xlab = "X", ylab = "Y",
     main = "Ecological fallacy: groups separate, individuals flat")
points(gmean$x, gmean$y, pch = 8, cex = 3, col = "black")

# ---- Activity questions (rActivity-230-6-1) ----
r_pooled <- cor(x, y)
r_within <- tapply(seq_along(x), group, function(i) cor(x[i], y[i]))
r_group  <- cor(gmean$x, gmean$y)
cat("Q1: within-group correlations:", paste(names(r_within), round(r_within, 3), collapse = ", "),
    "-> the individual-level relationship (flat); group-mean correlation =", round(r_group, 4),
    "-> the relationship between groups only.\n")

# Q2: decompose the variance of x and y into between-group and within-group parts
between_x <- var(gmean$x[match(group, gmean$group)]) ; within_x <- mean(tapply(x, group, var))
between_y <- var(gmean$y[match(group, gmean$group)]) ; within_y <- mean(tapply(y, group, var))
cat("Q2: pooled r =", round(r_pooled, 3), ". Between-group variance of x =", round(between_x, 2),
    "vs within-group", round(within_x, 2), "; for y", round(between_y, 2), "vs", round(within_y, 2),
    ". Between-group spread dominates, so the pooled r follows the diagonal of the group means",
    "rather than the flat within-cluster slopes.\n")
print(gmean)

cat("Q3: from the three stars alone one would infer a strong positive individual-level",
    "relationship (r =", round(r_group, 4), "); within every group the slope is ~0",
    "(r =", paste(round(r_within, 2), collapse = ", "), "), so that inference would be wrong:",
    "the ecological fallacy.\n")
within_slopes <- tapply(seq_along(x), group, function(i) unname(coef(lm(y[i] ~ x[i]))[2]))
cat("    Within-group regression slopes:", paste(names(within_slopes), round(within_slopes, 3), collapse = ", "),
    "; slope through the group means:", round(unname(coef(lm(y ~ x, data = gmean))[2]), 3), "\n")

# ==== Section 5: Final Assessment / Activity, Seeing the ecological fallacy in simulated data ====

set.seed(230)

# Three groups; means line up positively, individuals are flat within group.
group <- rep(c("A", "B", "C"), each = 50)
x <- c(rnorm(50, 2), rnorm(50, 5), rnorm(50, 8))
y <- c(rnorm(50, 3), rnorm(50, 6), rnorm(50, 9))

# Individual-level correlation (overall - mostly driven by group means)
cor(x, y)

# Individual-level correlation WITHIN each group (truth at the person level)
tapply(seq_along(x), group, function(i) cor(x[i], y[i]))

# Group-level (ecologic) correlation: nearly perfect by construction
gmean <- aggregate(cbind(x, y), list(group = group), mean)
cor(gmean$x, gmean$y)