| Title: | Ranking-Optimized Population Posterior Expected Percentiles |
|---|---|
| Description: | Implements the empirical Bayes ranking method described in Henderson and Hartman (2025) <doi:10.48550/arXiv.2511.16530>. The data structure of interest is a collection of cluster-specific effect size estimates, associated standard errors for these estimates, and cluster-level covariates. Estimates of the rankings of the cluster-specific effect sizes after adjusting for cluster-level covariates are generated. |
| Authors: | Nicholas Henderson [cre, aut], Nicholas Hartman [aut] |
| Maintainer: | Nicholas Henderson <[email protected]> |
| License: | GPL-2 |
| Version: | 0.2 |
| Built: | 2026-07-22 09:07:32 UTC |
| Source: | https://github.com/nchenderson/ropper |
Computes estimates of cluster-level posterior expected means given an estimate of the regression coefficients and an estimate of the random effects variance. For a fixed value of the random effects variance, the collection of posterior expected mean estimates is commonly referred to as the "best linear unbiased predictor", or BLUP.
BLUP(beta.hat, Y, X, ses, tau.sq)BLUP(beta.hat, Y, X, ses, tau.sq)
beta.hat |
The 1 x p vector of regression coefficient estimates. |
Y |
The K x 1 vector of cluster-specific effect estimates. |
X |
The K x p design matrix containing cluster-specific covariates. |
ses |
The K x 1 vector sampling standard errors squared. |
tau.sq |
Estimate of the random effects variance. This must be a scalar greater than 0. |
The BLUP for cluster k is defined as
, where
is the vector of regression coefficient estimates and
is defined as
A K x 1 vector of posterior mean estimates.
Nicholas C. Henderson and Nicholas Hartman
See also PEPP.
####################################################### ### Example of using PEPP with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Estimate regression coefficients by maximum likelihood ## estimation assuming that tau.sq = 0.5 beta.mle <- solve(crossprod(XX, (stderrs_sq + 0.5)*XX), crossprod(XX, (stderrs_sq + 0.5)*yy)) ## Extract the posterior means associated with the maximum ## likelihood estimate: blup.mle <- BLUP(beta.mle, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Estimate regression coefficients using ROPPER assuming ## that tau.sq = 0.5 rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Extract posterior means associated with the ROPPER estimate blup.rop <- BLUP(rfit$coefficients, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) plot(blup.mle, blup.rop, xlab="MLE Posterior Means", ylab="ROPPER Posterior Means")####################################################### ### Example of using PEPP with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Estimate regression coefficients by maximum likelihood ## estimation assuming that tau.sq = 0.5 beta.mle <- solve(crossprod(XX, (stderrs_sq + 0.5)*XX), crossprod(XX, (stderrs_sq + 0.5)*yy)) ## Extract the posterior means associated with the maximum ## likelihood estimate: blup.mle <- BLUP(beta.mle, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Estimate regression coefficients using ROPPER assuming ## that tau.sq = 0.5 rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Extract posterior means associated with the ROPPER estimate blup.rop <- BLUP(rfit$coefficients, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) plot(blup.mle, blup.rop, xlab="MLE Posterior Means", ylab="ROPPER Posterior Means")
Computes an estimate of the random effects variance using a nonparametric, k nearest neighbors approach. In this approach, the standard errors for the cluster-level effect estimates are assumed to be fixed.
KnnTausqEst(Y, X, ses, trim.quant = 0, k = 1)KnnTausqEst(Y, X, ses, trim.quant = 0, k = 1)
Y |
The K x 1 vector of cluster-specific effect estimates. |
X |
The K x p design matrix containing cluster-specific covariates. |
ses |
The K x 1 vector sampling standard errors squared. |
trim.quant |
The proportion of extreme estimates to exclude when estimating the random effects
variance. Cluster-specific estimates lower than the |
k |
The number of neighbors to include in the k nearest neighbors procedure. This should be an integer both greater than or equal to 1 and less than or equal to K. |
Calculation of the nearest neighbor estimator of the random effects variance involves a random
split of the data into a training set and a test set of approximately equal sizes.
The one-nearest neighbor estimator (i.e., when k=1) of is
then calculated as
, where is the number of observations in the test
set , and is the effect estimate associated
with the nearest neighbor
of in the training set.
The estimate of the random effects variance. A vector of length 1.
Nicholas C. Henderson and Nicholas Hartman
Devroye, L., Györfi, L., Lugosi, G., & Walk, H. (2018). A nearest neighbor estimate of the residual variance. Electronic Journal of Statistics, 12(1), 1752–1778. doi:10.1214/18-EJS1438
See also REMLTausqEst, ropper
########################################################## ### Example of kNN estimation with simulated data ##### ########################################################## n <- 500 ## 500 units ## True tau.sq is 0.5 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Run Knn estimation with a 1 nearest neighbor approach tausq_hat_knn <- KnnTausqEst(Y=yy, X=XX, ses=stderrs_sq, k=1) tausq_hat_knn########################################################## ### Example of kNN estimation with simulated data ##### ########################################################## n <- 500 ## 500 units ## True tau.sq is 0.5 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Run Knn estimation with a 1 nearest neighbor approach tausq_hat_knn <- KnnTausqEst(Y=yy, X=XX, ses=stderrs_sq, k=1) tausq_hat_knn
Data on the number of team wins during the 2025 major league baseball (MLB) season. The data includes the following additional team characteristics: payroll, strength of schedule, and ballpark difficulty.
data("mlb2025")data("mlb2025")
A data frame with 29 observations on the following 8 variables.
Teama character vector containing team names.
Winsa numeric vector giving the number of total wins over the 2025 MLB season.
Lossesa numeric vector giving the number of total losses over the 2025 MLB season.
WinPropa numeric vector with the proportion of wins for each team.
StdErrSqa numeric vector with the standard errors squared associated with the win proportions.
Payrolla numeric vector containing the projected payroll (in millions) at the start of the 2025 baseball season.
SOSa numeric vector measuring the 2025 strength of schedule.
ParkFactorsa numeric vector measuring the overall difficulty of each team's home ballpark. Park difficulty is averaged over the 2024-2026 seasons.
All MLB teams except the Athletics were included in this data.
Total wins, total losses, and 2025 strength of schedule were obtained from https://www.baseball-reference.com/. The park factor index was retrieved from https://baseballsavant.mlb.com/. Payroll projections for 2025 were obtained from https://blogs.fangraphs.com/the-2025-payrolls-and-beyond/.
data(mlb2025) str(mlb2025) ## Fit model with Payroll, SOS, ParkFactors as covariates X1 <- model.matrix(~ Payroll + SOS + ParkFactors, data=mlb2025) rank_mod <- ropper(Y = mlb2025$WinProp, X=X1, ses=mlb2025$StdErrSq) ## Look at estimated regression coefficients rank_mod$coef ## Look at ROPPER percentiles estimates ## Rank teams according to best percentiles mlb_estperc <- data.frame(Team=mlb2025$Team, percentile=round(rank_mod$pepp, 4)) mlb_estperc <- mlb_estperc[order(-mlb_estperc$percentile),] head(mlb_estperc)data(mlb2025) str(mlb2025) ## Fit model with Payroll, SOS, ParkFactors as covariates X1 <- model.matrix(~ Payroll + SOS + ParkFactors, data=mlb2025) rank_mod <- ropper(Y = mlb2025$WinProp, X=X1, ses=mlb2025$StdErrSq) ## Look at estimated regression coefficients rank_mod$coef ## Look at ROPPER percentiles estimates ## Rank teams according to best percentiles mlb_estperc <- data.frame(Team=mlb2025$Team, percentile=round(rank_mod$pepp, 4)) mlb_estperc <- mlb_estperc[order(-mlb_estperc$percentile),] head(mlb_estperc)
Computes estimates of cluster-level posterior expected population percentiles given an estimate of the regression coefficients and an estimate of the random effects variance.
PEPP(beta.hat, Y, X, ses, tau.sq)PEPP(beta.hat, Y, X, ses, tau.sq)
beta.hat |
The 1 x p vector of regression coefficient estimates. |
Y |
The K x 1 vector of cluster-specific effect estimates. |
X |
The K x p design matrix containing cluster-specific covariates. |
ses |
The K x 1 vector sampling standard errors squared. |
tau.sq |
Estimate of the random effects variance. This must be a scalar greater than 0. |
The posterior expected population percentile for cluster k is defined as
, where
is the cumulative distribution of a standard normal distribution,
is the vector of regression coefficient estimates, and
is defined as .
The parameter is the "population percentile" which is defined
as , where is the random effect
for the cluster-k effect estimate .
A K x 1 vector of posterior expected population percentile estimates.
Nicholas C. Henderson and Nicholas Hartman
Henderson, N. C., & Hartman, N. (2025). Targeted Parameter Estimation for Robust Empirical Bayes Ranking. arXiv preprint. doi:10.48550/arXiv.2511.16530
####################################################### ### Example of using PEPP with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Estimate regression coefficients by maximum likelihood ## estimation assuming that tau.sq = 0.5 beta.mle <- solve(crossprod(XX, (stderrs_sq + 0.5)*XX), crossprod(XX, (stderrs_sq + 0.5)*yy)) ## Extract the PEPP associated with the maximum ## likelihood estimate: pepp.mle <- PEPP(beta.mle, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Get PEPP using ROPPER rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq) pepp_rop1 <- rfit$pepp ## This should match using PPEP with ROPPER estimates of ## regression coefficients and tau.sq pepp_rop2 <- PEPP(rfit$coefficients, Y=yy, X=XX, ses=stderrs_sq, tau.sq=rfit$tau.sq.hat) all.equal(pepp_rop1, pepp_rop2)####################################################### ### Example of using PEPP with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes ## Estimate regression coefficients by maximum likelihood ## estimation assuming that tau.sq = 0.5 beta.mle <- solve(crossprod(XX, (stderrs_sq + 0.5)*XX), crossprod(XX, (stderrs_sq + 0.5)*yy)) ## Extract the PEPP associated with the maximum ## likelihood estimate: pepp.mle <- PEPP(beta.mle, Y=yy, X=XX, ses=stderrs_sq, tau.sq=0.5) ## Get PEPP using ROPPER rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq) pepp_rop1 <- rfit$pepp ## This should match using PPEP with ROPPER estimates of ## regression coefficients and tau.sq pepp_rop2 <- PEPP(rfit$coefficients, Y=yy, X=XX, ses=stderrs_sq, tau.sq=rfit$tau.sq.hat) all.equal(pepp_rop1, pepp_rop2)
Computes an estimate of the random effects variance by maximizing the restricted maximum likelihood objective function. The standard errors for the cluster-level effect estimates are assumed to be fixed.
REMLTausqEst(Y, X, ses, trim.quant = 0)REMLTausqEst(Y, X, ses, trim.quant = 0)
Y |
The K x 1 vector of cluster-specific effect estimates. |
X |
The K x p design matrix containing cluster-specific covariates. |
ses |
The K x 1 vector sampling standard errors squared. |
trim.quant |
The proportion of extreme estimates to exclude when estimating the random effects
variance. Cluster-specific estimates lower than the |
Estimate of the random effects variance is found by maximizing
the restricted maximum likelihood objective function
with respect to . Here, is the K x 1 vector of
cluster-level effect estimates, is the K x K diagonal matrix
with diagonal elements , and is
the K x K matrix defined as
The estimate of the random effects variance. A vector of length 1.
Nicholas C. Henderson and Nicholas Hartman
Jiang, J., & Nguyen, T. (2021). Linear and Generalized Linear Mixed Models and Their Applications (2nd ed.). Springer. doi:10.1007/978-1-0716-1282-8
See also ropper,
########################################################## ### Example of REML estimation with simulated data ##### ########################################################## n <- 500 ## 500 units ## True tau.sq is 0.5 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) mu.vec <- xx + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes tausq_hat_reml <- REMLTausqEst(Y=yy, X=XX, ses=stderrs_sq) tausq_hat_reml########################################################## ### Example of REML estimation with simulated data ##### ########################################################## n <- 500 ## 500 units ## True tau.sq is 0.5 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) mu.vec <- xx + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes tausq_hat_reml <- REMLTausqEst(Y=yy, X=XX, ses=stderrs_sq) tausq_hat_reml
Computes ranking of cluster-specific estimates after performing a covariate adjustment.
ropper(Y, X, ses, tau.sq = c("REML", "kNN"), H=1, opt.method=c("MM", "optim"), trim.quant=0, control = list())ropper(Y, X, ses, tau.sq = c("REML", "kNN"), H=1, opt.method=c("MM", "optim"), trim.quant=0, control = list())
Y |
The K x 1 vector of cluster-specific effect estimates. |
X |
The K x p design matrix containing cluster-specific covariates. |
ses |
The K x 1 vector sampling standard errors squared. |
tau.sq |
Method for estimating the random effects variance. This can be either the
restricted maximum likelihood estimation (REML) |
H |
The degree of the approximation used for the unbiased risk estimator. This must be an integer greater than or equal to 1. |
opt.method |
The optimization method used for computing the regression coefficient estimates.
This can be either an MM algorithm or the use of the |
trim.quant |
The proportion of extreme estimates to exclude when estimating the random effects
variance. Cluster-specific estimates lower than the |
control |
A list of control parameters specifying any changes to default values of the optimization procedures. Full names of control list elements must be specified, otherwise, user-specifications are ignored. See *Details*. |
Default values of control are:
maxiter=500,
tol=1e-06
maxiterAn integer denoting the maximum limit on the number
of MM iterations or optim iterations. Default value is 500.
tolA small, positive scalar that determines when iterations
should be terminated. Iterations are terminated when
.
Default is 1.e-06.
A list with the following components:
coefficients |
The length-p vector of ranking-targeted regression coefficient estimates |
pepp |
The length-K vector of estimated posterior expected population percentiles |
objfn.track |
A vector containing the value of the objective function at each iteration |
tau.sq.hat |
The estimate of the random effects variance |
Nicholas C. Henderson and Nicholas Hartman
Henderson, N. C., & Hartman, N. (2025). Targeted Parameter Estimation for Robust Empirical Bayes Ranking. arXiv preprint. doi:10.48550/arXiv.2511.16530.
####################################################### ### Example of using ropper with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) set.seed(2507) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq) ## Regression coefficient estimates rfit$coefficients ## ROPPER percentile estimates for first 5 units rfit$pepp[1:5] ## Plot true random effect rankings vs. ROPPer rankings: true_ranks <- rank(theta.vec) ropper_ranks <- rank(rfit$pepp) plot(true_ranks, ropper_ranks, las=1, xlab="True RE Ranking", ylab="ROPPER Ranking") abline(0, 1) ###### ## Example of using ropper with optim for optimization ## instead of MM ###### rfit.optim <- ropper(Y=yy, X=XX, ses=stderrs_sq, opt.method="optim") ## Compare coefficients obtained from the two optimization methods round(cbind(rfit$coefficients, rfit.optim$coefficients), 4) ## Compare first 5 percentiles obtained from the two optimization methods round(cbind(rfit$pepp[1:5], rfit.optim$pepp[1:5]), 4)####################################################### ### Example of using ropper with simulated data ##### ####################################################### n <- 50 tau.sq <- 0.5 cluster.sizes <- 1 + rpois(n, lambda=5) set.seed(2507) xx <- rnorm(n) theta.vec <- rnorm(n, sd=sqrt(tau.sq)) ## True fixed-effect function is pnorm(x) mu.vec <- pnorm(xx) + theta.vec yy <- mu.vec + rnorm(n)/sqrt(cluster.sizes) ## Fit model where it is assumed that fixed-effect ## function is beta0 + beta1*x XX <- cbind(rep(1, n), xx) ## define standard errors squared stderrs_sq <- 1/cluster.sizes rfit <- ropper(Y=yy, X=XX, ses=stderrs_sq) ## Regression coefficient estimates rfit$coefficients ## ROPPER percentile estimates for first 5 units rfit$pepp[1:5] ## Plot true random effect rankings vs. ROPPer rankings: true_ranks <- rank(theta.vec) ropper_ranks <- rank(rfit$pepp) plot(true_ranks, ropper_ranks, las=1, xlab="True RE Ranking", ylab="ROPPER Ranking") abline(0, 1) ###### ## Example of using ropper with optim for optimization ## instead of MM ###### rfit.optim <- ropper(Y=yy, X=XX, ses=stderrs_sq, opt.method="optim") ## Compare coefficients obtained from the two optimization methods round(cbind(rfit$coefficients, rfit.optim$coefficients), 4) ## Compare first 5 percentiles obtained from the two optimization methods round(cbind(rfit$pepp[1:5], rfit.optim$pepp[1:5]), 4)