{"id":2592,"date":"2022-07-23T19:56:12","date_gmt":"2022-07-23T19:56:12","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=2592"},"modified":"2022-11-08T03:34:15","modified_gmt":"2022-11-08T03:34:15","slug":"optimization-algorithms-in-r-returning-model-fit-metrics","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/optimization-algorithms-in-r-returning-model-fit-metrics\/","title":{"rendered":"Optimization Algorithms in R – returning model fit metrics"},"content":{"rendered":"

Introduction<\/strong><\/span><\/p>\n

A colleague had asked me if I knew of a way to obtain model fit metrics, such as AIC or r-squared, from the optim()<\/strong> function. First, optim()<\/strong> provides a general-purpose method of optimizing an algorithm to identify the best weights for either minimizing or maximizing whatever success metric you are comparing your model to (e.g., sum of squared error, maximum likelihood, etc.). From there, it continues until the model coefficients are optimal for the data.<\/p>\n

To make optim()<\/strong> work for us, we need to code the aspects of the model we are interested in optimizing (e.g., the regression coefficients) as well as code a function that calculates the output we are comparing the results to (e.g., sum of squared error).<\/p>\n

Before we get to model fit metrics, let’s walk through how optim()<\/strong> works by comparing our results to a simple linear regression. I’ll admit, optim()<\/strong> can be a little hard to wrap your head around (at least for me, it was), but building up a simple example can help us understand the power of this function and how we can use it later on down the road in more complicated analysis.<\/p>\n

Data<\/strong><\/span><\/p>\n

We will use data from the Lahman<\/strong> baseball data base. I’ll stick with all years from 2006 on and retain only players with a minimum of 250 at bats per season.<\/p>\n

\r\nlibrary(Lahman)\r\nlibrary(tidyverse)\r\n\r\ndata(Batting)\r\n\r\ndf <- Batting %>% \r\n  filter(yearID >= 2006,\r\n         AB >= 250)\r\n\r\nhead(df)\r\n<\/pre>\n
 <\/p>\n
Linear Regression<\/strong><\/span><\/p>\n
First, let’s just write a linear regression to predict HR from Hits, so that we have something to compare our optim()<\/strong> function against.<\/p>\n
\r\nfit_lm <- lm(HR ~ H, data = df)\r\nsummary(fit_lm)\r\n\r\nplot(df$H, df$HR, pch = 19, col = "light grey")\r\nabline(fit_lm, col = "red", lwd = 2)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
The above model is a simple linear regression model that fits a line of best fit based on the squared error of predictions to the actual values.<\/p>\n
(NOTE:<\/strong> We can see from the plot that this relationship is not really linear, but it will be okay for our purposes of discussion here.<\/em>)<\/p>\n
We can also use an optimizer to solve this.<\/p>\n
Optimizer<\/strong><\/span><\/p>\n
To write an optimizer, we need two functions:<\/p>\n
    \n
  1. \u00a0A function that allow us to model the relationship between HR and H and identify the optimal coefficients for the intercept and the slope that will help us minimize the residual between the actual and predicted values.<\/li>\n
  2. \u00a0A function that helps us keep track of the error between actual and predicted values as optim()<\/strong> runs through various iterations, attempting a number of possible coefficients for the slope and intercept.<\/li>\n<\/ol>\n
    Function 1: A linear function to predict HR from hits<\/strong><\/em><\/p>\n
    \r\nhr_from_h <- function(H, a, b){\r\n  return(a + b*H)\r\n}\r\n<\/pre>\n
    This simple function is just a linear equation and takes the values of our independent variable (H) and a value for the intercept (a) and slope (b). Although, we can plug in numbers and use the function right now, the values of a and b have not been optimized yet. Thus, the model will return a weird prediction for HR. Still, we can plug in some values to see how it works. For example, what if the person has 30 hits.<\/p>\n
    \r\nhr_from_h(H = 30, a = 1, b = 1)\r\n<\/pre>\n
    \"\"<\/a><\/p>\n
    Not a really good prediction of home runs!<\/p>\n
    Function 2: Sum of Squared Error Function<\/strong><\/em><\/p>\n
    Optimizers will try and identify the key weights in the function by either maximizing or minimizing some value. Since we are using a linear model, it makes sense to try and minimize the sum of the squared error.<\/p>\n
    The function will take 3 inputs:<\/p>\n
      \n
    1. \u00a0A data set of our dependent and independent variables<\/li>\n
    2. \u00a0A value for our intercept (a)<\/li>\n
    3. \u00a0A value for the slope of our model (b)<\/li>\n<\/ol>\n
      \r\nsum_sq_error <- function(df, a, b){\r\n  \r\n  # make predictions for HR\r\n  predictions <- with(df, hr_from_h(H, a, b))\r\n  \r\n  # get model errors\r\n  errors <- with(df, HR - predictions)\r\n  \r\n  # return the sum of squared errors\r\n  return(sum(errors^2))\r\n  \r\n}\r\n<\/pre>\n
       <\/p>\n
      Let’s make a fake data set of a few values of H and HR and then assign some weights for a<\/strong><\/em> and b<\/strong><\/em> to see if the sum_sq_error()<\/strong> produces a single sum of squared error value.<\/p>\n
      \r\nfake_df <- data.frame(H = c(35, 40, 55), \r\n                        HR = c(4, 10, 12))\r\n\r\nsum_sq_error(fake_df,\r\n             a = 3, \r\n             b = 1)\r\n<\/pre>\n
      \"\"<\/a><\/p>\n
      It worked! Now let’s write an optimization function to try and find the ideal weights for\u00a0a<\/strong><\/em> and b<\/strong><\/em> that minimize that sum of squared error. One way to do this is to create a large grid of values, write a for<\/strong> loop and let R plug along, trying each value, and then find the optimal values that minimize the sum of squared error. The issue with this is that if you have models with more independent variables, it will take really long. A more efficient way is to write an optimizer that can take care of this for us.<\/p>\n
      We will use the optim()<\/strong> function from base R.<\/p>\n
      The optim()<\/strong> function takes 2 inputs:<\/p>\n
        \n
      1. \u00a0A numeric vector of starting points for the parameters you are trying to optimize. These can be any values.<\/li>\n
      2. \u00a0A function that will receive the vector of starting points. This function will contain all<\/strong><\/em> of the parameters that we want to optimize. This function will take our sum_sq_error()<\/strong> function and it will get passed the starting values for a<\/strong> and b<\/strong> and then find their values that minimize the sum of squared error.<\/li>\n<\/ol>\n
        \r\noptimizer_results <- optim(par = c(0, 0),\r\n                           fn = function(x){\r\n                             sum_sq_error(df, x[1], x[2])\r\n                             }\r\n                           )\r\n<\/pre>\n
        Let’s have a look at the results.<\/p>\n
        \"\"<\/a><\/p>\n
        In the output:<\/p>\n