---
title: "Estimating District Means and SDs from Binned Test Scores, with the binest package"
author: "Paul T. von Hippel"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating District Means and SDs from Binned Test Scores, with the binest package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>"
)
```

# Overview

This vignette compares the package's three functions on real data:
district-level bin counts and reported mean scores from the 2017-18
State of Texas Assessments of Academic Readiness (STAAR) Grade 6
mathematics test. Because the state publishes each district's mean,
every estimate can be scored against a known truth, on accuracy and on
runtime.

All three functions fit the same heteroskedastic ordered probit
(HETOP) model (Reardon, Shear, Castellano & Ho 2017), in which each
district's scores are normally distributed around a mean and SD of
its own. They differ only in how
they fit it: one district at a time, jointly by maximum likelihood, or
by MCMC.

`fast_hetop()` is the recommended function; `mle_hetop()` and
`fh_hetop()` are deprecated and retained only for comparison. See
`?fast_hetop` for arguments and return value, and the accompanying
paper for the fuller empirical case.

# The Texas data

```{r}
library(binest)
data(tx_g6_math_2018)
dim(tx_g6_math_2018)
head(tx_g6_math_2018, 3)
```

The dataset has 1,151 districts, each with bin counts in four
proficiency categories and a reported average score that
serves as ground truth.

The published cut scores defining the bin boundaries are 1536,
1653, and 1772:

```{r}
ngk  <- with(tx_g6_math_2018,
             cbind(unsatisfactory, approaches, meets, masters))
cuts <- c(1536, 1653, 1772)
truth <- tx_g6_math_2018$reported_mean
```

# Method 1: `fast_hetop()` with known cutpoints

This uses the published cut scores directly. Output `est_raw$mean` is on the test-score scale; `est_std$mean`
is on a standardized scale (population-weighted state mean 0, total
state SD 1).

`scope` is required and has no default. It says whether the students
behind each district's counts are a sample from some larger population
or that district's entire population, which decides what the reported
standard errors mean. Texas reports counts for every tested student, so
a district's true mean is the actual mean of those students' scores —
there is no sampling error, and the only uncertainty is that we see four
bin counts instead of individual scores. That is `scope = "population"`,
and the returned SEs reflect the binning alone. Data from a survey or a
simulation would instead be `scope = "sample"`, whose SEs add the
variability from having drawn these units rather than others.

```{r}
t0 <- Sys.time()
fit_bm_known <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                           scope = "population")
t_bm_known   <- as.numeric(Sys.time() - t0, units = "secs")
cor(fit_bm_known$est_raw$mean, truth)
```

The two choices give the same point estimates and differ only in the
reported uncertainty:

```{r}
fit_bm_sample <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                            scope = "sample")
head(cbind(population = fit_bm_known$est_raw$mean_se,
           sample     = fit_bm_sample$est_raw$mean_se))
```

```{r, fig.width = 4.5, fig.height = 4.5}
plot(truth, fit_bm_known$est_raw$mean,
     pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "True district mean",
     ylab = "Estimated mean (test-score scale)",
     main = "fast_hetop (known cuts)")
abline(0, 1, col = "red", lty = 2)
```

# Method 2: `fast_hetop()` with cutpoints estimated from data

When cutpoints are not known, `fast_hetop()` derives them from pooled
state bin proportions via `qnorm(cumsum(pooled_props))`. Output is on
the standardized scale.

```{r}
t0 <- Sys.time()
fit_bm_null <- fast_hetop(ngk, scope = "population")
t_bm_null   <- as.numeric(Sys.time() - t0, units = "secs")
```

```{r, fig.width = 4.5, fig.height = 4.5}
plot(truth, fit_bm_null$est_std$mean,
     pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "True district mean",
     ylab = "Estimated mean (standardized)",
     main = "fast_hetop (cuts from data)")
```

# Method 3: `fast_hetop()` with empirical-Bayes shrinkage

The `estimator = "EB_shrunk"` option applies normal-normal empirical-Bayes
shrinkage to the per-district means and log-SDs, using a moment estimator
of the between-district variance and the Fisher-information sampling
variance of each district's MLE. The shrunken estimates come back in
the same `mean` and `sd` fields as the ML estimates, with
`estimator` recording which was used.

Shrinkage requires `scope = "sample"`, and `fast_hetop()` refuses it
under `scope = "population"`. The Texas data cover the whole
population, so the call below declares `scope = "sample"` purely to
demonstrate the option; for a real population analysis, use the
default `"ML"` estimator.

```{r}
fit_bm_eb <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                        scope = "sample", estimator = "EB_shrunk")
```

```{r, fig.width = 4.5, fig.height = 4.5}
plot(truth, fit_bm_eb$est_raw$mean,
     pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "True district mean",
     ylab = "EB-shrunk estimated mean",
     main = "fast_hetop (EB_shrunk, known cuts)")
