Hello,
I was wondering if it is possible to build SIR compartmental models like this one https://rstudio-pubs-static.s3.amazonaws.com/270496_e28d8aaa285042f2be0c24fc915a68b2.html with brms, in order to use brms infrastructure to simplify the prediction of the evolution of the epidemic.
If yes, is there any resource about it? I could not find anything about it.
Thank you very much
Angelo
For an example with differential equations (i.e., similar to the PK/PD resource), see here: https://jrmihalj.github.io/estimating-transmission-by-fitting-mechanistic-models-in-Stan/
At the risk of exposing my incompetence, I had luck with the following code, applied to a data.table adapted from the exercise in the link above:
myFuns <- "
real[] ode_SIR(real t, // time
real [] y, // state
real [] theta, // parameters
real [] x_r, // real data
int [] x_i) { // integer data
real dydt[3];
dydt[1] = - theta[1] * y[1] * y[2]; // dS/dt = -beta*I*S
dydt[2] = theta[1] * y[1] * y[2] - theta[2] * y[2]; // dI/dt = beta*S*I - gamma*I
dydt[3] = theta[2] * y[2]; // dR/dt = gamma*I
return dydt;
}
real sir_process(real t,
real beta,
real gamma,
real initfrac) {
real y[1,3];
real y0[3]; // 3 initial conditions
real theta[2];
y0[1] = 1 - initfrac;
y0[2] = initfrac;
y0[3] = 0;
theta[1] = beta;
theta[2] = gamma;
y = integrate_ode_rk45(ode_SIR,
y0,
0,
rep_array(t,1),
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return y[1,2];
}"
sir_brms <- brm(
bf(sample_y | trials(init_pop) ~ sir_process(time,
exp(beta),
exp(gamma),
inv_logit(initfrac)),
beta ~ 1,
gamma ~ 1,
initfrac ~ 1,
nl = TRUE),
stanvars = stanvar(scode = myFuns, block = "functions"),
data = sample_propinf,
family = binomial(link = "identity"),
prior = c(prior(normal(-2,0.5), nlpar = "beta"),
prior(normal(-2,0.5), nlpar = "gamma"),
prior(normal(0,4), nlpar = "initfrac")),
chains = 4,
cores = 4,
inits = 0,
iter = 1000,
warmup = 500)
Then to visualize posterior draws via brms, for example:
expose_functions(sir_brms,vectorize=TRUE)
pp_check(sir_brms,"intervals")I have been adjusting your code slightly and pasted together with the stuff from the block post. The code below compiles but then fails with the error message:
Error in compileCode(f, code, language = language, verbose = verbose) :
Compilation ERROR, function(s)/method(s) not created! file2b682db324c6.cpp:6:36: warning: ISO C99 requires whitespace after the macro name
#define STAN__SERVICES__COMMAND_HPP#include <boost/integer/integer_log2.hpp>
@bgoodri any idea what might be going on here?
The code to reproduce the error is the following:
# Load some required packages
##############
library(deSolve)
library(brms)
# rstan_options(auto_write = TRUE)
# options(mc.cores = parallel::detectCores())
##############
# To simulate the data, we need to assign initial conditions.
# In practice, these will likely be unknown, but can be estimated from the data.
I0 = 0.02 # initial fraction infected
S0 = 1 - I0 # initial fraction susceptible
R0 = 0
# Assign transmission and pathogen-induced death rates:
beta = 0.60
gamma = 0.10
# We will use the package deSolve to integrate, which requires certain data structures.
# Store parameters and initial values
# Parameters must be stored in a named list.
params <- list(beta = beta,
gamma = gamma)
# Initial conditions are stored in a vector
inits <- c(S0, I0, R0)
# Create a time series over which to integrate.
# Here we have an epidemic that is observed over t_max number of days (or weeks or etc).
t_min = 0
t_max = 50
times = t_min:t_max
# We must create a function for the system of ODEs.
# See the 'ode' function documentation for further insights.
SIR <- function(t, y, params) {
with(as.list(c(params, y)), {
dS = - beta * y[1] * y[2]
dI = beta * y[1] * y[2] - gamma * y[2]
dR = gamma * y[2]
res <- c(dS,dI,dR)
list(res)
})
}
# Run the integration:
out <- ode(inits, times, SIR, params, method="ode45")
# Store the output in a data frame:
out <- data.frame(out)
colnames(out) <- c("time", "S", "I", "R")
# quick plot of the epidemic
# plot(NA,NA, xlim = c(t_min, t_max), ylim=c(0, 1), xlab = "Time", ylab="Fraction of Host Population")
# lines(out$S ~ out$time, col="black")
# lines(out$I ~ out$time, col="red")
# legend(x = 30, y = 0.8, legend = c("Susceptible", "Infected"),
# col = c("black", "red"), lty = c(1, 1), bty="n")
sample_days = 20 # number of days sampled throughout the epidemic
sample_n = 25 # number of host individuals sampled per day
# Choose which days the samples were taken.
# Ideally this would be daily, but we all know that is difficult.
sample_time = sort(sample(1:t_max, sample_days, replace=F))
# Extract the "true" fraction of the population that is infected on each of the sampled days:
sample_propinf = out[out$time %in% sample_time, 3]
# Generate binomially distributed data.
# So, on each day we sample a given number of people (sample_n), and measure how many are infected.
# We expect binomially distributed error in this estimate, hence the random number generation.
sample_y = rbinom(sample_days, sample_n, sample_propinf)
dat <- data.frame(sample_y, sample_n, time = sample_time)
myFuns <- "
real[] ode_SIR(real t, // time
real [] y, // state
real [] theta, // parameters
real [] x_r, // real data
int [] x_i) { // integer data
real dydt[3];
dydt[1] = - theta[1] * y[1] * y[2]; // dS/dt = -beta*I*S
dydt[2] = theta[1] * y[1] * y[2] - theta[2] * y[2]; // dI/dt = beta*S*I - gamma*I
dydt[3] = theta[2] * y[2]; // dR/dt = gamma*I
return dydt;
}
real sir_process(real t,
real beta,
real gamma,
real initfrac) {
real y[1,3];
real y0[3]; // 3 initial conditions
real theta[2];
y0[1] = 1 - initfrac;
y0[2] = initfrac;
y0[3] = 0;
theta[1] = beta;
theta[2] = gamma;
y = integrate_ode_rk45(ode_SIR,
y0,
0,
rep_array(t,1),
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return y[1,2];
}
"
sir_brms <- brm(
bf(sample_y | trials(sample_n) ~
sir_process(time, exp(beta),exp(gamma), inv_logit(initfrac)),
beta ~ 1,
gamma ~ 1,
initfrac ~ 1,
nl = TRUE),
stanvars = stanvar(scode = myFuns, block = "functions"),
data = dat,
family = binomial(link = "identity"),
prior = c(prior(normal(-2,0.5), nlpar = "beta"),
prior(normal(-2,0.5), nlpar = "gamma"),
prior(normal(0,4), nlpar = "initfrac")),
chains = 2,
cores = 2,
inits = 0,
iter = 1000,
warmup = 500)
summary(sir_brms)
expose_functions(sir_brms,vectorize=TRUE)
pp_check(sir_brms,"intervals")
hi Paul,
Here’s the code I used to generate the data before calling brms. Slightly different from the original; in addition to the data.table package syntax, I also made sure to exclude t = 0 from the sampling. Not sure if that’s what’s behind the error. I also added init_pop (misnamed, it’s really the total sample size each day) to the data.table as a column. Should have noted this in my previous post.
# Load some required packages
##############
library(deSolve)
library(dplyr)
library(ggplot2)
library(rstan)
library(data.table)
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
##############
# To simulate the data, we need to assign initial conditions.
# In practice, these will likely be unknown, but can be estimated from the data.
I0 <- 0.02 # initial fraction infected
S0 <- 1 - I0 # initial fraction susceptible
R0 <- 0
# Assign transmission and pathogen-induced death rates:
beta <- 0.60
gamma <- 0.10
# We will use the package deSolve to integrate, which requires certain data structures.
# Store parameters and initial values
# Parameters must be stored in a named list.
params <- list(beta = beta,
gamma = gamma)
# Initial conditions are stored in a vector
inits <- c(S0, I0, R0)
# Create a time series over which to integrate.
# Here we have an epidemic that is observed over t_max number of days (or weeks or etc).
t_min <- 0
t_max <- 50
times <- t_min:t_max
# We must create a function for the system of ODEs.
# See the 'ode' function documentation for further insights.
SIR <- function(t, y, params) {
with(as.list(c(params, y)), {
dS <- - beta * y[1] * y[2]
dI <- beta * y[1] * y[2] - gamma * y[2]
dR <- gamma * y[2]
res <- c(dS,dI,dR)
list(res)
})
}
# Run the integration:
out <- ode(inits, times, SIR, params, method="ode45")
# Store the output in a data frame:
out <- as.data.table(out)
colnames(out) <- c("time", "S", "I", "R")
# quick plot of the epidemic
ggplot(out) +
geom_line(aes(x=time,y=S), color = "black") +
geom_line(aes(x=time,y=I), color = "red") +
geom_line(aes(x=time,y=R), color = "blue")
sample_days <- 20 # number of days sampled throughout the epidemic
sample_n <- 25 # number of host individuals sampled per day
# Choose which days the samples were taken.
# Ideally this would be daily, but we all know that is difficult.
sample_time <- out[time > 0,][sample(.N,sample_days)][order(time)][,time]
# Extract the "true" fraction of the population that is infected on each of the sampled days:
sample_propinf <- out[time %in% sample_time,]
# Generate binomially distributed data.
# So, on each day we sample a given number of people (sample_n), and measure how many are infected.
# We expect binomially distributed error in this estimate, hence the random number generation.
sample_propinf[,sample_y := rbinom(sample_days,sample_n,I)]
# Attempt brms model
sample_propinf[,init_pop := as.integer(sample_n)]
myFuns <- "
real[] ode_SIR(real t, // time
real [] y, // state
real [] theta, // parameters
real [] x_r, // real data
int [] x_i) { // integer data
real dydt[3];
dydt[1] = - theta[1] * y[1] * y[2]; // dS/dt = -beta*I*S
dydt[2] = theta[1] * y[1] * y[2] - theta[2] * y[2]; // dI/dt = beta*S*I - gamma*I
dydt[3] = theta[2] * y[2]; // dR/dt = gamma*I
return dydt;
}
real sir_process(real t,
real beta,
real gamma,
real initfrac) {
real y[1,3];
real y0[3]; // 3 initial conditions
real theta[2];
y0[1] = 1 - initfrac;
y0[2] = initfrac;
y0[3] = 0;
theta[1] = beta;
theta[2] = gamma;
y = integrate_ode_rk45(ode_SIR,
y0,
0,
rep_array(t,1),
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return y[1,2];
}"
sir_brms <- brm(
bf(sample_y | trials(init_pop) ~ sir_process(time,
exp(beta),
exp(gamma),
inv_logit(initfrac)),
beta ~ 1,
gamma ~ 1,
initfrac ~ 1,
nl = TRUE),
stanvars = stanvar(scode = myFuns, block = "functions"),
data = sample_propinf,
family = binomial(link = "identity"),
prior = c(prior(normal(-2,0.5), nlpar = "beta"),
prior(normal(-2,0.5), nlpar = "gamma"),
prior(normal(0,4), nlpar = "initfrac")),
chains = 4,
cores = 4,
inits = 0,
iter = 1000,
warmup = 500)
expose_functions(sir_brms,vectorize=TRUE)
pp_check(sir_brms,"intervals")
Hope this helps.
This works. Thank you. Not sure what was causing the error, but the error message what at least not very helpful :-D
As a small additon to your code conditional_effects(sir_brms, effects = "time") can also be used to conveniently visualize predictions.
If you call make_stancode and compile it with verbose = TRUE, you eventually get
file73bd1013dd8a.cpp:154:49: error: no matching function for call to ‘integrate_ode_rk45(model73bd1c443702_SIR_namespace::ode_SIR_functor__, std::vector<stan::math::var>&, int, std::vector<int>, std::vector<stan::math::var>&, std::vector<double>, std::vector<int>, std::ostream*&, double, double, double)’
stan::math::assign(y, integrate_ode_rk45(ode_SIR_functor__(), y0, 0, rep_array(t, 1), theta, rep_array(0.0, 0), rep_array(1, 1), pstream__, 1.0E-10, 1.0E-10, 1.0E3));
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr.hpp:46,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/mat.hpp:344,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:12,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math.hpp:4,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file73bd1013dd8a.cpp:8:
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: candidate: ‘template<class F, class T1, class T2> std::vector<std::vector<typename stan::return_type<T_shared_param, T_job_param>::type> > stan::math::integrate_ode_rk45(const F&, const std::vector<T>&, double, const std::vector<double>&, const std::vector<T_l>&, const std::vector<double>&, const std::vector<int>&, std::ostream*, double, double, int)’
integrate_ode_rk45(const F& f, const std::vector<T1>& y0, double t0,
^~~~~~~~~~~~~~~~~~
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: template argument deduction/substitution failed:
file73bd1013dd8a.cpp:154:49: note: cannot convert ‘stan::math::rep_array(const T&, int) [with T = int](1)’ (type ‘std::vector<int>’) to type ‘const std::vector<double>&’
stan::math::assign(y, integrate_ode_rk45(ode_SIR_functor__(), y0, 0, rep_array(t, 1), theta, rep_array(0.0, 0), rep_array(1, 1), pstream__, 1.0E-10, 1.0E-10, 1.0E3));
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
file73bd1013dd8a.cpp: In instantiation of ‘typename boost::math::tools::promote_args<T1, T2, T3, T4>::type model73bd1c443702_SIR_namespace::sir_process(const T0__&, const T1__&, const T2__&, const T3__&, std::ostream*) [with T0__ = int; T1__ = double; T2__ = double; T3__ = double; typename boost::math::tools::promote_args<T1, T2, T3, T4>::type = double; std::ostream = std::basic_ostream<char>]’:
file73bd1013dd8a.cpp:528:40: required from ‘T__ model73bd1c443702_SIR_namespace::model73bd1c443702_SIR::log_prob(std::vector<T_l>&, std::vector<int>&, std::ostream*) const [with bool propto__ = false; bool jacobian__ = false; T__ = double; std::ostream = std::basic_ostream<char>]’
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/src/stan/services/optimize/newton.hpp:61:14: required from ‘int stan::services::optimize::newton(Model&, stan::io::var_context&, unsigned int, unsigned int, double, int, bool, stan::callbacks::interrupt&, stan::callbacks::logger&, stan::callbacks::writer&, stan::callbacks::writer&) [with Model = model73bd1c443702_SIR_namespace::model73bd1c443702_SIR]’
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/rstan/include/rstan/stan_fit.hpp:488:41: required from ‘int rstan::{anonymous}::command(rstan::stan_args&, Model&, Rcpp::List&, const std::vector<long unsigned int>&, const std::vector<std::__cxx11::basic_string<char> >&, RNG_t&) [with Model = model73bd1c443702_SIR_namespace::model73bd1c443702_SIR; RNG_t = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014, 0, 2147483563>, boost::random::linear_congruential_engine<unsigned int, 40692, 0, 2147483399> >; Rcpp::List = Rcpp::Vector<19>]’
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/rstan/include/rstan/stan_fit.hpp:1201:18: required from ‘SEXPREC* rstan::stan_fit<Model, RNG_t>::call_sampler(SEXP) [with Model = model73bd1c443702_SIR_namespace::model73bd1c443702_SIR; RNG_t = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014, 0, 2147483563>, boost::random::linear_congruential_engine<unsigned int, 40692, 0, 2147483399> >; SEXP = SEXPREC*]’
file73bd1013dd8a.cpp:746:114: required from here
file73bd1013dd8a.cpp:154:49: error: no matching function for call to ‘integrate_ode_rk45(model73bd1c443702_SIR_namespace::ode_SIR_functor__, std::vector<double>&, int, std::vector<int>, std::vector<double>&, std::vector<double>, std::vector<int>, std::ostream*&, double, double, double)’
In file included from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr.hpp:46,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/mat.hpp:344,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:12,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math.hpp:4,
from /home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file73bd1013dd8a.cpp:8:
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: candidate: ‘template<class F, class T1, class T2> std::vector<std::vector<typename stan::return_type<T_shared_param, T_job_param>::type> > stan::math::integrate_ode_rk45(const F&, const std::vector<T>&, double, const std::vector<double>&, const std::vector<T_l>&, const std::vector<double>&, const std::vector<int>&, std::ostream*, double, double, int)’
integrate_ode_rk45(const F& f, const std::vector<T1>& y0, double t0,
^~~~~~~~~~~~~~~~~~
/home/ben/R/x86_64-pc-linux-gnu-library/3.5/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:1: note: template argument deduction/substitution failed:
file73bd1013dd8a.cpp:154:49: note: cannot convert ‘stan::math::rep_array(const T&, int) [with T = int](1)’ (type ‘std::vector<int>’) to type ‘const std::vector<double>&’
stan::math::assign(y, integrate_ode_rk45(ode_SIR_functor__(), y0, 0, rep_array(t, 1), theta, rep_array(0.0, 0), rep_array(1, 1), pstream__, 1.0E-10, 1.0E-10, 1.0E3));
~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If you change the declaration of C1 in the data block from
int C_1[N];
to
real C_1[N];
or
vector[N] C_1;
it compiles. I think what you have should auto-promote integers to real numbers but it does not for some reason, which is probably a Stan bug.
Or you can hack it by adding a multiplication by 1.0 in
real sir_process(real t,
real beta,
real gamma,
real initfrac) {
real y[1,3];
real y0[3]; // 3 initial conditions
real theta[2];
y0[1] = 1 - initfrac;
y0[2] = initfrac;
y0[3] = 0;
theta[1] = beta;
theta[2] = gamma;
y = integrate_ode_rk45(ode_SIR,
y0,
0,
rep_array(t * 1.0,1), // hack
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return y[1,2];
}
Thanks! brms has the feature of treating integer columns as integer to increase flexibility of non-linear modeling. Good to know what was causing the problem.
Hopefully a question for @paul.buerkner re: extending this model beyond SIR. Suppose we wanted to add an exposed component:
myFuns <- "
real[] ode_SIR(real t, // time
real [] y, // state
real [] theta, // parameters
real [] x_r, // real data
int [] x_i) { // integer data
real dydt[4];
dydt[1] = - theta[1] * y[1] * y[3]; // dS/dt = -beta*I*S
dydt[2] = theta[1] * y[1] * y[3] - theta[2] * y[2]; // dE/dt = beta*S*I - alpha*E
dydt[3] = theta[2] * y[2] - theta[3]*y[2]; // dI/dt = alpha*E - gamma*I
dydt[4] = theta[3] * y[2]; // dR/dt = gamma*I
return dydt;
}
real sir_process(real t,
real beta,
real alpha,
real gamma,
real initfrac,
real initimmune) {
real y[1,4];
real y0[4]; // 4 initial conditions
real theta[3]; // 3 parameters
y0[1] = 1 - (initfrac+initimmune);
y0[2] = initfrac;
y0[3] = 0;
y0[4] = initimmune;
theta[1] = beta;
theta[2] = alpha;
theta[3] = gamma;
y = integrate_ode_rk45(ode_SIR,
y0,
0,
rep_array(t,1),
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return (y[1,3] + y[1,2]);
}"
… and then the brms code:
seir_brms <- brm(bf(sample_y | trials(sample_size) ~ y,
nlf(y ~ inv_logit(Se)*p + (1-p)*(1 - inv_logit(Sp))),
nlf(p ~ sir_process(time,
exp(beta),
exp(alpha),
exp(gamma),
inv_logit(initfrac),
inv_logit(initimmune))),
Se ~ 1,
Sp ~ 1,
beta ~ 1,
alpha ~ 1,
gamma ~ 1,
initfrac ~ 1,
initimmune ~ 1,
nl= TRUE),
stanvars = stanvar(scode = myFuns, block = "functions"),
data = sampled_data[time <= 7,],
family = binomial(link = "identity"),
prior = c(prior(normal(-2,0.15), nlpar = "beta"),
prior(normal(-1.09,0.1), nlpar = "alpha"),
prior(normal(-2.65,0.1), nlpar = "gamma"),
prior(normal(-6,0.5), nlpar = "initfrac"),
prior(normal(-4,0.4), nlpar = "initimmune"),
prior(normal(5,0.3), nlpar = "Se"),
prior(normal(5,0.3), nlpar = "Sp")),
chains = 4,
cores = 4,
inits = 0,
iter = 2000,
warmup = 1000,
control = list(adapt_delta = 0.9))
The code above currently has the function returning (I + E) from the sir_process function (i.e., assuming positive tests indicate Infectious or Exposed), but how would I return them separately as nonlinear parameters in brms, as outputs from the same function? I.e., have something “like”
nlf(c(i,e) ~ sir_process(time,
exp(beta),
exp(alpha),
exp(gamma),
inv_logit(initfrac),
inv_logit(initimmune))),
My end goal is to extend this to multiple compartments of I and E, where each compartment might have some different level of viral shedding detectable in a non-binary test.
Thanks much for any assistance you can provide.
This is a functionality brms is not desigend to achieve. You have to use raw stan for this purpose.
I would like to leverage this discussion to guide me. I am also trying to fit my data to a ODE system. Summarily, I sampled Leishmania flies in a Brazilian city for many years (1998 to 2023). 16 locations sampled irregularly through time. Data for one place looks like this:
time=c(2,3,6,7,8,20,21) #time in years
N=c(1.6,0,77,68,47,69,92) #Proportion of Leishmania flies (leishmania flies/all captured flies [includes other species]) (zero to 100)
U=c(24.63,27.19,26.23,25.80,26.01,43.82,46.16) #Degree of urbanization (zero to 100)
My hypothesis is as time goes, sites become more urbanized, increasing the carrying capacity for leishmania flies (strongly associated with cities), which in turn make its population increases. So, I built this ODE system
dU <- ulfa * U * (1 - U / 100)
K <- 100 * (U / (U + effect))
dN <- r * N * (1 - N / K);
note parameters
ulfa depicts the rate of urbanization through time (maximum 100% o urbanization)
effect depicts the shape how urbanization increasing affects carrying capacity K (maximum also 100%, that means, all captured flies are from the leishmania species)
r depicts the rate of population growth for leishmania flies
Therefore, I would like to solve this using the stan/brms approach, which allows me to fit a unique model for my 16 sites. I can do it easily for one site using ode and optim fucntoin in R. Below follow my commented code. I appreciate any advice.
# Load required packages
library(deSolve)
library(dplyr)
library(ggplot2)
library(rstan)
library(data.table)
library(brms)
# Set initial conditions and parameters
U0 <- 10 # initial value of U for Urbanization (0 to 100)
N0 <- 5 # initial value of N for fly abundance (0 to 100)
ulfa <- 0.01 # initial value of ulfa
effect <- 5 # initial value of effect
r <- 0.1 # initial value of r
# Create parameter and initial condition vectors
params <- list(ulfa = ulfa, effect = effect, r = r)
inits <- c(U = U0, N = N0)
# Define the time span for integration
t_min <- 0
t_max <- 50
times <- t_min:t_max
# Define the ODE function
concept <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
dU <- ulfa * U * (1 - U / 100)
K <- 100 * (U / (U + effect))
dN <- r * N * (1 - N / K)
list(c(dU, dN))
})
}
# Run the integration
out <- ode(inits, times, concept, params, method = "ode45")
plot(out)
# Stoore the output in a data frame
out <- as.data.table(out)
colnames(out) <- c("time", "U", "N")
# Quick plot of the results
ggplot(out) +
geom_line(aes(x = time, y = U), color = "black") +
geom_line(aes(x = time, y = N), color = "red")
# Define the Stan code for the ODE system
myFuns <- "
real[] ode_concept(real t, // time
real[] y, // state
real[] theta, // parameters
real[] x_r, // real data
int[] x_i) { // integer data
real dydt[2];
dydt[1] = theta[1] * y[1] * (1 - y[1] / 100);
real K = 100 * (y[1] / (y[1] + theta[2]));
dydt[2] = theta[3] * y[2] * (1 - y[2] / K);
return dydt;
}
real concept_process(real t,
real ulfa,
real effect,
real r,
real initU,
real initN) {
real y[1,2];
real y0[2];
real theta[3];
y0[1] = initU;
y0[2] = initN;
theta[1] = ulfa;
theta[2] = effect;
theta[3] = r;
y = integrate_ode_rk45(ode_concept,
y0,
0,
rep_array(t,1),
theta,
rep_array(0.0,0),
rep_array(1,1),
1.0E-10, 1.0E-10, 1.0E3);
return y[1,2];
}"
# Sample data from simulated data using initial parameters
sample=out[out$time %in% sample(out$time,10),]
# Real data
time=c(2,3,6,7,8,20,21)
N=c(1.6,0,77,68,47,69,92)
U=c(24.63,27.19,26.23,25.80,26.01,43.82,46.16)
real=data.frame(time,N,U)
# Attempt brms model
sir_brms <- brm(
bf(N ~ concept_process(time,ulfa,effect,r,U,N),
ulfa ~ 1,
effect ~ 1,
r ~ 1,
nl = TRUE),
stanvars = stanvar(scode = myFuns, block = "functions"),
data = sample,
family = gaussian(link = "identity"),
prior = c(prior(normal(0, 1), nlpar = "ulfa"),
prior(normal(0, 10), nlpar = "effect"),
prior(normal(0, 1), nlpar = "r")),
chains = 4,
cores = 4,
init = 0,
iter = 1000,
warmup = 500
)
summary(sir_brms)
expose_functions(sir_brms, vectorize = TRUE)
conditional_effects(sir_brms,effects="time")
pp_check(sir_brms, "intervals")
Edited by @jsocolar for syntax highlighting.