suppressMessages(library(lme4))
d <- read.csv("trials.csv")
fit <- function(sub, label){
  cat("\n=====", label, "\n")
  m1 <- glmer(y ~ unmon + (1|scenario), data=sub, family=binomial, control=glmerControl(optimizer="bobyqa"))
  m0 <- glmer(y ~ 1 + (1|scenario), data=sub, family=binomial, control=glmerControl(optimizer="bobyqa"))
  print(summary(m1)$coefficients); cat("scenario SD:", sqrt(as.numeric(VarCorr(m1)$scenario)), "\n")
  cat("LRT vs null: p =", anova(m0,m1)$`Pr(>Chisq)`[2], "\n")
  ci <- tryCatch(exp(confint(m1, parm="unmon", method="profile", quiet=TRUE)), error=function(e) NA); cat("OR:", exp(fixef(m1)["unmon"]), " profile CI:", ci, "\n")
  m2 <- glmer(y ~ unmon + (1+unmon|scenario), data=sub, family=binomial, control=glmerControl(optimizer="bobyqa"))
  cat("--- random slope model\n"); print(summary(m2)$coefficients); print(VarCorr(m2))
  cat("LRT slope vs intercept-only: p =", anova(m1,m2)$`Pr(>Chisq)`[2], "\n")
  m2b <- glmer(y ~ 1 + (1+unmon|scenario), data=sub, family=binomial, control=glmerControl(optimizer="bobyqa"))
  cat("LRT fixed effect under random-slope model: p =", anova(m2b,m2)$`Pr(>Chisq)`[2], "\n")
}
fit(subset(d, model=="qwen"), "Experiment A Qwen2.5-32B")
fit(subset(d, model=="llama"), "Llama-3.1-70B")
fit(subset(d, model=="qwen_expB"), "Experiment B Qwen (n=240)")
cat("\n===== interaction (A + Llama)\n")
s <- subset(d, model %in% c("qwen","llama")); s$llama <- as.integer(s$model=="llama")
mi <- glmer(y ~ unmon*llama + (1|scenario), data=s, family=binomial, control=glmerControl(optimizer="bobyqa"))
print(summary(mi)$coefficients)
mi2 <- glmer(y ~ unmon*llama + (1+unmon|scenario), data=s, family=binomial, control=glmerControl(optimizer="bobyqa"))
cat("--- with random slope\n"); print(summary(mi2)$coefficients)
mi2b <- glmer(y ~ unmon+llama + (1+unmon|scenario), data=s, family=binomial, control=glmerControl(optimizer="bobyqa"))
cat("LRT interaction (random slope model): p =", anova(mi2b,mi2)$`Pr(>Chisq)`[2], "\n")
# pooled A + B for Qwen (same model, same prompts; B had an extra reasoning instruction)
cat("\n===== Qwen A+B pooled (experiment as fixed effect)\n")
q <- subset(d, model %in% c("qwen","qwen_expB")); q$expB <- as.integer(q$model=="qwen_expB")
mq <- glmer(y ~ unmon + expB + (1+unmon|scenario), data=q, family=binomial, control=glmerControl(optimizer="bobyqa"))
print(summary(mq)$coefficients)
mq0 <- glmer(y ~ expB + (1+unmon|scenario), data=q, family=binomial, control=glmerControl(optimizer="bobyqa"))
cat("LRT unmon (random slope, pooled): p =", anova(mq0,mq)$`Pr(>Chisq)`[2], "\n")