Data Fitting To ODE Models

[SAHEED]

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}

[aephidayatuloh]

@SAHEED

paras = c(beta = par[1]

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?

[SAHEED]

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

[aephidayatuloh]

@SAHEED

lfn2 = function(p,

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.

[SAHEED]

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.

[aephidayatuloh]

What is the remaining error?

[SAHEED]

Thanks, I was able to fix the other errors.