# Definitive Table 1: exclusive dummies, recovery breakpoint 26.
# This specification reproduces every published coefficient except B1's gamma0,
# which appears in the manuscript with a flipped sign.
suppressPackageStartupMessages({library(readxl); library(sandwich); library(lmtest)})

df <- read_excel("reponse_all.xlsx")
df$Type <- gsub("ai_ls", "ai-ls", df$Type)
df$W <- df$WOA1
df$H <- as.integer(grepl("human", df$Type))
df$S <- as.integer(grepl("ls",    df$Type))
df$C1 <- as.integer(df$Row >= 19 & df$Row < 26)
df$C2 <- as.integer(df$Row >= 26)
df$pid <- rep(1:120, each = 30)

star <- function(p) if (p < .001) "***" else if (p < .01) "**" else if (p < .05) "*" else "n.s."

fits <- list(
  A1 = W ~ C1,
  A2 = W ~ H + C1 + H:C1,
  A3 = W ~ S + C1 + S:C1,
  B1 = W ~ C1 + C2,
  B2 = W ~ H + C1 + C2 + H:C1 + H:C2,
  B3 = W ~ S + C1 + C2 + S:C1 + S:C2
)

for (nm in names(fits)) {
  m  <- lm(fits[[nm]], data = df)
  hc <- coeftest(m, vcov. = vcovHC(m, type = "HC1"))
  cl <- coeftest(m, vcov. = vcovCL(m, cluster = df$pid))
  cat("\n#####", nm, " F =", sprintf("%.3f", summary(m)$fstatistic[1]),
      " R2 =", sprintf("%.4f", summary(m)$r.squared), "\n")
  for (i in seq_len(nrow(hc)))
    cat(sprintf("  %-14s b=%8.4f | robustSE=%6.4f p=%.4f %-4s | clusterSE=%6.4f p=%.4f %s\n",
                rownames(hc)[i], hc[i,1], hc[i,2], hc[i,4], star(hc[i,4]),
                cl[i,2], cl[i,4], star(cl[i,4])))
}