Package {binest}

Estimation of Group Means and SDs from Binned Count Data

Description

Estimates group-level means and standard deviations from binned (coarsened) count data, where the within-bin scores are unobserved. All three functions fit the same heteroskedastic ordered probit (HETOP) model of Reardon, Shear, Castellano and Ho (2017), in which each group's values are normally distributed around a mean and SD of its own. They differ only in how they fit it, and they share a common output structure:

• fast_hetop(): fits each group separately, solving closed-form truncated-normal score equations. Runs in time linear in the number of groups times the number of bins. This is the preferred function in the package.

• mle_hetop(): maximizes the likelihood over all groups at once. Returns the same estimates as fast_hetop(estimator = "ML"), far more slowly. Deprecated in favor of fast_hetop().

• fh_hetop(): fits the model by MCMC, placing a Fay-Herriot hyperprior over the group parameters and reporting posterior means (Lockwood, Castellano and Shear 2018). Deprecated in favor of fast_hetop().

The mle_hetop() and fh_hetop() functions are forked from the HETOP package by J. R. Lockwood (CRAN, last released 2019). mle_hetop() has been modified to speed up its runtime via a vectorized inner loop and to remove two user-facing arguments (fixedcuts and svals) that some users found confusing; cutpoints and starting values are now derived internally from the data. fh_hetop() has been modified the same way, removing its fixedcuts argument for the same reason: it was never a channel for supplying cutpoints known in advance on the test-score scale, only a coordinate-system normalization, and supplying native-scale cutpoints could send the MCMC sampler into a degenerate region rather than producing an error.

mle_hetop() and fh_hetop() are deprecated in favor of fast_hetop() and remain in the package only for comparison purposes. See vignette("binest") for an empirical comparison on Texas STAAR Grade 6 mathematics data.

Bundled data

The package ships with tx_g6_math_2018, a district-level dataset of bin counts and reported mean scores from the 2017-18 administration of the State of Texas Assessments of Academic Readiness (STAAR) Grade 6 mathematics test. See the vignette for usage.

Author(s)

Paul T. von Hippel ph3828@eid.utexas.edu, David J. Hunter, and J. R. Lockwood.

References

Fisher, R. A. (1922). On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London A, 222, 309-368.

Lockwood, J. R., Castellano, K. E., and 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.

Reardon, S. F., Shear, B. R., Castellano, K. E., and 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.

Sheppard, W. F. (1898). On the calculation of the most probable values of frequency-constants for data arranged according to equidistant divisions of a scale. Proceedings of the London Mathematical Society, 29, 353-380.

Fast Estimation of Group Means and SDs from Binned Counts

Description

Education agencies often report school or district score distributions as the number of students scoring in each of several score ranges, or bins, separated by threshold scores, or cuts. The functions in this package translate those bin counts into estimates of the mean and standard deviation (SD). They do so using the heteroskedastic ordered probit (HETOP) model, which assumes that scores follow a normal distribution within each school or district, each of which has its own mean and SD.

fast_hetop() is the preferred function in this package. It can produce estimates for hundreds or thousands of schools or districts in less than a second. It can provide either maximum likelihood estimates (estimator = "ML") or shrunken empirical Bayes estimates (estimator = "EB_shrunk"). It also has several features that the other package functions lack, including

• a test for whether scores really follow a normal distribution, as the HETOP model assumes;

• the ability to accept published cut scores;

• the ability to adjust standard errors downward when the data represent a population rather than a sample.

fast_hetop(estimator = "ML") provides identical estimates to mle_hetop(), but runs much more quickly, especially when the number of schools or districts is large.

fast_hetop(estimator = "EB_shrunk") provides very similar estimates to fh_hetop(), and again runs much more quickly.

mle_hetop() and fh_hetop() were written by J. R. Lockwood; mle_hetop() has since been edited slightly by Paul T. von Hippel. Both are included in the package for comparison, but because of their relative slowness they are no longer maintained or recommended.

Usage

