I’ve made an SIRD model using the deSolve package for the spread of COVID-19 in early 2020 in Mexico. It involves 4 parameters and I am now trying to estimate the values of these parameters using it’s predicted values of D and an optimising sums of squared distances to the true data approach. I tried making a function that takes parameter values from my model as an input and returns the sum of squared distances from the true data for these parameter values, I then tried defining a new function in terms of this first one that vectorises the parameters so I can use the optim() function on it. (The function “sird_1” runs the model for the parameters specified in it’s argument and I know it works as intended) This is my code:

   OWID <- read.csv("owid_covid_data.csv")
   DEATHS <- OWID[65675:65875 , 8]

   squaresum <- function(beta, sigma, theta, mu, data = DEATHS, N = 126000000){
   I0 <- DEATHS[1]
   times <- time_values 
   predictions <- sird_1(beta = beta, sigma = sigma, theta = theta, 
   mu = mu, S0 = N - I0, I0 = I0, R0 = 0, D0 = 0, times = time_values)
   sum((predictions$D[-1] - DEATHS[-1])^2, na.rm = T)
   
   }
   squaresum(beta = 0.160101, sigma = 0.1055556, theta = 0.00959596, mu = 0) 


  ss2 <- function(x) {
  squaresum(beta = x[1], sigma = x[2], theta = x[3], mu = x[4])
  }
  ss2(c(0.004, 0.5, 0.2, 0))

  starting_param_val <- c(0.5, 0.01, 0.04, 0)
  ss_optim <- optim(starting_param_val, ss2)

  ss_optim

  ss_optim$par

When I try and run this code it outputs [1]0 for any value of the squaresum() I try and then when I run the second half ss_optimjust returns my intial parameter values. The second error is clearly a consequence of the first however I don’t understand the error I have made when defining the squaresum() function and would really appreciate help.
(Note: ‘DEATHS’ contains N/A cells, I tried to ignore these using na.rm = T )