elomov: Margin of Victory Elo Rating Systems

#' Multiplicative MOV Elo Ratings
#'
#' This function calculates MOV Elo ratings using a multiplicative model that combines information for both the win result and MOV result.
#'
#' @param winners. Character vector or formula specifying the winners of each result
#' @param losers. Character vector or formula specifying the losers of each result
#' @param margin. Numeric vector vector or formula specifying the margin of victory, given as winner score - loser score
#' @param k.win. Numeric value of the learning rate to be applied to the win result
#' @param scale.margin. Numeric scaling factor applied in the expectation step for the MOV
#' @param scale.win. Numeric scaling factor applied in the expectation step for the win prediction
#' @param alpha. Numeric base rate for the logistic multiplication factor.
#' @param data. Data frame containing winner, loser, and margin variables if using a data/formula specification.
#' @param default. Numeric value of the initial rating to assign to new competitors
#'
#' @return A data frame with Elo ratings before and after each event result.
#'
#' @section Details:
#' Datasets should be ordered from first game result to last.
#' Competitors must be uniquely and consistently identified in the winner and loser vectors.
#' Missing values in the MOV variable will be omitted and will throw a warning. 
#' 
#' The E-step for the multiplicative model is the same as standard Elo:
#' \deqn{\hat{W} = \frac{1}{1+10^({R_j - R_i}{\sigma_{win}})}}.
#' The U-step for the multiplicative model involves updates based on the observed MOV. In terms of the \eqn{i}th player, 
#' \deqn{R_{i+1} = R_i + K_{win} (1 + |\hat{MOV}_{ij}/\sigma_{margin}|)^\alpha (W_{ij} - \hat{W}_{ij})}.
#' Where \eqn{W_{ij}} is a 0-1 indicator for the win result. The unknown parameters are \eqn{\sigma_{margin}},  \eqn{\sigma_{win}}, the power \eqn{\alpha} and learning rate \eqn{K_{win}}. Typical values for \eqn{K_{win}} and \eqn{\sigma_{win}} can be guided by the standard Elo system. For, \eqn{\sigma_{margin}}, reasonable choices would be the standard deviation (or some multiple) of the MOV variable. The \eqn{\alpha} modifies the rate of the MOV multiplier and choices between positive value between 0.2 and 2, for an SD-scaled MOV would ensure that the magnification for a single update was usually no more than double the margin difference in SD units. 
#' @examples
#' # Grand Slam MOV Elo Rating
#' ratings <- multiplicative(~ winner, ~loser, ~ game_margin, data = atp_games, alpha = 1, k.win = 24, scale.margin = 2, scale.win = 400)
#' @export
multiplicative <- function(winners, losers, margin, k.win, scale.margin, scale.win, alpha,  data, default = 1500) {		

	   ratings <- list()
	   
	   if(missing(data)){
		   winners <- as.character(winners)
		   losers <- as.character(losers)
		}
		else{
			v <- variables(winners, losers, margin, data)						
			winners <- as.character(v[[1]])
			losers <- as.character(v[[2]])
			margin <- v[[3]]
		}
	   
	   if(any(is.na(margin))){
	   	
	   		warning("Missing values in MOV found and will be excluded.")
	   	
	   		exclude <- is.na(margin)
	   		winners <- winners[!exclude]
	   		losers <- losers[!exclude]
	   		margin <- margin[!exclude]
	   }
	   
	   nlength <- length(winners)
		
	   results <- data.frame(
	   						winner = winners,
	   						loser = losers,
	   						winner_margin = margin,
	   						winner_before_elo = numeric(nlength),
	   						loser_before_elo = numeric(nlength),
	   						win_prediction = numeric(nlength),
	   						winner_elo = numeric(nlength), 
	              			loser_elo = numeric(nlength),
	              			stringsAsFactors = F
	              			)
              			
	
	   for (i in 1:nlength) {
		
	    cur_winner <- winners[i]
	    cur_loser <- losers[i]
	
	    winner_elo <- lookup(cur_winner, ratings, default)
	    loser_elo <- lookup(cur_loser, ratings, default)
	
	    winner_margin <- abs(margin[i] / scale.margin)
	    winner_prob <- 1/(1 + 10^(-1 * (winner_elo - loser_elo)/scale.win))
	    	    
	    winner_update <- k.win * (1 + winner_margin)^alpha * (1 - winner_prob)
	    loser_update <-  -winner_update
	    	     	   
	    ratings[[cur_winner]] = winner_elo + winner_update 
	    ratings[[cur_loser]] = loser_elo + loser_update 	   
	    results$win_prediction[i] <- winner_prob
	    results$winner_elo[i] <- ratings[[cur_winner]]
	    results$loser_elo[i] <- ratings[[cur_loser]]
	    results$winner_before_elo[i] <- winner_elo
	    results$loser_before_elo[i] <- loser_elo	    
	    
	  }
	  
results		
}