Maximum Likelihood Method for an ageing SIR model

Question

I am trying to use MLE to fit an age structured, deterministic SIR style model for varicella to data. Independently the model runs fine and gives feasible results. This is being used as a more robust fit to data. But... I get an error code.

var.data <- read.csv("/Users/laurenadams/Documents/gitrepo/Varicella/Input CSV/Age Specific Incidence.csv")

source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Equations/No Vax.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Functions.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Force of Infection.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Parameters.R")

#set initial conditions
M_inits     <- P[1]
Sv0_inits   <- P[2] - 0.000049*P[2] - 0.5*P[2] - 0.000012*P[2] - 0.20*P[2] - 0.20*P[2]      
Sv_inits    <- P_new[-1]-0.000049*P_new[2:100] - 0.5*P[2:100] - 0.000012*P_new[2:100] -      0.20*P_new[2:100] - 0.20*P_new[2:100]
Iv0_inits   <- 0.000049*P[2]
Iv_inits    <- 0.000049*P_new[2:100] 
Sz0_inits   <- 0.5*P[2]-0.000012*P[2]   
Sz_inits    <- 0.5*P[2:100]-0.000012*P_new[2:100]
Iz0_inits   <- 0.000012*P[2]   
Iz_inits    <- 0.000012*P_new[2:100] 
Rz0_inits   <- 0.40*P[2]
Rz_inits    <- 0.40*P_new[2:100] 
B10_inits   <- 0    
B1_inits    <- 0*P_new[2:100]
VInc0_inits <- 0
VInc_inits  <- c(rep(0,99))
ZInc0_inits <- 0
ZInc_inits  <- c(rep(0,99))

inits_all <- c(M=M_inits, Sv0=Sv0_inits, Sv=Sv_inits,
           Iv0=Iv0_inits,Iv=Iv_inits,
           Sz0 = Sz0_inits,Sz=Sz_inits,
           Iz0=Iz0_inits,Iz=Iz_inits,
           Rz0=Rz0_inits,Rz=Rz_inits,
           B10=B10_inits, B1=B1_inits,
           VInc0=VInc0_inits, VInc=VInc_inits,
           ZInc0=ZInc0_inits, ZInc=ZInc_inits)


# mle.sir - Maximum likelihood estimation function for closed SIR model
mle.sir <- function(b) {
  t <- seq(0, 365)    
  beta <- exp(b)

  results <- as.data.frame(ode(y= inits_all, 
                           times=t, 
                           varicella_fun_novax,
                           parms=c(parms)), row.names=F)
  Y <- results[365,]
  Y <- Y[103:202]
 nll <- -sum(dpois(x=var.data$Case_Proportion, lambda=Y, log=TRUE))
 return(nll)
 #return(results)
}

# initial - Initial estimates of beta and gamma
initial <- list(b=beta)
# fit0 - preliminary fit of model to data using the initial estimates
fit0 <- mle2(mle.sir, start=initial)

Here is the error code i get...

Error in dpois(x = var.data$Case_Proportion, lambda = Y, log = TRUE) : 
  Non-numeric argument to mathematical function

I believe it's a problem with how i'm trying to read in the lambda argument. But I can't think of another way to do it to make it numeric?

Can you please give us a [mcve]? – Ben Bolker Dec 12 '22 at 16:48 — Ben Bolker, Dec 12 '22 at 16:48
only if i share all my code? – Lauren Adams Dec 12 '22 at 18:32 — Lauren Adams, Dec 12 '22 at 18:32

score 0 · Answer 1 · answered Dec 12 '22 at 22:32

If you take a single row of a data frame it won't be a numeric vector and can't be handled by dpois().

dd <- data.frame(x = 1:3, y = 1:3, z = 1:3)
dpois(dd[3,], lambda = 1)

Error in dpois(dd[3, ], lambda = 1) : Non-numeric argument to mathematical function

It will work if you turn the data frame back into a matrix (or, in your case, refrain from making it into a data frame in the first place):

dpois(as.matrix(dd)[3, ], lambda = 1)

or use unlist() or as.numeric() to convert the list into a numeric vector:

dpois(as.numeric(dd[3,]), lambda = 1)

Maximum Likelihood Method for an ageing SIR model

1 Answers1

Linked