Data Fitting To ODE Models

Hi, I am new in the use of R. I am fitting to an ODE Model to get some estimated values, but I did not get it right. It was full of errors. I need assistance on this.

library("deSolve")
require(deSolve)

datAAA <- data.frame(
  time=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
         13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
         24, 25, 26, 27, 28, 29, 30),
  cases= c(214934, 307403, 417550, 541812, 674511, 808728, 940842,
           1065300, 1177623, 1271363, 1347177, 1406166, 1449357, 1479819,
           1500481, 1515892, 1527636, 1542107, 1558937, 1581336,
           1609292, 1638694, 1670713, 1707410, 1752498, 1797982,
           1841027, 1882445, 1922997, 1963044))

sir_epidemicmodnew = function(t, y, parms) {

  S = y[1]
  L = y[2]
  I = y[3]
  J = y[4]
  K = y[5]
  A = y[6]
  M = y[7]
  with(as.list(params), {
    dS <-   618880 - ((beta * (L + eta1*I + eta2*J + eta3*K + eta4*A) * S)/(S + L + I + J + K + A)) - mu*S
    dL <- ((epsilon * beta * (L + eta1*I + eta2*J + eta3*K + eta4*A) * S)/(S + L + I + J + K + A)) - (kappa+mu)*L + phi*K 
    dI <- (((1 - epsilon) * beta * (L + eta1*I + eta2*J + eta3*K + eta4*A) * S)/(S + L + I + J + K + A)) + (1 - omega)*kappa*L - (gamma + mu)*I
    dJ <- omega*kappa*L - (tau + mu)*J + gamma*I
    dK <- (1 - alpha)*tau*J - (phi + mu)*K
    dA <-  alpha*J -(mu+delta)*A
    dM <- J + K + A + mu*(J + K) + (mu + delta)*A
    res = c(dS, dL, dI, dJ, dK, dA, dM)
    list(res)
  })
}

lfn2 = function(p, M, eta2, eta3, eta4, epsilon, kappa, phi, omega, gamma, tau, alpha, mu, delta) {
  times = seq(1, 30, by = 1)
  start = c(S = 94997521, L=0, I = 21493,J= 0, K = 0, A=0, M=214934)
  paras = c(beta = par[1], eta1 = par[2], eta2=eta2, eta3=eta3, eta4=eta4, epsilon=epsilon, kappa=kappa, phi=phi, omega=omega, gamma=gamma, tau=tau, alpha=alpha, mu=mu, delta=delta)
  out = as.data.frame(ode(start, times = times,sir_epidemicmodnew, paras))
  n = length(M)
  rss = sum((M - out$M)^2)
  return(log(rss)*(n/2) - n*(log(n) - log(2*pi) - 1)/2)
}


paras0 = c(0.752, 1.5)
sirfit = optim(paras0,lfn2,M = datAAA$cases,eta2=1.4, eta3=1.3, eta4=0.015, epsilon=0.8, kappa=0.01, phi=0.01, omega=0.1, gamma=0.028475, tau=0.99, alpha=0.1, mu=0.022, delta=0.33,hessian = TRUE)
#> Error in par[1]: object of type 'closure' is not subsettable



times = seq(1, 30, by=.1)
start = c(S = 94997521, L=0, I = 21493,J= 0, K = 0, A=0, M=214934)
paras=c(beta=sirfit$par[1], eta1=sirfit$par[2], eta2=1.4, eta3=1.3, eta4=0.015, epsilon=0.8, kappa=0.01, phi=0.01, omega=0.1, gamma=0.028475, tau=0.99, alpha=0.1, mu=0.022, delta=0.33)
#> Error in eval(expr, envir, enclos): object 'sirfit' not found
out = as.data.frame(ode(start, times=times,
                        sir_epidemicmodnew, paras))
#> Error in as.list(params): object 'params' not found
plot(out$time, out$M, ylab="Prevalence",
     xlab="Day", type="l")
#> Error in plot(out$time, out$M, ylab = "Prevalence", xlab = "Day", type = "l"): object 'out' not found
points(datAAA$time, datAAA$cases)
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

Created on 2022-11-21 with reprex v2.0.2

When defining the lfn2 function, you call object named par for paras object. But there is no argument called par in lfn2 argument. There is only p argument. Are refering to this argument as par?

I have changed par to parms, the code still gave errors.

Why to parms? Argument you used in lfn2 are p, M, and so on. May be you need to change paras = c(beta = p[1], eta1 = p[2], ....) in your lfn2 function?

And FYI, par is a R's reserved name for a function. So,if you don't declare it as an argument of a function or an object you created before call it then par is refer to the R {graphics}'s par function. That's why you get error Error in par[1] : object of type 'closure' is not subsettable.

Yes, when I changed par[1] to p[1], the error concerning it went off, but I still have some error messages as before. The code did not work.

What is the remaining error?

Thanks, I was able to fix the other errors.