--- title: "A Primer on Location-Varying Gating for Quantile Mixtures" author: "Kailas Venkitasubramanian, University of North Carolina at Charlotte" output: rmarkdown::html_vignette: toc: true toc_depth: 2 bibliography: mixqrgate.bib vignette: > %\VignetteIndexEntry{A Primer on Location-Varying Gating for Quantile Mixtures} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r 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) ``` This primer is for applied researchers (in education, economics, sociology, public health) who suspect their population is a **mixture of latent groups**, who care about the **whole distribution** of an outcome rather than its average, and who want to know whether **group membership itself changes across that distribution**. No background in mixture models or quantile regression is assumed. By the end you will be able to fit the model, read every number it returns, judge whether to trust it, and report it responsibly. ```{r libs} library(mixqrgate) library(ggplot2) pal <- c("Supportive" = "#1b6ca8", "Under-resourced" = "#e07b39") theme_set(theme_minimal(base_size = 12)) ``` # 1. The question Imagine test scores from many schools. Two stories could be true at once: * **One population, unequal returns.** Family socio-economic status (SES) raises scores, more steeply at the bottom of the distribution than at the top. * **Two kinds of school.** Some schools are *under-resourced*: scores track SES steeply, so advantaged children pull far ahead. Others are *supportive*: the SES gradient is flatter and the baseline higher, so the school compensates. If the second story holds, two further questions follow. **Which schools are which, and what predicts it?** And -- the question this package is built for -- **does the mix change across the score distribution?** A school might look supportive among its low scorers and under-resourced among its high scorers. That is *location-varying mixing* [@furno2025], and until now it had no inferential tool in R. # 2. The model in one page We model the conditional **quantile** of the outcome, not its mean. The \eqn{\tau}-th quantile (\eqn{\tau = 0.5} is the median) of score given SES is a mixture of \eqn{G} latent regimes [@wuyao2016]: * In regime \eqn{k}, the \eqn{\tau}-th quantile is linear in SES, \eqn{x_i'\beta_k(\tau)}, estimated by minimising the **check loss** [@koenkerbassett1978], the workhorse of quantile regression [@koenker2005]. * A school belongs to regime \eqn{k} with probability \eqn{\pi_{ik}}. In an ordinary mixture these mixing weights are constants. Here they follow a **multinomial-logit gate** on *gating covariates* \eqn{z} (e.g. school funding): \deqn{\pi_{ik} = \frac{\exp(z_i'\gamma_k)}{\sum_l \exp(z_i'\gamma_l)}.} This is the concomitant-variable idea of @grunleisch2008, here attached to *quantile* components. * The new move: the gate coefficients may carry a \eqn{\tau} argument, \eqn{\gamma_k(\tau)}, so the mix can shift across the distribution. Estimation is by EM [@mclachlanpeel2000]. The gate step is a weighted multinomial logistic regression of the soft group labels on \eqn{z}, the principled replacement for Furno's reweighting heuristic, and what lets us put **standard errors on the gate** (Section 6 discusses making those errors trustworthy). Setting `gating = ~1` removes the gate and recovers an ordinary mixture of quantile regressions. # 3. An illustrative dataset `sim_gate2()` generates this structure. We rename its columns to the schools story: `score`, `ses`, and `funding` (the gating covariate). Two regimes differ in their SES gradient; funding raises the odds of the supportive regime. ```{r 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)) ``` The two regimes overlap heavily in the raw scatter, which is the whole problem. You cannot eyeball the groups; you need a model to recover them. ```{r 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") ``` # 4. Fitting your first model One call. We fit two regimes at the median, gate on `funding`, and ask for the Louis standard errors (Section 6 explains the choice). ```{r fit} fit <- mixqrgate(score ~ ses, data = schools, gating = ~ funding, G = 2, tau = 0.5, variance = "louis") fit ``` `print` shows the class-average gate probabilities and the component coefficients. Components are **ordered by their SES slope**, so component 1 is always the flatter-gradient *supportive* regime and component 2 the steeper *under-resourced* one. The estimated slopes recover the two true gradients: the flatter near 2, the steeper near 6--7. # 5. Interpreting the gate `summary()` adds the piece Furno's method cannot: the gate coefficients with standard errors. ```{r summary} summary(fit) ``` Read the gate block as a logistic regression for *being in component 2 (under-resourced) rather than component 1 (supportive)*. The coefficient on `funding` is negative and significant: **better-funded schools are less likely to be in the steep-gradient regime**. Exponentiate it for an odds ratio. ```{r 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 ``` The class-average gate probabilities (`fit$gate_prob`) summarise the mix; the per-school probabilities live in `predict(fit, type = "prob")`. But before we trust that coefficient, the standard error behind it has to be sound, which in a mixture is less obvious than it looks. # 6. Getting the uncertainty right A mixture never knows for certain which group an observation belongs to. A naive standard error that ignores this -- the default `"sandwich"`, computed *as if* memberships were known -- is too small and under-covers. `mixqrgate` offers two classification-aware alternatives: * `variance = "louis"` -- an analytic observed-information correction [@louis1982] that subtracts the "missing information" from the latent labels. Fast, and the recommended default for inference. * `variance = "stochEM"` -- a multiple-imputation estimate [@littlerubin2002]: draw labels, refit the gate, combine. Slower, a useful cross-check. The difference is not cosmetic. The conditional sandwich under-covers (in a Monte-Carlo check it held a nominal 95% interval only about 80% of the time), and the Louis correction is designed to give back exactly the coverage the sandwich loses (in the same check it returned the interval to roughly nominal). ```{r 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")) ``` The Louis standard errors are larger: the real cost of not knowing the labels. # 7. Does membership shift across the distribution? Now the headline question. We need data in which the answer is *yes*: where a school's regime depends on where its students sit in the score distribution. `sim_gate2(loc_vary = 3)` builds that: membership is tied to the quantile rank, so the under-resourced regime is scarce among low scorers and common higher up. (With `loc_vary = 0`, the package default, the mix is constant across \eqn{\tau}; the tool should -- and does -- report no shift there.) ```{r 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) ``` Pass a grid of quantiles with `vary_gating = "discrete"`; the gate is refit at each. **Each gate carries its own uncertainty**, so "does the mix change with the quantile?" is answered with inference, not by eyeballing a wiggly line. ```{r 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) ``` We need uncertainty on each point. At every \eqn{\tau} we draw many plausible gates from that fit's covariance and record the spread of the implied class share, so the ribbon shows sampling noise, not just the point estimate. ```{r 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") ``` The under-resourced regime is least common in the lower tail and more common through the middle and upper range: the location-varying mixing built into the data, recovered by the model. Two cautions keep the reading sober. First, the per-quantile gates are fit independently and are noisy (the "classification ambiguity across \eqn{\tau}" of @wuyao2016), so neighbouring points wobble; lean on the bands, not the dots. Second, on data with a *constant* true gate (`loc_vary = 0`) this same picture comes out flat within its uncertainty; the method does not manufacture a trend. The value of the tool is that it lets you *test* the question with standard errors rather than over-reading a curve; `summary(fitv)` reports the gate coefficients with their SEs at each \eqn{\tau} for that. # 8. Reading the fit in pictures Two pictures earn their place here. (The regimes themselves are just the two fitted median gradients -- slopes near 2 and 6--7, already printed in Section 4 -- so we do not re-draw the scatter.) **Who is in which group?** Posterior probabilities express the (ir)reducible uncertainty of assignment. ```{r 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") ``` **The gate as a curve.** Because the gate is a logistic model on funding, we can draw the fitted regime probability against funding directly. ```{r 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") ``` # 9. Diagnostics: can you trust it? A few checks before believing any mixture fit. * **Convergence.** `fit$converged` should be `TRUE` at every \eqn{\tau}. `mixqrgate` warns otherwise. * **No collapsed component.** If a regime's average gate probability is tiny, the mixture has effectively one component; the package warns. Inspect `fit$gate_prob`. * **A well-conditioned gate.** Highly collinear gating covariates make the gate coefficients unstable; `fit$gate_cond` reports the conditioning and the package warns when it is large. * **Distinct components.** The regimes must differ (Theorem 2.1 of @wuyao2016); near-identical slopes are weakly identified. ```{r diag} c(converged = all(fit$converged), min_gate_prob = round(min(fit$gate_prob), 3), gate_condition = round(fit$gate_cond, 1)) ``` Multi-start (the `nstart` control) guards against local optima; the mixture likelihood is multimodal, so several starts are essential. # 10. Practical guidance and pitfalls * **Start at the median, gate on a few substantively chosen covariates.** A gate with many collinear predictors is hard to identify. * **Report Louis (or stochEM) standard errors,** never the conditional sandwich, whenever you make an inferential claim about the gate. * **Treat the discrete \eqn{\tau}-gate as exploratory.** Per-quantile gates are noisy; look at the uncertainty bands and resist over-reading a trend. * **Choose \eqn{G} for interpretability, not fit alone.** Two or three substantively meaningful regimes usually beat a larger number chosen by an information criterion. * **The component labels are slope-ordered,** so they are stable across runs; describe them by what their slopes mean. # 11. How mixqrgate relates to other tools * **`flexmix`** has the concomitant multinomial gate but no quantile / check-loss components, so it cannot produce \eqn{\tau}-quantile regimes [@grunleisch2008]. * **`brms`** can gate on covariates and fix a quantile per component, but cannot make the gate a function of \eqn{\tau} or fit several \eqn{\tau} jointly. * **Furno's (2025) Stata routine** is the only prior location-varying-gating tool; it is a reweighting heuristic with no standard errors on the gate. `mixqrgate` turns that idea into a likelihood/EM object with inference [@furno2025]. * **`mixqr`** is the parent: an ordinary mixture of quantile regressions. `mixqrgate` with `gating = ~1` reproduces it (to EM tolerance). # 12. Reporting and reproducibility A compact, reproducible report of the headline fit: ```{r 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) ``` Report: the number of regimes and why; the component gradients; the gate coefficients with their (classification-aware) standard errors and the odds-ratio interpretation; the standard-error method; and, if you used the \eqn{\tau}-grid, the gate-vs-\eqn{\tau} picture *with* its uncertainty. That last clause is the whole point. The schools were never labelled; the model recovered two gradients, told us funding predicts which a school follows, and -- when the data warranted it -- showed the mix shifting across the score distribution. Each of those is a claim with a standard error behind it. The discipline the primer asks for is to read those pictures with their uncertainty, and to test rather than eyeball. # References