The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Estimating District Means and SDs from Binned Test Scores, with the binest package

Paul T. von Hippel

2026-09-22

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

library(binest)
data(tx_g6_math_2018)
dim(tx_g6_math_2018)
#> [1] 1151    8
head(tx_g6_math_2018, 3)
#>   district_id district_name n_tested unsatisfactory approaches meets masters
#> 1           1    CAYUGA ISD       40              8         13    13       6
#> 2           2   ELKHART ISD       85             12         30    35       8
#> 3           3 FRANKSTON ISD       63              3         26    23      11
#>   reported_mean
#> 1          1657
#> 2          1653
#> 3          1675

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:

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.

t0 <- Sys.time()
fit_bm_known <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                           scope = "population")
#> fast_hetop: scope = "population". The reported SEs and 95% CIs reflect
#>   binning error only -- the coarsening from observing bin counts
#>   rather than individual scores; each group's units are treated as its
#>   whole population, so there is no sampling error,
#>   estimated from the model-implied Fisher information, which assumes
#>   within-group normality holds.
t_bm_known   <- as.numeric(Sys.time() - t0, units = "secs")
cor(fit_bm_known$est_raw$mean, truth)
#> [1] 0.9681113

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

fit_bm_sample <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                            scope = "sample")
#> fast_hetop: scope = "sample". The reported SEs and 95% CIs reflect
#>   sampling and binning error combined -- both drawing these units
#>   rather than others, and observing only bin counts rather than
#>   individual scores,
#>   estimated from the model-implied Fisher information, which assumes
#>   within-group normality holds.
head(cbind(population = fit_bm_known$est_raw$mean_se,
           sample     = fit_bm_sample$est_raw$mean_se))
#>      population   sample
#> [1,]   7.329076 21.19847
#> [2,]   4.090688 11.51668
#> [3,]   4.712282 12.82530
#> [4,]   8.680751 25.09820
#> [5,]   3.583806 10.12020
#> [6,]   7.301079 16.14469

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.

t0 <- Sys.time()
fit_bm_null <- fast_hetop(ngk, scope = "population")
#> fast_hetop: scope = "population". The reported SEs and 95% CIs reflect
#>   binning error only -- the coarsening from observing bin counts
#>   rather than individual scores; each group's units are treated as its
#>   whole population, so there is no sampling error,
#>   estimated from the model-implied Fisher information, which assumes
#>   within-group normality holds.
#>   Cutpoints were estimated from the pooled bin proportions; est_std is
#>   reported on the scale where the pooled distribution has mean 0 and SD 1 (the default
#>   standardized scale).
t_bm_null   <- as.numeric(Sys.time() - t0, units = "secs")

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.

fit_bm_eb <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cuts,
                        scope = "sample", estimator = "EB_shrunk")
#> fast_hetop: scope = "sample". The reported SEs and 95% CIs reflect
#>   sampling and binning error combined -- both drawing these units
#>   rather than others, and observing only bin counts rather than
#>   individual scores,
#>   estimated from the model-implied Fisher information, which assumes
#>   within-group normality holds.

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.

set.seed(1)
sub <- sample(nrow(ngk), 50)
t0 <- Sys.time()
fit_mle <- mle_hetop(ngk[sub, ], iterlim = 200)
#> Warning in mle_hetop(ngk[sub, ], iterlim = 200): optimization algorithm may not
#> have converged properly; see 'nlmdetails' element of object
t_mle <- as.numeric(Sys.time() - t0, units = "secs")
cor(fit_mle$est_star$mug, truth[sub])
#> [1] 0.9716946

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.

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.

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)	0.2 s	1151 / 1151
`fast_hetop()` (cutpoints from data)	0.3 s	1151 / 1151
`mle_hetop()` (50-district subsample)	2.3 s	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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.