## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE,
  fig.width = 7, fig.height = 4.3, dpi = 150, fig.align = "center"
)
set.seed(2025)

## ----libs---------------------------------------------------------------------
library(mixqrgate)
library(ggplot2)

pal <- c("Supportive" = "#1b6ca8", "Under-resourced" = "#e07b39")
theme_set(theme_minimal(base_size = 12))

## ----data---------------------------------------------------------------------
raw <- sim_gate2(n = 800, gamma = c(-0.3, 1.3),
                 b1 = c(48, 7), b2 = c(55, 2), sigma = c(6, 7))
schools <- data.frame(score = raw$y, ses = raw$x, funding = raw$z)
head(round(schools, 2))

## ----rawplot, fig.alt = "Scatter of score against SES, two latent regimes overlapping."----
# class 1 is, by construction, the steep-gradient (under-resourced) regime
schools$truth <- ifelse(raw$class == 1, "Under-resourced", "Supportive")
ggplot(schools, aes(ses, score, colour = truth)) +
  geom_point(alpha = 0.5, size = 1.3) +
  scale_colour_manual(values = pal, name = "True regime") +
  labs(x = "Student SES (standardised)", y = "Test score",
       title = "Two regimes are present but entangled") +
  theme(legend.position = "top")

## ----fit----------------------------------------------------------------------
fit <- mixqrgate(score ~ ses, data = schools, gating = ~ funding,
                 G = 2, tau = 0.5, variance = "louis")
fit

## ----summary------------------------------------------------------------------
summary(fit)

## ----oddsratio----------------------------------------------------------------
g <- coef(fit, "gating")["funding", 1, 1]
cat(sprintf("odds ratio per 1 SD of funding: %.2f\n", exp(g)))
confint(fit)        # gate-coefficient confidence intervals

## ----vcompare-----------------------------------------------------------------
se <- function(v) sqrt(diag(vcov(mixqrgate(score ~ ses, data = schools,
        gating = ~ funding, G = 2, tau = 0.5, variance = v,
        control = mixqrgate_control(seed = 1)))))
rbind(sandwich = se("sandwich"), louis = se("louis"))

## ----lvdata-------------------------------------------------------------------
lv <- sim_gate2(n = 800, gamma = c(-0.3, 1.3),
                b1 = c(48, 7), b2 = c(55, 2), sigma = c(6, 7), loc_vary = 3)
schools_lv <- data.frame(score = lv$y, ses = lv$x, funding = lv$z)

## ----vary---------------------------------------------------------------------
grid <- c(0.1, 0.25, 0.5, 0.75, 0.9)
fitv <- mixqrgate(score ~ ses, data = schools_lv, gating = ~ funding,
                  G = 2, tau = grid, vary_gating = "discrete",
                  variance = "louis")
round(fitv$gate_prob, 3)

## ----gatevtau, fig.alt = "Average gate probability against the quantile level with uncertainty bands."----
band <- do.call(rbind, lapply(seq_along(grid), function(g) {
  V <- fitv$gate_vcov[[g]]; gam <- as.numeric(fitv$gamma[, , g])
  L <- chol(V + 1e-8 * diag(nrow(V)))
  d <- replicate(500, {
    gd <- matrix(gam + as.numeric(crossprod(L, rnorm(length(gam)))),
                 length(fitv$znames))
    mean(mixqrgate:::gate_predict(gd, fitv$z)[, 2])
  })
  data.frame(tau = grid[g], p = mean(d),
             lo = quantile(d, .025), hi = quantile(d, .975))
}))

ggplot(band, aes(tau, p)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = pal[2], alpha = 0.2) +
  geom_line(colour = pal[2], linewidth = 1.1) +
  geom_point(colour = pal[2], size = 2.3) +
  ylim(0, 1) +
  labs(x = expression(tau), y = "P(under-resourced regime)",
       title = "The mix shifts across the score distribution",
       subtitle = "Per-quantile gate estimates, with simulated uncertainty")

## ----classify, fig.alt = "Posterior probability of the under-resourced regime across the data."----
schools$p_under <- fit$posterior[, 2, 1]
ggplot(schools, aes(ses, score, colour = p_under)) +
  geom_point(size = 1.4) +
  scale_colour_gradient(low = pal[1], high = pal[2],
                        name = "P(under-resourced)") +
  labs(x = "Student SES", y = "Test score",
       title = "Soft classification: confident at the edges, unsure in the middle")

## ----gatecurve, fig.alt = "Fitted probability of each regime as a function of school funding."----
nd <- data.frame(funding = seq(-2.5, 2.5, length.out = 100))
pp <- predict(fit, newdata = nd, type = "prob")
gd <- data.frame(funding = rep(nd$funding, 2),
                 prob = c(pp[, 1], pp[, 2]),
                 regime = rep(c("Supportive", "Under-resourced"), each = 100))
ggplot(gd, aes(funding, prob, colour = regime)) +
  geom_line(linewidth = 1.2) +
  scale_colour_manual(values = pal, name = NULL) +
  labs(x = "School funding (standardised)", y = "Gate probability",
       title = "The gate: funding shifts the regime mix") +
  theme(legend.position = "top")

## ----diag---------------------------------------------------------------------
c(converged = all(fit$converged),
  min_gate_prob = round(min(fit$gate_prob), 3),
  gate_condition = round(fit$gate_cond, 1))

## ----report-------------------------------------------------------------------
set.seed(2025)
fit_final <- mixqrgate(score ~ ses, data = schools, gating = ~ funding,
                       G = 2, tau = 0.5, variance = "louis",
                       control = mixqrgate_control(seed = 2025))
list(slopes = round(fit_final$beta["ses", , 1], 2),
     gate_funding = round(coef(fit_final, "gating")["funding", 1, 1], 3),
     gate_prob = round(fit_final$gate_prob[, 1], 3),
     se_method = fit_final$se_method)