fast_hetop(ngk, cutpoints_known = FALSE, cutpoints = NULL,
          pooled_mean = 0, pooled_sd = 1,
          scope,
          estimator = "ML",
          tol = 1e-4, maxit = 100,
          conf.level = 0.95,
          estimate_unidentified_districts = TRUE)

Arguments

Value

A list with the following components:

Author(s)

Paul T. von Hippel and David J. Hunter.

References

Reardon S., Shear B.R., Castellano K.E. and 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.

Lockwood J.R., Castellano K.E. and 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.

Examples

set.seed(1001)
G   <- 10

## Let means and SDs vary across the groups.
mug       <- seq(from = -2.0, to = 2.0, length = G)
sigmag    <- seq(from =  2.0, to = 0.8, length = G)
cutpoints <- c(-1.0, 0.0, 0.8)
ng        <- rep(1000, G)
ngk       <- gendata_hetop(G, K = 4, ng, mug, sigmag, cutpoints)

## Here the counts were simulated by sampling ng scores per group, so
## scope = "sample" and the reported SEs will include sampling as well
## as binning error.
##
## Cutpoints known: both est_raw (test-score scale) and est_std
## (standardized scale) are returned.
bm <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cutpoints,
                 scope = "sample")
print(cbind(true = mug,       est = bm$est_raw$mean))
print(cbind(true = sigmag,    est = bm$est_raw$sd))
print(cbind(true = cutpoints, est = bm$est_raw$cutpoints))
print(cbind(est = bm$est_raw$mean,
            se  = bm$est_raw$mean_se,
            lo  = bm$est_raw$mean_ci_lower,
            hi  = bm$est_raw$mean_ci_upper))

## If the data represented the whole population, then scope =
## "population", and the reported standard errors are smaller because
## they reflect only binning error.
bm_pop <- fast_hetop(ngk, cutpoints_known = TRUE, cutpoints = cutpoints,
                     scope = "population")
print(cbind(sample_se     = bm$est_raw$mean_se,
            population_se = bm_pop$est_raw$mean_se))

## If cutpoints_known = FALSE, then cutpoints are unknown and are
## estimated on a scale with mean and SD given by pooled_mean
## (default 0) and pooled_sd (default 1).
bm2 <- fast_hetop(ngk, scope = "sample")
print(bm2$est_std$cutpoints)
print(bm2$est_std$mean)

## If the cut scores are unknown but the pooled mean and SD are known,
## you can pass the pooled mean and SD and the estimates will come back
## on that scale.
bm3 <- fast_hetop(ngk, pooled_mean = 1640.2, pooled_sd = 138.6,
                  scope = "sample")
print(bm3$est_std$cutpoints)

Fit Fay-Herriot Heteroskedastic Ordered Probit (FH-HETOP) Model using JAGS

Description

Fits the FH-HETOP model described by Lockwood, Castellano and Shear (2018) using the jags function in the suggested package R2jags. Requires JAGS (a system binary, not an R package) to be installed; see https://sourceforge.net/projects/mcmc-jags/.

This function previously required the caller to supply a fixedcuts argument, two cutpoints used to identify the location and scale of the group parameters. That argument has been removed: it was never a channel for incorporating genuinely known cutpoints (e.g. published cut scores on the native test-score scale), only a coordinate-system normalization required for any identified fit, and supplying cutpoints on the native test-score scale reliably sent the MCMC sampler into a degenerate region it could not escape rather than producing an error. Cutpoints are now derived internally from the pooled bin proportions, matching mle_hetop(), which was fixed the same way for the same reason.

Usage

fh_hetop(ngk, p, m, gridL, gridU, Xm=NULL, Xs=NULL,
seed=12345, modelfileonly = FALSE, modloc=NULL, ...)

Arguments

Details

The function is basically a wrapper for R2jags::jags, building model code depending on the specification of the Efron priors and any covariates for the group means and group standard deviations. Details on the FH-HETOP model are provided by Lockwood, Castellano and Shear (2018).

