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?
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)