I am using MCMCpack to fit an age structured SIR model to seroprevalence data. I have the following code which includes an SIR model coded in Odin.
source("C:/Users/lra36/OneDrive - University of Cambridge/Documents/Varicella/Code/Odin/Model Run.R")
pop <- read.csv("C:/Users/lra36/OneDrive - University of Cambridge/Documents/Varicella/Code/population/Model_population.csv")$x
population <- pop[3:23]
observations <- read.csv("C:/Users/lra36/OneDrive - University of Cambridge/Documents/Varicella/Code/input CSVs/fitting data/UKHSA Seroprevalance.csv")$seropos
observations <- as.numeric(observations)/100
sigma <- read.csv("C:/Users/lra36/OneDrive - University of Cambridge/Documents/Varicella/Code/input CSVs/fitting data/UKHSA Seroprevalance.csv")$sd
sigma <- as.numeric(sigma)
B10 <- 0.95
B20 <- 0.99
B30 <- 0.15
params <- c(B10,B20,B30)
ModelSero <- function(params){
B1 <- exp(params[1])
B2 <- exp(params[2])
B3 <- exp(params[3])
Varicella_Final <- Varicella_Model$new(N_age=46,
contact_matrix = as.matrix(Contact_Rate),
n_births = births,
ageing = ageing,
beta = c(rep(B1,8), rep(B2,6), rep(B3,8), beta[23:46]),
rho = rho,
mort = mort,
Sv_ini =Sv0,
Ev_ini =Ev0,
Iv_ini =Iv0,
Sz_ini =Sz0,
Iz_ini =Iz0,
Rz_ini =Rz0,
B_ini =B0)
# Runs the model for 100 years and creates a vector of seroprevalence that has the same structure as the UKHSA data
model_results <- as.data.frame(Varicella_Final$run(100*365)[,])
age <- c(1,1.5,2:20)
model_results <- model_results[36501,]
S <- as.numeric(model_results[6:26])
E <- as.numeric(model_results[52:72])
I <- as.numeric(model_results[98:118])
seropos <- population - S - E - I
seroprev <- vector()
for (i in 1:21) {
seroprev[i]<- (seropos[i]/population[i])
}
return(seroprev)
}
prior_prob <- function(params){
#creates a beta vector with random values for the values I am trying to fit, keeps the others fixed.
B1 <- params[1]
B2 <- params[2]
B3 <- params[3]
dnorm(B1, mean=0, sd=1, log = TRUE) + dnorm(B2, mean=0, sd=1, log = TRUE) + dnorm(B3, mean=0, sd=1, log = TRUE)
}
logliklihood <- function(params, data){
predictions <- ModelSero(params)
log_likelihood <- vector()
print(params)
for (i in 1:21){
log_likelihood[i] = observations[i] * log(predictions[i]) + (1 - observations[i]) * log(1 - predictions[i])
}
return(sum(log_likelihood))
}
tic <- Sys.time() #store start time of simulation
print(Sys.time())
Post_Sample <- MCMCmetrop1R(fun = function(params){logliklihood(params,observations) + prior_prob(params)}, theta.init = log(c(B10,B20,B30)), mcmc = 27, burnin=0)
(Sys.time()-tic) #show duration of simulation run
beep(3)
I initially ran the model for 10 iterations and then tried to run it for 50 to get an idea of how long it was going to take. At this point I got the following error.
Error in MCMCmetrop1R(fun = function(params) { : `fun' returned NaN or NA
I understand that the chain is proposing parameter values that are resulting in NaN or NA values. I have done some investigating and found that after a certain amount of time the SIR function stops working and outputs INF and NA values. But now I don't know how to fix it?
I think I might have problems with my exp and logs but I can't see an error. Any ideas would be most welcome.