Covariates to predict the group means and group log standard deviations are optional. However, Xm and Xs must both be either NULL, or specified; the current version of this function cannot use covariates to predict one set of parameters but not use any covariates to predict the other set. While covariates in general must be present or absent simultaneously for the two sets of parameters, it is not necessary that the same covariates be used to predict the two sets of parameters. All covariates must be centered so that they sum to zero across groups.

The location and scale of the group means are identified for the purpose of conducting the estimation by fixing two of the cutpoints. This function derives the two fixed cutpoints internally from the pooled bin proportions via qnorm(cumsum(colSums(ngk)/sum(ngk))[1:2]), the same expression mle_hetop() uses, which places them on the same standardized scale as the model's internal latent parameters. This is a normalization needed for any identified fit, not a way to supply cutpoints that are known in advance on the test-score scale; see mle_hetop()'s documentation for further discussion.

Value

A object of class rjags, with additional information specific to the FH-HETOP model. The additional information is stored as a list called fh_hetop_extras with the following components:

Deprecated

fh_hetop() is deprecated in favor of fast_hetop(), which runs much faster, requires no external dependencies (no JAGS installation), and on real data has been found to be at least as accurate. fh_hetop() is retained in the package only for comparison purposes and receives no further development.

Author(s)

J.R. Lockwood jrlockwood@ets.org (original implementation); Paul T. von Hippel ph3828@eid.utexas.edu (removal of the fixedcuts argument in favor of internal derivation, matching mle_hetop()).

References

Efron B. (2016). “Empirical Bayes deconvolution estimates,” Biometrika 103(1):1–20.

Lockwood J.R., Castellano K.E. and 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.

See Also

R2jags::jags

Examples

## Not run: 
## fh_hetop() requires JAGS, an external system binary; see
## https://sourceforge.net/projects/mcmc-jags/.  The example below
## is wrapped in \dontrun{} so that it is not executed by R CMD
## check, but should run interactively once JAGS is installed.

set.seed(1001)

## define mean-centered covariates
G  <- 12
z1 <- sample(c(0,1), size=G, replace=TRUE)
z2 <- 0.5*z1 + rnorm(G)
Z  <- cbind(z1 - mean(z1), z2 = z2 - mean(z2))

## define true parameters dependent on covariates
beta_m    <- c(0.3,  0.8)
beta_s    <- c(0.1, -0.1)
mug       <- Z[,1]*beta_m[1] + Z[,2]*beta_m[2] + rnorm(G, sd=0.3)
sigmag    <- exp(0.3 + Z[,1]*beta_s[1] + Z[,2]*beta_s[2] + 0.2*rt(G, df=7))
cutpoints <- c(-1.0, 0.0, 1.2)

## generate data
ng   <- rep(200,G)
ngk  <- gendata_hetop(G, K = 4, ng, mug, sigmag, cutpoints)
print(ngk)

## fit FH-HETOP model including covariates
## NOTE: using an extremely small number of iterations for testing,
##       so that convergence is not expected
m <- fh_hetop(ngk, p = c(10,10),
              m = c(100, 100), gridL = c(-5.0, log(0.10)),
              gridU = c(5.0, log(5.0)), Xm = Z, Xs = Z,
              n.iter = 100, n.burnin = 50)

print(m)
print(names(m$fh_hetop_extras))

s <- m$BUGSoutput$summary
print(data.frame(truth = c(beta_m, beta_s), s[grep("beta", rownames(s)),]))

print(cor(mug,    s[grep("mu",    rownames(s)),"mean"]))
print(cor(sigmag, s[grep("sigma", rownames(s)),"mean"]))

## manual calculation of WAIC (see help file for waic_hetop)
tmp <- waic_hetop(ngk, m$BUGSoutput$sims.matrix)
identical(tmp, m$fh_hetop_extras$waicinfo)

## End(Not run)

Generate count data from Heteroskedastic Ordered Probit (HETOP) Model

