R/BTdecayLassoC.R

In BTdecayLasso: Bradley-Terry Model with Exponential Time Decayed Log-Likelihood and Adaptive Lasso

#' Bradley-Terry Model with Exponential Decayed weighted likelihood and weighted Lasso with AIC or BIC criteria
#' 
#' 
#' Model selection via AIC or BIC criteria. For Lasso estimators, the degree of freedom is the number of distinct groups of estimated abilities.
#' 
#' @param dataframe Generated using \code{\link{BTdataframe}} given raw data.
#' @param ability A column vector of teams ability, the last row is the home parameter.
#' The row number is consistent with the team's index shown in dataframe. It can be generated using \code{\link{BTdataframe}} given raw data.
#' @param weight Weight for Lasso penalty on different abilities
#' @param criteria "AIC" or "BIC"
#' @param type "HYBRID" or "LASSO"
#' @param model An Lasso path object with class wlasso or swlasso. If NULL, the whole lasso path will be run.
#' @param decay.rate The exponential decay rate. Usually ranging from (0, 0.01), A larger decay rate weights more
#' importance to most recent matches and the estimated parameters reflect more on recent behaviour.
#' @param fixed A teams index whose ability will be fixed as 0. The worstTeam's index
#' can be generated using \code{\link{BTdataframe}} given raw data.
#' @param thersh Threshold for convergence
#' @param max Maximum weight for w_{ij} (weight used for Adaptive Lasso)
#' @param iter Number of iterations used in L-BFGS-B algorithm.
#' @details
#' This function is usually run after the run of whole Lasso path. "model" parameter is obtained by whole
#' Lasso pass's run using \code{\link{BTdecayLasso}}. If no model is provided, this function will run Lasso path first (time-consuming).
#' 
#' Users can select the information score added to HYBRID Lasso's likelihood or original Lasso's likelihood. ("HYBRID" is recommended)
#' 
#' summary() function can be applied to view the outputs.
#' @return
#' \item{Score}{Lowest AIC or BIC score}
#' \item{Optimal.degree}{The degree of freedom where lowest AIC or BIC score is achieved}
#' \item{Optimal.ability}{The ability where lowest AIC or BIC score is achieved}
#' \item{ability}{Matrix contains all abilities computed in this algorithm}
#' \item{Optimal.lambda}{The lambda where lowest score is attained}
#' \item{Optimal.penalty}{The penalty (1- s/\eqn{\max(s)}) where lowest score is attained}
#' \item{type}{Type of model selection method}
#' \item{decay.rate}{Decay rate of this model}
#' @seealso \code{\link{BTdataframe}} for dataframe initialization, 
#' \code{\link{BTdecayLasso}} for obtaining a whole Lasso path
#' @references 
#' Masarotto, G. and Varin, C.(2012) The Ranking Lasso and its Application to Sport Tournaments. 
#' *The Annals of Applied Statistics* **6** 1949--1970.
#' 
#' Zou, H. (2006) The adaptive lasso and its oracle properties. 
#' *J.Amer.Statist.Assoc* **101** 1418--1429.
#' @examples 
#' ##Initializing Dataframe
#' x <- BTdataframe(NFL2010)
#' 
#' ##The following code runs the main results
#' \donttest{
#' ##Model selection through AIC
#' z <- BTdecayLassoC(x$dataframe, x$ability, weight = NULL, fixed = x$worstTeam,
#'                    criteria = "AIC", type = "LASSO")
#' summary(z)
#' 
#' ##If the whole Lasso path is run, we use it's result for model selection (recommended)
#' ##Note that the decay.rate used in model selection should be consistent with
#' ##the one which is used in whole Lasso path's run (keep the same model)
#' y1 <- BTdecayLasso(x$dataframe, x$ability, lambda = 0.1, 
#'                    decay.rate = 0.005, fixed = x$worstTeam)
#'                    
#' z1 <- BTdecayLassoC(x$dataframe, x$ability, weight = NULL, model = z1,
#'                     decay.rate = 0.005,
#'                     fixed = x$worstTeam, criteria = "BIC", type = "HYBRID")
#' }
#' 
#' @export

BTdecayLassoC <- function(dataframe, ability, weight = NULL, criteria = "AIC", type = "HYBRID", model = NULL, decay.rate = 0, 
                          fixed = 1, thersh = 1e-5, iter = 100, max = 100) {
  
 
  if (is.null(weight)) {
    weight <- BTLasso.weight(dataframe, ability, decay.rate = decay.rate, fixed = fixed, thersh = thersh, max = max, iter = iter)
  }
  
  if (is.null(model)) {
    Lp <- BTdecayLasso(dataframe, ability, lambda = NULL, weight = weight, path = TRUE, decay.rate = decay.rate,
                       fixed = fixed, thersh = thersh, max = max, iter = iter)
  } else if (class(model) == "wlasso" || class(model) == "swlasso") {
    Lp <- model
  } else {
    stop("Please provide a model contains whole lasso path generated by BTdecayLasso")
  }
  
  n1 <- nrow(dataframe)
  if (criteria == "AIC") {
    mul <- 2
  } else if (criteria == "BIC") {
    mul <- log(n1)
  } else {
    stop("criteria should either be AIC or BIC")
  }
  y <- Lp$df.path * mul
  
  x <- 2 * Lp$HYBRID.likelihood.path + y
  ind <- which(x == min(x))
  ind <- max(ind)
  dg <- Lp$df.path[ind]

  if (type == "HYBRID") {
    
    output <- list(Score = min(x), Optimal.degree = dg, Optimal.ability = as.matrix(Lp$HYBRID.ability.path[, ind]), 
                   Optimal.lambda = Lp$Lambda.path[ind], Optimal.penalty = Lp$penalty.path[ind], type = type, decay.rate = decay.rate)
    
  } else if (type == "LASSO") {
    
    m0 <- Lp$Lambda.path[ind]
    m1 <- Lp$Lambda.path[ind + 1]
    k <- m0 - m1
    j <- 1
    
    while (k > thersh * 10) {
      k <- k * 0.5
      m1 <- m0 - k
      BT <- BTdecayLasso.step2(dataframe, ability, m1, weight, decay.rate = decay.rate, fixed = fixed, thersh = thersh, iter = iter)
      if (dg == BT$df) {
        m0 <- m1
      }
      j <- j + 1
    }
    
    if (j == 1) {
      output <- list(Score = 2 * Lp$likelihood.path[ind] + dg * mul, Optimal.degree = dg, Optimal.ability = as.matrix(Lp$ability.path[, ind]),
                     Optimal.lambda = m0, Optimal.penalty = Lp$penalty.path[ind], type = type, decay.rate = decay.rate)
    } else {
      l <- BTLikelihood(dataframe, BT$ability, decay.rate = decay.rate)
      output <- list(Score = 2 * l + dg * mul, Optimal.degree = dg, Optimal.ability = BT$ability,
                     Optimal.lambda = m1, Optimal.penalty = BT$penalty, type = type, decay.rate = decay.rate)
    }
    
  } else {
    stop("Please provide a selection type HYBRID or LASSO")
  }
  class(output) <- "BTC"
  output
}