```

# Method 4: `mle_hetop()` on a 50-district subsample

`mle_hetop()` maximizes the likelihood over all districts at once,
rather than one district at a time. On the full 1,151-district
dataset it does not converge within a reasonable runtime, so we
illustrate it on a representative random subsample of 50 districts.

```{r}
set.seed(1)
sub <- sample(nrow(ngk), 50)
t0 <- Sys.time()
fit_mle <- mle_hetop(ngk[sub, ], iterlim = 200)
t_mle <- as.numeric(Sys.time() - t0, units = "secs")
cor(fit_mle$est_star$mug, truth[sub])
```

```{r, fig.width = 4.5, fig.height = 4.5}
plot(truth[sub], fit_mle$est_star$mug,
     pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "True district mean",
     ylab = "Estimated mean (standardized)",
     main = "HETOP MLE (50-district subsample)")
```

# Method 5: `fh_hetop()` on the full dataset

`fh_hetop()` fits the model by MCMC (Lockwood, Castellano & Shear
2018), placing a prior over the district parameters and reporting
posterior means. Each Gibbs
iteration is linear in the number of groups; on 1,151 districts the
model fits in about nine minutes. The runtime in this vignette is deliberately short to keep
the package build tractable; for production use, longer chains are
recommended.

```{r, eval = FALSE}
t0 <- Sys.time()
fit_fh <- fh_hetop(
  ngk       = ngk,
  p         = c(10, 10),
  m         = c(100, 100),
  gridL     = c(-5.0, log(0.10)),
  gridU     = c( 5.0, log(5.0)),
  n.iter    = 2000,
  n.burnin  = 1000,
  seed      = 3142
)
t_fh <- as.numeric(Sys.time() - t0, units = "secs")
cor(fit_fh$fh_hetop_extras$est_star_mug$theta_pm, truth)
```

Because this chunk takes about nine minutes, the package ships with
the per-district HETOP-Bayes posterior means pre-computed and stored
in `inst/extdata/fh_hetop_means_tx_g6_math_2018.rds`. The scatterplot
below uses those cached values.

```{r, fig.width = 4.5, fig.height = 4.5}
fh_means_file <- system.file("extdata",
                             "fh_hetop_means_tx_g6_math_2018.rds",
                             package = "binest")
fh_means <- readRDS(fh_means_file)

plot(truth, fh_means,
     pch = 16, cex = 0.5, col = rgb(0, 0, 0, 0.3),
     xlab = "True district mean",
     ylab = "Estimated mean (standardized)",
     main = "HETOP Bayes (full data, cached)")
```

# Summary

| Method | Runtime | Districts |
|--------|---------|-----------|
| `fast_hetop()` (known cutpoints) | `r sprintf("%.1f s", t_bm_known)` | `r sum(!is.na(fit_bm_known$est_raw$mean))` / 1151 |
| `fast_hetop()` (cutpoints from data) | `r sprintf("%.1f s", t_bm_null)` | `r sum(!is.na(fit_bm_null$est_std$mean))` / 1151 |
| `mle_hetop()` (50-district subsample) | `r sprintf("%.1f s", t_mle)` | 50 / 50 |
| fh_hetop (full data, not run above) | 450 s | 1151 / 1151 |

Note: `fast_hetop()` ran on all 1,151 districts in a fraction of a second --- roughly 1,000 times faster than `fh_hetop()`, and roughly 10 times faster than `mle_hetop()`, even though `mle_hetop()` ran on a 50-district subsample rather than the full data.

# References

Lockwood, J. R., Castellano, K. E., & Shear, B. R. (2018). Flexible Bayesian models for inferences from coarsened, group-level achievement data. *Journal of Educational and Behavioral Statistics, 43*(6), 663–692. https://doi.org/10.3102/1076998618795124

Reardon, S. F., Shear, B. R., Castellano, K. E., & Ho, A. D. (2017). Using heteroskedastic ordered probit models to recover moments of continuous test score distributions from coarsened data. *Journal of Educational and Behavioral Statistics, 42*(1), 3–45. https://doi.org/10.3102/1076998616666279