Package {mvglmmRank}


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 ORCID iD [cre, aut], Jennifer Broatch [aut]
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:

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 mvglmmRank().

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 mvglmmRank().

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 mvglmmRank().

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 mvglmmRank().

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.

home

a factor

Game.Date

a POSIXlt date variable

away

a factor

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

neutral.site

a numeric vector

partition

a 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.

home

a factor

Game.Date

a POSIXlt date variable

away

a factor

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

neutral.site

a numeric vector

partition

a 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.

home

a factor

Game.Date

a POSIXlt

away

a factor

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

neutral.site

a numeric vector

partition

a 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.

home

a factor

Game.Date

a POSIXlt

away

a factor

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

neutral.site

a numeric vector

partition

a 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.

home

a factor

Game.Date

a POSIXlt

away

a factor

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

neutral.site

a numeric vector

partition

a 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 "mvglmmRank" returned by mvglmmRank.

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 TRUE and res was fit with home.field = TRUE using data that included neutral-site games, use the neutral-site score mean and omit the binary home-field effect in the printed prediction. If res was fit with home.field = FALSE, this argument has no additional effect because no home-field fixed effect was modeled.

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

mvglmmRank

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 home, away, home.response, and away.response. The optional column binary.response contains home-team outcomes coded as 0 (away win), 1 (home win), or 2 (tie). If omitted, mvglmmRank() creates a home-win indicator from the score columns. The optional column neutral.site should be 1 for neutral-site games and 0 otherwise; if omitted, all games are treated as non-neutral. Scores must be finite and, for Poisson methods, nonnegative integers. Method B needs no score columns when binary.response is supplied. If OT.flag = TRUE, game.data should also contain a numeric OT column with no missing values.

method

Character string naming the model to fit. Choices are "B", "P0", "P1", "N", "NB", "PB0", "PB1", "NB.mov", and "N.mov".

first.order

Logical. If TRUE, use only the first-order Laplace approximation. If FALSE, continue to the fully exponential Laplace approximation for non-Gaussian models. Gaussian methods use exact calculations and ignore this argument.

home.field

Logical. If TRUE, include home-field fixed effects: score models use separate home and away mean scores, plus a neutral-site mean score when neutral-site games are present, and binary models include a home-field effect. If FALSE, score models use a single mean score and binary models use no fixed home-field effect. For margin-of-victory methods, home.field = FALSE omits the fixed home margin effect.

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 first.order = TRUE.

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 \max(1, |l_k|) is less than tol.n, where l_k is the log-likelihood at iteration k.

verbose

Logical. If TRUE, print iteration information.

OT.flag

Logical. If TRUE, include the numeric OT column as a score-model fixed effect, with or without home-field effects. Use 0 for games without overtime. This option is not supported for B, N.mov, or NB.mov. The OT column must vary independently of the other fixed effects.

Hessian

Logical. If TRUE, approximate the Hessian of the model parameters by central differences of the negative log-likelihood. For non-Gaussian methods this is the first-order Laplace log-likelihood, even when fully exponential moment corrections are used for estimation. For Gaussian REML, the Hessian contains covariance parameters only; fixed-effect uncertainty is in N.output$ybetas_asycov.

REML.N

Logical. If TRUE, use REML estimation for method = "N" and method = "N.mov".

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.defense

Normal-model offensive and defensive ratings, or NULL.

p.ratings.offense, p.ratings.defense

Poisson-model offensive and defensive ratings, or NULL.

b.ratings

Binary/probit win-propensity ratings, or NULL.

n.ratings.mov

Normal margin-of-victory ratings, or NULL. N.mov retains its two-column data frame (team, rating); NB.mov returns a named numeric vector.

n.mean, p.mean, b.mean

Estimated fixed-effect means or home-field effects for the fitted model components.

G, G.cor

Random-effects covariance and correlation matrices.

R, R.cor

Normal-model error covariance and correlation matrices, or NULL.

home.field

Logical indicating whether a home-field effect was modeled.

Hessian

Numerical Hessian if requested, otherwise NULL.

parameters

Vector of fitted model parameters.

actual, pred, sresid

Observed values, fitted values, and scaled residuals where available.

N.output

Additional normal-model matrices and covariance output for method = "N" and method = "N.mov".

fixed.effect.model.output

Additional fixed-effect margin-of-victory output for method = "N.mov".

logLik

Final log-likelihood: exact ML/REML for Gaussian methods, first-order Laplace for the other methods.

iterations, converged, approximation

Number of iterations, whether the requested convergence criterion was met, and the approximation actually reached. Check converged after a fit with a small iteration limit.

method

The 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

game.pred

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.

Date

a factor

away

a factor

home

a factor

OT

a factor

partition

a numeric vector

neutral.site

a numeric vector

ot.count

a numeric vector

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a 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.

date

a factor

away

a factor

home

a factor

neutral.site

a numeric vector

partition

a numeric vector

home_win

a numeric vector

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a 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.

Date

a factor

away

a factor

home

a factor

neutral.site

a numeric vector

home.response

a numeric vector

home.score

a numeric vector

away.response

a numeric vector

away.score

a numeric vector

partition

a 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 mvglmmRank().

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.