MCMC code comes into an error when increasing iterations

28 views Asked by At

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.

0

There are 0 answers