Description

Generates count data for G groups and K ordinal categories under a heteroskedastic ordered probit model, given the total number of units in each group and parameters determining the category probabilities for each group.

Usage

gendata_hetop(G, K, ng, mug, sigmag, cutpoints)

Arguments

Details

For each group g, the function generates ng IID normal random variables with mean mug[g] and standard deviation sigmag[g], and then assigns each to one of K ordered groups, depending on cutpoints. The resulting data for a group is a table of category counts summing to ng[g].

Value

A G x K matrix where column k of row g provides the number of simulated units from group g falling into category k.

Author(s)

J.R. Lockwood jrlockwood@ets.org

References

Reardon S., Shear B.R., Castellano K.E. and 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.

Lockwood J.R., Castellano K.E. and 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.

Examples

set.seed(1001)

## define true parameters
G         <- 10
mug       <- seq(from= -2.0, to= 2.0, length=G)
sigmag    <- seq(from=  2.0, to= 0.8, length=G)
cutpoints <- c(-1.0, 0.0, 0.8)

## generate data with large counts
ng   <- rep(100000,G)
ngk  <- gendata_hetop(G, K = 4, ng, mug, sigmag, cutpoints)
print(ngk)

## compare theoretical and empirical cell probabilities
phat  <- ngk / ng
ptrue <- t(sapply(1:G, function(g){
    tmp <- c(pnorm(cutpoints, mug[g], sigmag[g]), 1)
    c(tmp[1], diff(tmp))
}))
print(max(abs(phat - ptrue)))

Maximum Likelihood Estimation of Heteroskedastic Ordered Probit (HETOP) Model

Description

Computes MLEs of G group means and standard deviations using count data from K ordinal categories under a heteroskedastic ordered probit model. Estimation is conducted conditional on two fixed cutpoints, and additional constraints on group parameters are imposed if needed to achieve identification in the presence of sparse counts.

This implementation is forked from the HETOP package by J. R. Lockwood (CRAN, last released 2019). We have modified the original code in two ways: (1) the inner cell-probability loop is vectorized, which substantially speeds up the runtime per likelihood evaluation; and (2) the user-facing arguments fixedcuts and svals have been removed, because some users found them confusing and supplying incompatible values caused silent optimization failures. Cutpoints and starting values are now derived internally from the data.

Usage

mle_hetop(ngk, iterlim = 1500, ...)

Arguments

Details

This function requires K >= 3. If ngk has all nonzero counts, all model parameters are identified. Alternatively, arbitrary identification rules are required to ensure the existence of the MLE when there are one or more groups with nonzero counts in fewer than three categories. This function adopts the following rules. For any group with nonzero counts in fewer than three categories, the log of the group standard deviation is constrained to equal the mean of the log standard deviations for the remaining groups. Further constraints are imposed to handle groups for which all data fall into either the lowest or highest category. Let S be the set of groups for which it is not the case that all data fall into an extreme category. Then for any group with all data in the lowest category, the mean for that group is constrained to be the minimum of the group means over S. Similarly, for any group with all data in the highest category, the mean for that grou is constrained to be the maximum of the group means over S.

The location and scale of the group means are identified for the purpose of conducting the estimation by fixing two of the cutpoints. This function derives the two fixed cutpoints internally from the pooled bin proportions via qnorm(cumsum(colSums(ngk)/sum(ngk))[1:2]), which places them on the same standardized scale as the internal starting values for the group means and log standard deviations. However in practice it may be desirable to express the group means and standard deviations on a scale that is more easily interpreted; see Reardon et al. (2017) for details. This function reports estimates on four different scales: (1) the original estimation scale with two fixed cutpoints; (2) a scale defined by forcing the group means and log group standard deviations each to have weighted mean of zero, where weights are proportional to the total count for each group; (3) a scale where the population mean of the latent variable is zero and the population standard deviation is one; and (4) a scale similar to (3) but where a bias correction is applied. See Reardon et al. (2017) for details on this bias correction.

