richfitz
commented
10 months ago You have out of date versions of one of the packages. There was quite a bit of movement in these earlier in the year with breaking changes as we finalise the interface. Please me sure you have the most recent versions of everything and it should work:
- dust 0.14.10
- odin.dust 0.3.7
- mcstate 0.9.17
- odin 1.5.6
Install it all with
install.packages(
c("odin", "odin.dust", "mcstate", "dust"),
repos = c("https://mrc-ide.r-universe.dev", "https://cloud.r-project.org"))Usually we update the minimum required versions when these incompatibilities arise, but R does enforce this in all situations.
whzhuscu
Maybe a silly question. I can call the function using the "odin.dust::odin_dust" statement in the R software, as shown in Figure 1. However, when running "mcstate::particle_filter$new" I encountered an error: "Error in model$public_methods$time_type() : attempt to apply non-function" as shown in Figure 2-3. I also tried using other computers, but the errors seemed to be more frequent, possibly due to the lack of C++ environment configuration, as shown in Figure 4-5.
But if I want to change the initial values of the SIR model or the dynamics of the transmission model (such as SEIR), it seems necessary to write the model myself and call it through "odin.dust::odin_dust". My code is provided below, and I would appreciate your guidance on how to make modifications.
` rm(list = ls()) library(mcstate) library(odin) library("odin.dust")
incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate")) incidence plot(cases ~ day, incidence, type = "o", xlab = "Day", ylab = "New cases", pch = 19)
history <- readRDS("C:/sir_true_history.rds") history <- t(drop(history)) colnames(history) <- c("S", "I", "R", "cases") history <- data.frame(t = seq_len(nrow(history)) - 1L, history) history matplot(history$t, history[c("S", "I", "R")], type = "l", lty = 1, xlab = "Day", ylab = "Variable")
sir <- odin.dust::odin_dust({ update(S) <- S - n_SI update(I) <- I + n_SI - n_IR update(R) <- R + n_IR
n_IR <- rbinom(I, 1 - exp(-beta * I / (S + I + R))) n_SI <- rbinom(S, 1 - exp(-gamma))
initial(S) <- S_ini initial(I) <- I_ini initial(R) <- 0
S_ini <- user(1000) I_ini <- user(10) beta <- user(0.2) gamma <- user(0.1) })
dt <- 0.25 data <- mcstate::particle_filter_data(incidence, time = "day", rate = 1 / dt)
compare <- function(state, observed, pars = NULL) { if (is.na(observed$cases)) { return(NULL) } exp_noise <- 1e6 incidence_modelled <- state[1, , drop = TRUE] incidence_observed <- observed$cases lambda <- incidence_modelled + rexp(n = length(incidence_modelled), rate = exp_noise) dpois(x = incidence_observed, lambda = lambda, log = TRUE) }
index <- function(info) { list(run = 5L, state = 1:5) }
proposal_kernel <- rbind(c(0.00057, 0.00034), c(0.00034, 0.00026)) pars <- mcstate::pmcmc_parameters$new( list(mcstate::pmcmc_parameter("beta", 0.2, min = 0, max = 1, prior = function(p) log(1e-10)), mcstate::pmcmc_parameter("gamma", 0.1, min = 0, max = 1, prior = function(p) log(1e-10))), proposal = proposal_kernel) pars
n_steps <- 500 n_particles <- 100 n_threads <- dust::dust_openmp_threads()
p1 <- mcstate::particle_filter$new(data, sir, n_particles, compare, index = index, n_threads = n_threads) control1 <- mcstate::pmcmc_control(n_steps, save_trajectories = TRUE, save_state = TRUE, save_restart = 40) res1 <- mcstate::pmcmc(pars, p1, control = control1)
t <- 0:100 matplot(t, t(res1$trajectories$state[2, , ]), type = "l", lty = 1, col = "#00000011", xlab = "Day", ylab = "Infected") points(I ~ t, history, pch = 19, col = "blue")
matplot(t, t(res1$trajectories$state[5, , ]), type = "l", lty = 1, col = "#00000005", xlab = "Day", ylab = "Daily incidence") points(cases ~ day, incidence, col = "blue", pch = 19) ` code.txt
