# =====================================================================
# Talk to the Ficus (TTF) - analysis script
# FICTIONAL DEMO STUDY - created for an EHPS SIG Open Science webinar.
# No participants, plants or staplers were harmed.
#
# Runs in base R (>= 4.0). Optional sensitivity analysis needs 'mice'.
# Input: TTF_data_v2_clean.csv (same folder as this script)
# Follows the analysis plan in the preregistration (Analysis Plan, q17-q22).
# =====================================================================

d <- read.csv("TTF_data_v2_clean.csv", stringsAsFactors = FALSE)
d$group       <- factor(d$group, levels = c("stapler", "ficus"))  # stapler = reference
d$plant_named <- factor(d$plant_named, levels = c("no", "yes"))

# --- Sanity checks (preregistered: PSS-10 range 0-40) -------------------
stopifnot(all(d$pss_pre  >= 0 & d$pss_pre  <= 40))
stopifnot(all(d$pss_post >= 0 & d$pss_post <= 40, na.rm = TRUE))

# --- Descriptives ---------------------------------------------------------
cat("\nRandomised per group:\n"); print(table(d$group))
cat("\nWithdrawn (no post-test) per group:\n"); print(tapply(is.na(d$pss_post), d$group, sum))
cat("\nMeans (SD) by group:\n")
print(aggregate(cbind(pss_pre, pss_post, talk_minutes) ~ group, data = d,
                FUN = function(x) round(c(M = mean(x, na.rm = TRUE), SD = sd(x, na.rm = TRUE)), 2),
                na.action = na.pass))

# --- H1: ANCOVA, post-test PSS ~ group + baseline PSS (complete cases) ----
cc <- subset(d, !is.na(pss_post))
m1 <- lm(pss_post ~ group + pss_pre, data = cc)
cat("\n=== H1: ANCOVA ===\n"); print(summary(m1))
b1  <- coef(m1)["groupficus"]
ci1 <- confint(m1)["groupficus", ]
cat(sprintf("Adjusted difference (ficus - stapler): %.2f, 95%% CI [%.2f, %.2f]\n", b1, ci1[1], ci1[2]))
cat(sprintf("Standardised (b / SD of baseline PSS): d = %.2f\n", b1 / sd(d$pss_pre)))

# --- H2: non-inferiority of leaf loss (margin = 2 leaves) -----------------
d$leaf_change <- d$leaves_post - d$leaves_pre
m2   <- lm(leaf_change ~ group, data = d)
ci90 <- confint(m2, level = 0.90)["groupficus", ]
cat("\n=== H2: Leaf non-inferiority ===\n")
cat(sprintf("Difference in leaf change (ficus - stapler): %.2f, 90%% CI [%.2f, %.2f]\n",
            coef(m2)["groupficus"], ci90[1], ci90[2]))
cat(if (ci90[1] > -2) "Lower bound > -2: the ficus is NOT harmed by being talked to.\n"
    else "Lower bound <= -2: non-inferiority not shown. Apologise to the ficus.\n")

# --- Exploratory: does naming the plant matter? ---------------------------
m3 <- lm(pss_post ~ group * plant_named + pss_pre, data = cc)
cat("\n=== Exploratory: group x plant_named ===\n"); print(summary(m3)$coefficients)

# --- Sensitivity: multiple imputation (m = 20, PMM) -----------------------
if (requireNamespace("mice", quietly = TRUE)) {
  imp <- mice::mice(d[, c("group", "pss_pre", "pss_post", "age", "phd_year")],
                    m = 20, method = "pmm", seed = 2026, printFlag = FALSE)
  fit <- with(imp, lm(pss_post ~ group + pss_pre))
  cat("\n=== Sensitivity: multiple imputation ===\n"); print(summary(mice::pool(fit)))
} else {
  cat("\n(Package 'mice' not installed: sensitivity analysis skipped.)\n")
}