The function also returns an estimated intracluster correlation (ICC) of the latent variable, defined as the ratio of the between-group variance of the latent variable to its marginal variance. Scales (1)-(3) above lead to the same estimated ICC; scale (4) uses a bias-corrected estimate of the ICC which will not in general equal the estimate from scales (1)-(3).

Value

A list with the following components:

Deprecated

mle_hetop() is deprecated in favor of fast_hetop(), which runs much faster, produces an estimate for every identified group, and on real data has been found to be at least as accurate. mle_hetop() is retained in the package only for comparison purposes and receives no further development.

Author(s)

J. R. Lockwood (original implementation); David J. Hunter and Paul T. von Hippel ph3828@eid.utexas.edu (vectorization and API simplification).

References

Reardon S., Shear B.R., Castellano K.E. and 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.

Lockwood J.R., Castellano K.E. and 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.

Examples

set.seed(1001)

## define true parameters
G         <- 10
mug       <- seq(from= -2.0, to= 2.0, length=G)
sigmag    <- seq(from=  2.0, to= 0.8, length=G)
cutpoints <- c(-1.0, 0.0, 0.8)

## generate data with large counts
ng   <- rep(100000,G)
ngk  <- gendata_hetop(G, K = 4, ng, mug, sigmag, cutpoints)
print(ngk)

## compute MLE and check parameter recovery (cutpoints derived from data):
m    <- mle_hetop(ngk)
print(cbind(true = mug,       est = m$est_fc$mug))
print(cbind(true = sigmag,    est = m$est_fc$sigmag))
print(cbind(true = cutpoints, est = m$est_fc$cutpoints))

## estimates on other scales:
p    <- ng/sum(ng)
print(sum(p * m$est_zero$mug))
print(sum(p * log(m$est_zero$sigmag)))

print(sum(p * m$est_star$mug))
print(sum(p * (m$est_star$mug^2 + m$est_star$sigmag^2)))

## dealing with sparse counts
ngk_sparse <- matrix(rpois(G*4, lambda=5), ncol=4)
ngk_sparse[1,] <- c(5,8,0,0)
ngk_sparse[2,] <- c(0,10,10,0)
ngk_sparse[3,] <- c(12,0,0,0)
ngk_sparse[4,] <- c(0,0,0,10)
print(ngk_sparse)

m    <- mle_hetop(ngk_sparse)
print(m$pstatus)
print(unique(m$est_fc$sigmag[1:4]))
print(exp(mean(log(m$est_fc$sigmag[5:10]))))
print(m$est_fc$mug[3])
print(min(m$est_fc$mug[-3]))
print(m$est_fc$mug[4])
print(max(m$est_fc$mug[-4]))

Shen and Louis (1998) Triple Goal Estimators

Description

triple_goal() implements the “Triple Goal” estimates of Shen and Louis (1998) for a vector of parameters given a sample from the posterior distribution of those parameters. Also computes “constrained Bayes” estimators of Ghosh (1992).

Usage

triple_goal(s, stop.if.ties = FALSE, quantile.type = 7)

Arguments

Details

In typical applications, the matrix s will be a sample of size n from the joint posterior distribution of a vector of K group-specific parameters. Both the triple goal and constrained Bayes estimators are designed to mitigate problems arising from underdispersion of posterior means; see references.

Value

A dataframe with K rows with fields:

Author(s)

J.R. Lockwood jrlockwood@ets.org

References

Shen W. and Louis T.A. (1998). “Triple-goal estimates in two-stage hierarchical models,” Journal of the Royal Statistical Society, Series B 60(2):455-471.

Ghosh M. (1992). “Constrained Bayes estimation with applications,” Journal of the American Statistical Association 87(418):533-540.

Examples

set.seed(1001)
.K <- 50
.nsamp <- 500
.theta_true <- rnorm(.K)
.s <- matrix(.theta_true, ncol=.K, nrow=.nsamp, byrow=TRUE) +
      matrix(rnorm(.K*.nsamp, sd=0.4), ncol=.K, nrow=.nsamp)
