R/mvglmmRank.R

In mvglmmRank: Multivariate Generalized Linear Mixed Models for Ranking Sports Teams

#' 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.
#'
#' @param game.data A data frame with columns \code{home}, \code{away},
#'   \code{home.response}, and \code{away.response}. The optional column
#'   \code{binary.response} should contain binary home-team outcomes; if it is
#'   omitted, \code{mvglmmRank()} creates a home-win indicator from the score
#'   columns. The optional column \code{neutral.site} should be \code{1} for
#'   neutral-site games and \code{0} otherwise; if omitted, all games are
#'   treated as non-neutral. If \code{OT.flag = TRUE}, \code{game.data} should
#'   also contain a numeric \code{OT} column with no missing values.
#' @param method Character string naming the model to fit. Choices are
#'   \code{"B"}, \code{"P0"}, \code{"P1"}, \code{"N"}, \code{"NB"},
#'   \code{"PB0"}, \code{"PB1"}, \code{"NB.mov"}, and \code{"N.mov"}.
#' @param first.order Logical. If \code{TRUE}, use only the first-order Laplace
#'   approximation. If \code{FALSE}, continue to the fully exponential Laplace
#'   approximation where supported.
#' @param home.field Logical. If \code{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 \code{FALSE}, score models use a single mean score
#'   and binary models use no fixed home-field effect. For margin-of-victory
#'   methods, \code{home.field = FALSE} omits the fixed home margin effect.
#' @param max.iter.EM Maximum number of EM iterations.
#' @param tol1 Convergence tolerance for the first-order Laplace
#'   approximation, based on the maximum relative parameter change.
#' @param tol2 Convergence tolerance for the fully exponential Laplace
#'   approximation. Not used when \code{first.order = TRUE}.
#' @param tolFE Intermediate convergence tolerance for the fully exponential
#'   approximation. Corrections to the random-effects covariance matrix begin
#'   after this tolerance is reached.
#' @param tol.n Convergence tolerance for the normal models. Convergence is
#'   declared when \eqn{(l_k-l_{k-1})/l_k < tol.n}, where \eqn{l_k} is the
#'   log-likelihood at iteration \eqn{k}.
#' @param verbose Logical. If \code{TRUE}, print iteration information.
#' @param OT.flag Logical. If \code{TRUE}, include the numeric \code{OT} column
#'   as a score-model fixed effect. Use \code{0} for games without overtime.
#' @param Hessian Logical. If \code{TRUE}, approximate the Hessian of the model
#'   parameters by central differences.
#' @param REML.N Logical. If \code{TRUE}, use REML estimation for
#'   \code{method = "N"} and \code{method = "N.mov"}.
#'
#' @details
#' The available methods are:
#' \describe{
#'   \item{\code{"B"}}{Binary/probit model for home win/loss indicators.}
#'   \item{\code{"P0"}}{Poisson score model without a game-level random
#'     effect.}
#'   \item{\code{"P1"}}{Poisson score model with a game-level random effect.}
#'   \item{\code{"N"}}{Normal score model with an unstructured within-game
#'     error covariance matrix.}
#'   \item{\code{"NB"}}{Joint normal score and binary/probit win/loss model.}
#'   \item{\code{"PB0"}}{Joint Poisson score and binary/probit win/loss model
#'     without a game-level random effect.}
#'   \item{\code{"PB1"}}{Joint Poisson score and binary/probit win/loss model
#'     with a game-level random effect.}
#'   \item{\code{"NB.mov"}}{Joint normal margin-of-victory and binary/probit
#'     win/loss model.}
#'   \item{\code{"N.mov"}}{Normal margin-of-victory model.}
#' }
#'
#' Neutral-site games are represented in \code{game.data$neutral.site}. Use
#' \code{1} for neutral-site games and \code{0} otherwise. For neutral-site
#' games, the teams may be assigned to the \code{home} and \code{away} columns
#' arbitrarily. With \code{home.field = TRUE}, score models estimate a
#' neutral-site mean score when neutral-site games are present. With
#' \code{home.field = FALSE}, the home/away and neutral-site mean structure is
#' suppressed.
#'
#' Setting \code{first.order = TRUE} yields the first-order Laplace
#' approximation. A partial fully exponential Laplace approximation can be
#' obtained by setting \code{tol1 > tol2} and \code{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.
#'
#' The \code{"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.
#'
#' @return
#' An object of class \code{"mvglmmRank"}. The object is a list whose
#' components depend on \code{method} and may include:
#' \describe{
#'   \item{\code{n.ratings.offense}, \code{n.ratings.defense}}{Normal-model
#'     offensive and defensive ratings, or \code{NULL}.}
#'   \item{\code{p.ratings.offense}, \code{p.ratings.defense}}{Poisson-model
#'     offensive and defensive ratings, or \code{NULL}.}
#'   \item{\code{b.ratings}}{Binary/probit win-propensity ratings, or
#'     \code{NULL}.}
#'   \item{\code{n.ratings.mov}}{Normal margin-of-victory ratings, or
#'     \code{NULL}.}
#'   \item{\code{n.mean}, \code{p.mean}, \code{b.mean}}{Estimated fixed-effect
#'     means or home-field effects for the fitted model components.}
#'   \item{\code{G}, \code{G.cor}}{Random-effects covariance and correlation
#'     matrices.}
#'   \item{\code{R}, \code{R.cor}}{Normal-model error covariance and
#'     correlation matrices, or \code{NULL}.}
#'   \item{\code{home.field}}{Logical indicating whether a home-field effect
#'     was modeled.}
#'   \item{\code{Hessian}}{Numerical Hessian if requested, otherwise
#'     \code{NULL}.}
#'   \item{\code{parameters}}{Vector of fitted model parameters.}
#'   \item{\code{actual}, \code{pred}, \code{sresid}}{Observed values,
#'     fitted values, and scaled residuals where available.}
#'   \item{\code{N.output}}{Additional normal-model matrices and covariance
#'     output for \code{method = "N"} and \code{method = "N.mov"}.}
#'   \item{\code{fixed.effect.model.output}}{Additional fixed-effect
#'     margin-of-victory output for \code{method = "N.mov"}.}
#'   \item{\code{method}}{The model method supplied by the user.}
#' }
#'
#' @references
#' Broatch, J.E. and Karl, A.T. (2018). Multivariate Generalized Linear Mixed
#' Models for Joint Estimation of Sporting Outcomes. \emph{Italian Journal of
#' Applied Statistics}, 30(2), 189-211. Also available from
#' \url{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. \emph{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. \emph{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. \emph{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. \emph{Journal of Quantitative Analysis in Sports}, 8(3).
#' \doi{10.1515/1559-0410.1471}.
#'
#' @seealso \code{\link{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")
#'
#' \donttest{
#' result <- mvglmmRank(nfl2012, method = "PB0", first.order = TRUE,
#'                      verbose = FALSE)
#' print(result)
#' game.pred(result, home = "Denver Broncos", away = "Green Bay Packers")
#' }
#'
#' @export
#' @keywords regression
mvglmmRank <-
function (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) 
{
    if (!inherits(class("method"), "character")) {
        cat("*Error: method must be a character string (using quotation marks)")
        flush.console()
        return(0)
    }
    if (first.order) 
        control <- list(iter.EM = max.iter.EM, tol1 = tol1, tol2 = tol2, 
            tolFE = tolFE, verbose = verbose,OT.flag=OT.flag,Hessian=Hessian)
    if (!first.order) 
        control <- list(iter.EM = max.iter.EM, tol1 = tol1, tol2 = tol2, 
            tolFE = tolFE, verbose = verbose,OT.flag=OT.flag, Hessian=Hessian)
    control.n <- list(iter.EM = max.iter.EM, tol1 = tol.n, verbose = verbose,OT.flag=OT.flag,Hessian=Hessian,REML.N=REML.N)
    Z_mat <- game.data
    Z_mat$Score.For <- Z_mat$home.response
    Z_mat$Score.Against <- Z_mat$away.response
    if(is.null(Z_mat$binary.response)){
    Z_mat$home_win <- as.numeric(Z_mat$Score.For > Z_mat$Score.Against)
    if(sum(Z_mat$Score.For == Z_mat$Score.Against)>0){
    Z_mat[Z_mat$Score.For == Z_mat$Score.Against,]$home_win<-2
    }
    }else{
    Z_mat$home_win<-Z_mat$binary.response
    }
    Z_mat$home <- as.character(Z_mat$home)
    Z_mat$away <- as.character(Z_mat$away)
    if (is.null(Z_mat$neutral.site)) 
        Z_mat$neutral.site <- 0
    if (method == "B") {
        res <- binary_cre(Z_mat = Z_mat, first.order = first.order, home.field=home.field, 
            control = control)
    }
    else if (method == "P0") {
        res <- poisson_cre(Z_mat = Z_mat, first.order = first.order,home.field=home.field, 
            control = control, game.effect = FALSE)
    }
    else if (method == "P1") {
        res <- poisson_cre(Z_mat = Z_mat, first.order = first.order,home.field=home.field, 
            control = control, game.effect = TRUE)
    }
    else if (method == "N") {
        res <- normal_cre(Z_mat = Z_mat, first.order = first.order,home.field=home.field, 
            control = control.n)
    }
    else if (method == "NB") {
        res <- NB_cre(Z_mat = Z_mat, first.order = first.order, home.field=home.field,
            control = control)
    }
    else if (method == "PB0") {
        res <- PB_cre(Z_mat = Z_mat, first.order = first.order, home.field=home.field, 
            control = control, game.effect = FALSE)
    }
    else if (method == "PB1") {
        res <- PB_cre(Z_mat = Z_mat, first.order = first.order,  home.field=home.field,
            control = control, game.effect = TRUE)
    }
    else if (method == "NB.mov") {
        res <- NB_mov(Z_mat = Z_mat, first.order = first.order, home.field=home.field,
            control = control)
    }
    else if (method == "N.mov") {
      res <- N_mov(Z_mat = Z_mat, first.order = TRUE, home.field=home.field,
                    control = control.n)
    }
    else {
        cat("Error in specification of method. This field is case sensitive.\n")
        return(0)
    }
    res<-c(res,method=method)
    class(res) <- "mvglmmRank"
    return(res)
}