| Type: | Package |
| Title: | Multivariate Generalized Linear Mixed Models for Ranking Sports Teams |
| Version: | 1.2-6 |
| Depends: | R (≥ 3.2.0), Matrix |
| Imports: | numDeriv, methods, stats, utils, MASS |
| Date: | 2026-09-18 |
| Description: | Maximum likelihood estimates are obtained via an EM algorithm with either a first-order or a fully exponential Laplace approximation as documented by Broatch and Karl (2018) <doi:10.48550/arXiv.1710.05284>, Karl, Yang, and Lohr (2014) <doi:10.1016/j.csda.2013.11.019>, and by Karl (2012) <doi:10.1515/1559-0410.1471>. Karl and Zimmerman <doi:10.1016/j.jspi.2020.06.004> use this package to illustrate how the home field effect estimator from a mixed model can be biased under nonrandom scheduling. |
| ByteCompile: | yes |
| LazyLoad: | yes |
| LazyData: | yes |
| License: | GPL-2 |
| Encoding: | UTF-8 |
| NeedsCompilation: | no |
| Author: | Andrew T. Karl |
| Maintainer: | Andrew T. Karl <akarl@asu.edu> |
| Config/roxygen2/version: | 8.1.0 |
| Suggests: | testthat (≥ 3.0.0) |
| Config/testthat/edition: | 3 |
| Packaged: | 2026-09-18 17:18:53 UTC; andre |
| Repository: | CRAN |
| Date/Publication: | 2026-09-18 17:40:02 UTC |
mvglmmRank: Multivariate generalized linear mixed models for ranking sports teams
Description
The package fits multivariate generalized linear mixed models for team scores, win/loss indicators, and margin-of-victory responses. Maximum likelihood estimates are obtained by an EM algorithm using either a first-order or fully exponential Laplace approximation.
Details
See mvglmmRank for the fitting interface and
game.pred for printed game predictions from fitted models.
Author(s)
Maintainer: Andrew T. Karl akarl@asu.edu (ORCID)
Authors:
Andrew T. Karl akarl@asu.edu (ORCID)
Jennifer Broatch
References
Broatch, J.E. and Karl, A.T. (2018). Multivariate Generalized Linear Mixed Models for Joint Estimation of Sporting Outcomes. Italian Journal of Applied Statistics, 30(2), 189-211. Also available from https://arxiv.org/abs/1710.05284.
Karl, A.T. and Zimmerman, D.L. (2021). A Diagnostic for Bias in Linear Mixed Model Estimators Induced by Dependence Between the Random Effects and the Corresponding Model Matrix. Journal of Statistical Planning and Inference, 211, 107-118. doi:10.1016/j.jspi.2020.06.004.
Karl, A.T., Yang, Y. and Lohr, S. (2014). Computation of Maximum Likelihood Estimates for Multiresponse Generalized Linear Mixed Models with Non-nested, Correlated Random Effects. Computational Statistics & Data Analysis, 73, 146-162. doi:10.1016/j.csda.2013.11.019.
Karl, A.T. (2012). The Sensitivity of College Football Rankings to Several Modeling Choices. Journal of Quantitative Analysis in Sports, 8(3). doi:10.1515/1559-0410.1471.
Internal NB_cre model fitting routine
Description
Implements the EM/Laplace calculations selected by mvglmmRank().
Usage
NB_cre(
Z_mat = Z_mat,
first.order = first.order,
home.field = home.field,
control = control
)
Arguments
Z_mat |
Prepared game data from |
first.order |
Use the first-order Laplace approximation. |
home.field |
Include home-field fixed effects. |
control |
Validated iteration, tolerance, and output controls. |
Value
A list of fitted parameters, ratings, and diagnostics.
Internal NB_mov model fitting routine
Description
Implements the EM/Laplace calculations selected by mvglmmRank().
Usage
NB_mov(
Z_mat = Z_mat,
first.order = first.order,
home.field = home.field,
control = control
)
Arguments
Z_mat |
Prepared game data from |
first.order |
Use the first-order Laplace approximation. |
home.field |
Include home-field fixed effects. |
control |
Validated iteration, tolerance, and output controls. |
Value
A list of fitted parameters, ratings, and diagnostics.
Internal normal margin-of-victory model
Description
Internal normal margin-of-victory model
Usage
N_mov(Z_mat, first.order = TRUE, home.field, control)
Arguments
Z_mat |
Validated game data with internal score and outcome columns. |
first.order |
Logical; unused for this exact Gaussian calculation. |
home.field |
Whether to include location fixed effects. |
control |
List of iteration, tolerance, overtime, Hessian and REML controls. |
Value
A fitted-model list used by mvglmmRank().
Internal PB_cre model fitting routine
Description
Implements the EM/Laplace calculations selected by mvglmmRank().
Usage
PB_cre(
Z_mat = Z_mat,
first.order = first.order,
home.field = home.field,
control = control,
game.effect = game.effect
)
Arguments
Z_mat |
Prepared game data from |
first.order |
Use the first-order Laplace approximation. |
home.field |
Include home-field fixed effects. |
control |
Validated iteration, tolerance, and output controls. |
game.effect |
Include an independent game-level random effect. |
Value
A list of fitted parameters, ratings, and diagnostics.
Internal binary_cre model fitting routine
Description
Implements the EM/Laplace calculations selected by mvglmmRank().
Usage
binary_cre(
Z_mat = Z_mat,
first.order = first.order,
home.field,
control = control
)
Arguments
Z_mat |
Prepared game data from |
first.order |
Use the first-order Laplace approximation. |
home.field |
Include home-field fixed effects. |
control |
Validated iteration, tolerance, and output controls. |
Value
A list of fitted parameters, ratings, and diagnostics.
2008 FBS College Football Regular Season Data
Description
2008 FBS College Football Regular Season Data
Usage
f2008
Format
A data frame with 772 observations on the following 9 variables.
homea factor
Game.Datea POSIXlt date variable
awaya factor
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
neutral.sitea numeric vector
partitiona numeric vector
Source
http://web1.ncaa.org/mfb/download.jsp?year=2008&div=IA
Examples
data(f2008)
str(f2008)
2009 FBS College Football Regular Season Data
Description
2009 FBS College Football Regular Season Data
Usage
f2009
Format
A data frame with 772 observations on the following 7 variables.
homea factor
Game.Datea POSIXlt date variable
awaya factor
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
neutral.sitea numeric vector
partitiona numeric vector
Source
http://web1.ncaa.org/mfb/download.jsp?year=2009&div=IA
Examples
data(f2009)
str(f2009)
2010 FBS College Football Regular Season Data
Description
2010 FBS College Football Regular Season Data
Usage
f2010
Format
A data frame with 770 observations on the following 9 variables.
homea factor
Game.Datea POSIXlt
awaya factor
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
neutral.sitea numeric vector
partitiona numeric vector
Source
http://web1.ncaa.org/mfb/download.jsp?year=2010&div=IA
Examples
data(f2010)
str(f2010)
2011 FBS College Football Regular Season Data
Description
2011 FBS College Football Regular Season Data
Usage
f2011
Format
A data frame with 781 observations on the following 9 variables.
homea factor
Game.Datea POSIXlt
awaya factor
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
neutral.sitea numeric vector
partitiona numeric vector
Source
http://web1.ncaa.org/mfb/download.jsp?year=2011&div=IA
Examples
data(f2011)
str(f2011)
2012 FBS College Football Regular Season Data
Description
2012 FBS College Football Regular Season Data
Usage
f2012
Format
A data frame with 809 observations on the following 9 variables.
homea factor
Game.Datea POSIXlt
awaya factor
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
neutral.sitea numeric vector
partitiona numeric vector
Source
http://web1.ncaa.org/mfb/download.jsp?year=2012&div=IA
Examples
data(f2012)
str(f2012)
Print predictions for a future game
Description
Uses a fitted mvglmmRank object to print predicted scores,
win probability, and/or margin of victory for a specified matchup.
Usage
game.pred(res, home, away, neutral.site = FALSE)
Arguments
res |
An object of class |
home |
Character string naming the home team. The name should match a team name in the fitted object. |
away |
Character string naming the away team. The name should match a team name in the fitted object. |
neutral.site |
Logical. If |
Details
Neutral-site predictions require the training data supplied to
mvglmmRank to contain neutral.site = 1 games. If a
fitted score model has no neutral-site mean, neutral-site score predictions
produce an error. Predictions refer to games with OT = 0; no
overtime increment is added. Unknown team names also produce an error.
Value
Prints predictions and returns NULL invisibly.
References
Broatch, J.E. and Karl, A.T. (2018). Multivariate Generalized Linear Mixed Models for Joint Estimation of Sporting Outcomes. Italian Journal of Applied Statistics, 30(2), 189-211. Also available from https://arxiv.org/abs/1710.05284.
Karl, A.T., Yang, Y. and Lohr, S. (2014). Computation of Maximum Likelihood Estimates for Multiresponse Generalized Linear Mixed Models with Non-nested, Correlated Random Effects. Computational Statistics & Data Analysis, 73, 146-162. doi:10.1016/j.csda.2013.11.019.
Karl, A.T. (2012). The Sensitivity of College Football Rankings to Several Modeling Choices. Journal of Quantitative Analysis in Sports, 8(3). doi:10.1515/1559-0410.1471.
See Also
Examples
data(nfl2012)
fit <- mvglmmRank(nfl2012, method = "PB0", first.order = TRUE,
max.iter.EM = 1, verbose = FALSE)
game.pred(fit, home = "Denver Broncos", away = "Green Bay Packers")
Exact Gaussian mixed-model computations
Description
Shared ML and REML calculations for normal scores and margins. REML uses Henderson prediction-error covariances in variance updates; ML uses conditional random-effect covariances. See Karl (2026), equations (6)–(10).
Fit multivariate generalized linear mixed models for sports rankings
Description
Fits one of several generalized linear mixed models for team scores, win/loss indicators, or margin of victory. The fitted random effects are used as team ratings.
Usage
mvglmmRank(
game.data,
method = "PB0",
first.order = FALSE,
home.field = TRUE,
max.iter.EM = 1000,
tol1 = 1e-04,
tol2 = 1e-04,
tolFE = 0,
tol.n = 1e-07,
verbose = TRUE,
OT.flag = FALSE,
Hessian = FALSE,
REML.N = TRUE
)
Arguments
game.data |
A data frame with columns |
method |
Character string naming the model to fit. Choices are
|
first.order |
Logical. If |
home.field |
Logical. If |
max.iter.EM |
Maximum number of EM iterations. |
tol1 |
Convergence tolerance for the first-order Laplace approximation, based on the maximum relative parameter change. |
tol2 |
Convergence tolerance for the fully exponential Laplace
approximation. Not used when |
tolFE |
Intermediate convergence tolerance for the fully exponential approximation. Corrections to the random-effects covariance matrix begin after this tolerance is reached. |
tol.n |
Convergence tolerance for the normal models. Convergence is
declared when the absolute log-likelihood change divided by
|
verbose |
Logical. If |
OT.flag |
Logical. If |
Hessian |
Logical. If |
REML.N |
Logical. If |
Details
The available methods are:
"B"Binary/probit model for home win/loss indicators.
"P0"Poisson score model without a game-level random effect.
"P1"Poisson score model with a game-level random effect.
"N"Normal score model with an unstructured within-game error covariance matrix.
"NB"Joint normal score and binary/probit win/loss model.
"PB0"Joint Poisson score and binary/probit win/loss model without a game-level random effect.
"PB1"Joint Poisson score and binary/probit win/loss model with a game-level random effect.
"NB.mov"Joint normal margin-of-victory and binary/probit win/loss model.
"N.mov"Normal margin-of-victory model.
Neutral-site games are represented in game.data$neutral.site. Use
1 for neutral-site games and 0 otherwise. For neutral-site
games, the teams may be assigned to the home and away columns
arbitrarily. With home.field = TRUE, score models estimate a
neutral-site mean score when neutral-site games are present. With
home.field = FALSE, the home/away and neutral-site mean structure is
suppressed.
Setting first.order = TRUE yields the first-order Laplace
approximation. A partial fully exponential Laplace approximation can be
obtained by setting tol1 > tol2 and tolFE = 0. This applies
fully exponential corrections to the vector of team ratings, but not to the
covariance matrix of this vector. Karl, Yang, and Lohr (2014) show that this
approach produces a large portion of the benefit of the fully exponential
Laplace approximation in only a fraction of the time.
A tied binary outcome contributes one home-win and one away-win
observation, retaining the package's historical treatment of ties. This
is a paired binary contribution, not a model with a separate tie probability.
With exclusively neutral-site games, methods containing a binary or margin
component require home.field = FALSE because a home effect cannot
be estimated. The normal and Poisson score-only methods can estimate a
single neutral-site score mean.
The "NB" label denotes a joint normal/binary model, not a
negative-binomial distribution. REML.N applies only to the two
Gaussian methods; the remaining methods use approximate ML.
The "PB1" method is the least scalable, as its memory and
computational requirements are at least quadratic in the number of teams
plus the number of games.
Value
An object of class "mvglmmRank". The object is a list whose
components depend on method and may include:
n.ratings.offense,n.ratings.defenseNormal-model offensive and defensive ratings, or
NULL.p.ratings.offense,p.ratings.defensePoisson-model offensive and defensive ratings, or
NULL.b.ratingsBinary/probit win-propensity ratings, or
NULL.n.ratings.movNormal margin-of-victory ratings, or
NULL.N.movretains its two-column data frame (team, rating);NB.movreturns a named numeric vector.n.mean,p.mean,b.meanEstimated fixed-effect means or home-field effects for the fitted model components.
G,G.corRandom-effects covariance and correlation matrices.
R,R.corNormal-model error covariance and correlation matrices, or
NULL.home.fieldLogical indicating whether a home-field effect was modeled.
HessianNumerical Hessian if requested, otherwise
NULL.parametersVector of fitted model parameters.
actual,pred,sresidObserved values, fitted values, and scaled residuals where available.
N.outputAdditional normal-model matrices and covariance output for
method = "N"andmethod = "N.mov".fixed.effect.model.outputAdditional fixed-effect margin-of-victory output for
method = "N.mov".logLikFinal log-likelihood: exact ML/REML for Gaussian methods, first-order Laplace for the other methods.
iterations,converged,approximationNumber of iterations, whether the requested convergence criterion was met, and the approximation actually reached. Check
convergedafter a fit with a small iteration limit.methodThe model method supplied by the user.
References
Karl, A.T. (2026). Motivating REML via Prediction-Error Covariances in EM Updates for Linear Mixed Models. https://arxiv.org/abs/2602.09247.
Broatch, J.E. and Karl, A.T. (2018). Multivariate Generalized Linear Mixed Models for Joint Estimation of Sporting Outcomes. Italian Journal of Applied Statistics, 30(2), 189-211. Also available from https://arxiv.org/abs/1710.05284.
Karl, A.T. and Zimmerman, D.L. (2021). A Diagnostic for Bias in Linear Mixed Model Estimators Induced by Dependence Between the Random Effects and the Corresponding Model Matrix. Journal of Statistical Planning and Inference, 211, 107-118. doi:10.1016/j.jspi.2020.06.004.
Karl, A.T., Yang, Y. and Lohr, S. (2013). Efficient Maximum Likelihood Estimation of Multiple Membership Linear Mixed Models, with an Application to Educational Value-Added Assessments. Computational Statistics and Data Analysis, 59, 13-27.
Karl, A.T., Yang, Y. and Lohr, S. (2014). Computation of Maximum Likelihood Estimates for Multiresponse Generalized Linear Mixed Models with Non-nested, Correlated Random Effects. Computational Statistics & Data Analysis, 73, 146-162. doi:10.1016/j.csda.2013.11.019.
Karl, A.T. (2012). The Sensitivity of College Football Rankings to Several Modeling Choices. Journal of Quantitative Analysis in Sports, 8(3). doi:10.1515/1559-0410.1471.
See Also
Examples
data(nfl2012)
fit <- mvglmmRank(nfl2012, method = "PB0", first.order = TRUE,
max.iter.EM = 1, verbose = FALSE)
game.pred(fit, home = "Denver Broncos", away = "Green Bay Packers")
result <- mvglmmRank(nfl2012, method = "PB0", first.order = TRUE,
verbose = FALSE)
print(result)
game.pred(result, home = "Denver Broncos", away = "Green Bay Packers")
2013 NBA Data
Description
2013 NBA Data
Usage
nba2013
Format
A data frame with 1229 observations on the following 11 variables.
Datea factor
awaya factor
homea factor
OTa factor
partitiona numeric vector
neutral.sitea numeric vector
ot.counta numeric vector
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
Source
http://masseyratings.com/data.php
Examples
data(nba2013)
str(nba2013)
2012 NCAA Division I Basketball Results
Description
2012 NCAA Division I Basketball Results
Usage
ncaab2012
Format
A data frame with 5253 observations on the following 10 variables.
datea factor
awaya factor
homea factor
neutral.sitea numeric vector
partitiona numeric vector
home_wina numeric vector
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
Source
http://masseyratings.com/data.php
Examples
data(ncaab2012)
str(ncaab2012)
2012 NFL Regular Season Data
Description
2012 NFL Regular Season Data
Usage
nfl2012
Format
A data frame with 256 observations on the following 9 variables.
Datea factor
awaya factor
homea factor
neutral.sitea numeric vector
home.responsea numeric vector
home.scorea numeric vector
away.responsea numeric vector
away.scorea numeric vector
partitiona numeric vector
Source
http://masseyratings.com/data.php
Examples
data(nfl2012)
str(nfl2012)
Internal normal score model
Description
Internal normal score model
Usage
normal_cre(Z_mat, first.order = TRUE, home.field, control)
Arguments
Z_mat |
Validated game data with internal score and outcome columns. |
first.order |
Logical; unused for this exact Gaussian calculation. |
home.field |
Whether to include location fixed effects. |
control |
List of iteration, tolerance, overtime, Hessian and REML controls. |
Value
A fitted-model list used by mvglmmRank().
Numerical helpers for the model fitting routines
Description
Internal utilities for stable probit derivatives, design matrices, convergence comparisons, and numerical differentiation.
Internal poisson_cre model fitting routine
Description
Implements the EM/Laplace calculations selected by mvglmmRank().
Usage
poisson_cre(
Z_mat = Z_mat,
first.order = first.order,
control = control,
game.effect = game.effect,
home.field = home.field
)
Arguments
Z_mat |
Prepared game data from |
first.order |
Use the first-order Laplace approximation. |
control |
Validated iteration, tolerance, and output controls. |
game.effect |
Include an independent game-level random effect. |
home.field |
Include home-field fixed effects. |
Value
A list of fitted parameters, ratings, and diagnostics.