.e <- triple_goal(.s)
str(.e)
head(.e)

Texas STAAR Grade 6 Mathematics, 2017-18: District-Level Bin Counts

Description

District-level counts of students in each of four proficiency categories on the Texas State of Texas Assessments of Academic Readiness (STAAR) Grade 6 mathematics test, 2017-18 administration. For each district the dataset also reports the average score across all tested students, which can be used as ground truth for evaluating estimators that recover district means from binned counts.

Usage

data(tx_g6_math_2018)

Format

A data frame with 1151 rows and 8 columns:

district_id

Sequential integer identifier (1 to 1151).

district_name

District name (character).

n_tested

Total students tested in the district.

unsatisfactory

Students scoring below 1536 (proficiency category "Did Not Meet Grade Level").

approaches

Students scoring in [1536, 1653) ("Approaches Grade Level").

meets

Students scoring in [1653, 1772) ("Meets Grade Level").

masters

Students scoring >= 1772 ("Masters Grade Level").

reported_mean

District average score, computed by the Texas Education Agency from individual student scores.

Details

The three published cut scores defining the bin boundaries are 1536, 1653, and 1772. The administrative floor of the STAAR scale is 1062 and the ceiling is 2143. Of the 1151 districts, 1014 have nonzero counts in all four bins, 120 have nonzero counts in three bins, and 17 have nonzero counts in two bins.

Source

Texas Education Agency, Academic Performance Reports (TAPR), 2017-18. Compiled by D.\ J.\ Hunter and P.\ T.\ von Hippel.

Examples

data(tx_g6_math_2018)
str(tx_g6_math_2018)

## Recover district means using fast_hetop with known cutpoints.
ngk <- with(tx_g6_math_2018,
            cbind(unsatisfactory, approaches, meets, masters))
## These are administrative counts covering every tested student, so
## each district's students are its whole population: scope = "population".
fit <- fast_hetop(ngk, cutpoints_known = TRUE,
                  cutpoints = c(1536, 1653, 1772),
                  scope = "population")

## Correlation with reported truth on the test-score scale.
## (The 17 districts with fewer than three populated bins cannot
## identify both a mean and an SD; by default fast_hetop() estimates
## their means using a borrowed SD.  Pass
## estimate_unidentified_districts = FALSE to get NA instead, in which
## case use = "complete.obs" is needed here.)
cor(fit$est_raw$mean, tx_g6_math_2018$reported_mean)

WAIC for FH-HETOP model

Description

Computes the Watanabe-Akaike information criterion (WAIC) for the FH-HETOP model using the data and posterior samples of the group means, group standard deviations and cutpoints.

Usage

waic_hetop(ngk, samps)

Arguments

Details

Although this function can be called directly by the user, it is primarily intended to be used to compute WAIC as part of the function fh_hetop(). Details on the WAIC calculation are provided by Vehtari and Gelman (2017).

Value

A list with the following components:

Author(s)

J.R. Lockwood jrlockwood@ets.org

References

Lockwood J.R., Castellano K.E. and 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.

Vehtari A., Gelman A. and Gabry J. (2017). “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC,” Statistics and Computing. 27(5):1413–1432.

Examples

if (requireNamespace("R2jags", quietly = TRUE)) {
  set.seed(42)
  G <- 10
  ngk <- gendata_hetop(G = G, K = 4, ng = rep(50, G),
                       mug = rnorm(G), sigmag = exp(rnorm(G, 0, 0.2)),
                       cutpoints = c(-1, 0, 1))
  m <- fh_hetop(ngk,
                p = c(10, 10), m = c(100, 100),
                gridL = c(-5, log(0.10)), gridU = c(5, log(5.0)),
                n.iter = 200, n.burnin = 100, seed = 1)
  waic <- waic_hetop(ngk, m$BUGSoutput$sims.matrix)
  print(waic)
}