SIR models in brms

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

4 Likes

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")
2 Likes

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")
1 Like

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.

2 Likes

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.

1 Like

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.

1 Like

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]; 
}
2 Likes

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.

1 Like