Discrete Time Transmission Models in STAN
Arie Voorman
April 2017
Simple SIR Model
Lets start with a simple SIR model. Suppose we observe whether someone is Susceptable (S), Infectious (I), or Recovered (R) for individuals \(j=1 \dots n\) at times \(t=1 \dots T\). Denote \(S_t\), \(I_t\), and \(R_t\) as the set of individuals in each state at time \(t\). We need models for the transmission of a person \(S\rightarrow I\), and \(I \rightarrow R\).
Starting with the transition \(S\rightarrow I\), let the \(\lambda\) be the probability that a single infectious person in \(I\) infects a susceptable person over 1 unit of time. If the population is well-mixed, then the probability of a susceptable person becoming infected at time \(t\) given infectios status of the community at \(I_{t-1}\) is \[Pr(j \in I_t \,|\,\, j \in S_{t-1},\, I_{t-1} ) = 1- (1- \lambda)^{|I_{t-1}|},\] where \(|I_{t}|\) is the number of people who are infected at a particular time. The expression on the r.h.s. is just 1 - probability of not being infected from each person in \(I_{t-1}\). NB: This is an exact model for small population, though I believe some models will use the approximation $(1-)^{|I_{t-1}|} e^{-|I_{t-1}|} $ for large-ish populations.
Now, for the transition \(I\rightarrow R\), lets assume that mean shedding duration is \(\gamma\) units of time, and that shedding duration follows a geometric distribution (i.e. the probability an individual stops shedding over 1-unit of time is \(1/\gamma\))
Simulation example
The code below simulates data from this model. After the simulation, I’ll write up a model in STAN that estimates the parameters \(\lambda\) and \(\gamma\) from observed data, and then use them to estimate \(R_0\).
##################### simulate data ###################
set.seed(42)
N = 300 # number of subjects
t = 200 # number of time points
lambda = 0.001 # transmission probability (per person per 'close contact', per day)
I0 = 25 # number of children initially infected
gamma = 7 # mean shedding duration for an individual
#data frame to store results
data = expand.grid(subject = 1:N, time = 0:t) %>%
tbl_df()
#wdata is the 'wide' version of the data frame
wdata = data %>% mutate(shedding = NA_character_) %>%
spread(subject,shedding) %>% setNames(names(.) %>% make.names)
# keep track of immune and naive individuals
infected = c(rep(1,I0),rep(0,N-I0)) %>% as.numeric
susceptible = 1-infected
recovered = rep(0,N)
wdata[1,-1][susceptible ==1] = 'S'
wdata[1,-1][infected == 1] = 'I'
wdata[1,-1][recovered == 1] = 'R'
#run through the simulation:
for(i in 1:t){
#find new infected individuals
#note: multiplying by susceptible zeros out all those who are not susceptible
prob.infected = (1-(1-lambda)^sum(infected))*susceptible
new.infected = rbinom(N,1,prob.infected)
#probability of recovering 1= 1/gamma
new.recovered = rbinom(N,1,1/gamma)
new.recovered = new.recovered*infected
#update SIR
susceptible = if_else(new.infected==1,0,susceptible)
infected = if_else(new.infected==1,1,infected)
infected = if_else(new.recovered==1,0,infected)
recovered = if_else(new.recovered==1,1,recovered)
wdata[i+1,-1][susceptible ==1] = 'S'
wdata[i+1,-1][infected == 1] = 'I'
wdata[i+1,-1][recovered == 1] = 'R'
}
data = gather(wdata,subject,shedding, -time) %>% tbl_df()
data = data %>% group_by(subject) %>% mutate(start_shed = as.numeric(min(time[shedding == 'I']))) %>%
ungroup %>%
mutate(subject = reorder(subject,-start_shed), shedding = factor(shedding, levels = c('S','I','R'),ordered = T))
#plot the SIR status for each individual over time:
ggplot(data,aes(x=time,y=as.numeric(subject),fill = shedding)) +
geom_tile() +
scale_fill_manual('Status',values = c('S'='green3','I'='indianred','R'='blue3')) +
ylab('Subject') + xlab('Time')
Fitting the model in STAN
The code below specifies the model in (STAN)[http://mc-stan.org/] using the ‘rstan’ library. For the simple SIR case, I’ll use the fact that \(S \rightarrow I\) transition are identical bernoulli random variables, and can be aggregated to a binomial distribution.
I’ll also use the probability model to estimate \(R_0\) by simulating SIR model with a single infected individual. This simulation uses the estimated parameters an naturally propogates their uncertainty to give a confidence interval.
sircode =
"
data {
//declare variables
int<lower=0> N; //number of individuals
int<lower=0> T; //number of time points
int<lower = 0, upper =N> S[T]; //SIR
int<lower = 0, upper =N> I[T];
int<lower = 0, upper =N> R[T];
real<lower=0> lambda_max; //maximum value of lambda (sometimes useful to help STAN converge)
}
transformed data {
//Calculate transitions, based on SIR status
int<lower=0, upper = N> z_si[T-1];
int<lower=0, upper = N> z_ir[T-1];
for(i in 1:(T-1)){
z_si[i] = S[i] - S[i+1];
z_ir[i] = R[i+1] - R[i];
}
}
parameters {
real<lower=1> gamma;
real<lower=0, upper = lambda_max> lambda;
}
model {
lambda ~ uniform(0,lambda_max);
for(i in 1:(T-1)){
if(I[i] > 0){ //only define z_si when there are infections - otherwise distribution is degenerate and STAN has trouble
z_si[i] ~ binomial(S[i], 1-(1-lambda)^I[i]);
}
z_ir[i] ~ binomial(I[i], 1/gamma);
}
}
generated quantities {
int<lower=0> r_0; // simulate a single infected individual in the population, and count secondary infections (i.e. R_0).
r_0 = 0;
while(1){
r_0 = r_0 + binomial_rng(N-r_0,lambda);
if(bernoulli_rng(1/gamma)) break;
}
}
"
mod = stan_model(model_name = 'sir',model_code = sircode)## In file included from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,
## from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file2da052929dd.cpp:8:
## C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"# format data for STAN
data_stan = data %>% group_by(time) %>%
summarise(S = sum(shedding =='S') %>% as.integer,
I = sum(shedding =='I') %>% as.integer,
R = sum(shedding =='R') %>% as.integer) %>% select(-time) %>%
as.list %>%
c(N = as.integer(N),T = as.integer(t)+1,lambda_max = 1, .)
fit = sampling(mod,data_stan, iter = 1000, chains = 2)##
## SAMPLING FOR MODEL 'sir' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%] (Sampling)
## Elapsed Time: 0.082 seconds (Warm-up)
## 0.068 seconds (Sampling)
## 0.15 seconds (Total)
##
##
## SAMPLING FOR MODEL 'sir' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2, Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 2, Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 2, Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 2, Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 2, Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 2, Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2, Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 2, Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 2, Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 2, Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 2, Iteration: 1000 / 1000 [100%] (Sampling)
## Elapsed Time: 0.074 seconds (Warm-up)
## 0.069 seconds (Sampling)
## 0.143 seconds (Total)stan_dens(fit)
print(fit,digits = 4)## Inference for Stan model: sir.
## 2 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=1000.
##
## mean se_mean sd 2.5% 25% 50%
## gamma 6.7349 0.0147 0.4047 6.0215 6.4563 6.7196
## lambda 0.0010 0.0000 0.0001 0.0008 0.0009 0.0010
## r_0 2.0260 0.0759 2.3140 0.0000 0.0000 1.0000
## lp__ -1698.3973 0.0488 1.0041 -1701.4818 -1698.7709 -1698.1042
## 75% 97.5% n_eff Rhat
## gamma 6.977 7.5694 759 0.9997
## lambda 0.001 0.0011 1000 0.9987
## r_0 3.000 8.0000 929 1.0026
## lp__ -1697.665 -1697.4300 424 0.9981
##
## Samples were drawn using NUTS(diag_e) at Tue May 09 13:04:32 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).So, we see that the estimate of \(\lambda\) is 0.0010, with 95% credible interval of (0.0008, 0.0011), encompasing the true value of 0.001. Likewise, the CI for \(\gamma\) is (6,7.5), encolsing the true value of 7.
Finally, \(R_0\), which we did not know in advance is estimated at 1.9. Interpreting the confidence intervals is a little tricky. Since \(R_0\) comes from simulation, the above CI represents monte-carlo error, rather than uncertainty in parameters. One could get around this either by repeating the simulation for each draw from the posterior, which may be easiest for more complex models, or by explicitely calculated the expectation of \(R_0\).
SIR with Covariates
Making the model a little more complicated, lets assume that instead of a well-mixed population, there are children and adults in families. In the simulation below, the probability of getting infected from someone within a family is much higher than getting infected from someone in the community. In addition, the duration of shedding is much shorter for adults than for children.
Here, we use the probability model for infection: \[Pr(j \in I_t \,|\,\, j \in S_{t-1},\, I_{t-1} ) = 1- (1- \lambda_F)^{|I^F_{t-1}|}(1- \lambda_C)^{|I^C_{t-1}|},\] where \(\lambda_F\) and \(\lambda_C\) are transmission rates within families and the community, respectively, and \(I^F_t\) and \(I^C_t\) are the set of infected individuals in the family and community at time \(t\), respectively. While I’m not using the large-sample approximation \((1- \lambda_F)^{|I^F_{t-1}|}(1- \lambda_C)^{|I^C_{t-1}|} \approx e^{-\lambda_F|I^F_{t-1}| - \lambda_C |I^C_{t-1}|}\), it makes the probability log-linear, and is conducive to adding other covariates as well.
Likewise, we also now include \(\gamma_{adult}\) and \(\gamma_{child}\), which are just as in the model above.
set.seed(42)
N = 300 # number of subjects
t = 200 # number of time points
lambda_family = 0.1 # transmission probability (per person per 'close contact', per day)
lambda_community = 0.0005# transmission probability (per person per 'person in community', per day)
I0 = 25 # number of children initially infected
gamma_child = 10 # mean shedding duration for child (naive child)
gamma_adult = 3 # mean shedding duration for an adult (prior immunity)
fmat = bdiag(replicate(N/4,matrix(1,4,4),simplify = F)) #matrix indicating family membership
diag(fmat) = 0
adult = rep(c(1,1,0,0),N/4) # two adults and two children per HH
data = expand.grid(subject = 1:N, time = 0:t) %>%
tbl_df()
wdata = data %>% mutate(shedding = NA_character_) %>%
spread(subject,shedding) %>% setNames(names(.) %>% make.names)
# keep track of immune and naive individuals
infected = c(rep(c(0,0,0,1),I0), rep(0,N-I0*4)) %>% as.numeric
susceptible = 1-infected
recovered = rep(0,N)
wdata[1,-1][susceptible ==1] = 'S'
wdata[1,-1][infected == 1] = 'I'
wdata[1,-1][recovered == 1] = 'R'
for(i in 1:t){
#probability of infection = 1- probability of not getting infected, from either within HH, or within community
#set to zero for those that aren't susceptible
prob.infection = 1-((1-lambda_family)^(fmat%*%infected))*((1-lambda_community)^(sum(infected) - fmat%*%infected))
prob.infection = as.numeric(prob.infection*susceptible)
new.infection = rbinom(N,1,prob.infection)
#probability of recovering 1= 1/gamma
new.recovered = rbinom(N,1,adult*(1/gamma_adult) + (1-adult)*(1/gamma_child))
new.recovered = new.recovered*infected
#update SIR
susceptible = if_else(new.infection==1,0,susceptible)
recovered = if_else(new.recovered==1,1,recovered)
infected = if_else(new.infection==1,1,infected)
infected = if_else(new.recovered==1,0,infected)
wdata[i+1,-1][susceptible ==1] = 'S'
wdata[i+1,-1][infected == 1] = 'I'
wdata[i+1,-1][recovered == 1] = 'R'
}
data = gather(wdata,subject,shedding, -time) %>% tbl_df()
data = data %>% group_by(subject) %>% mutate(start_shed = as.numeric(min(time[shedding == 'I']))) %>%
ungroup %>%
mutate(subject = reorder(subject,-start_shed), shedding = factor(shedding, levels = c('S','I','R'),ordered = T))
ggplot(data,aes(x=time,y=as.numeric(subject),fill = shedding)) +
geom_tile() +
scale_fill_manual('Status',values = c('S'='green3','I'='indianred','R'='blue3')) +
ylab('Subject') + xlab('Time')
You can see here that the time of shedding is much more heterogeneous. Some individuals are infected for long periods, while others are infected for relatively short periods.
Fitting the model in STAN:
Fitting the model is much the same in STAN, with the notable exception that we have to model individuals as independent non-identical bernoulli random variables, rather than collapsing them into a Binomial random variable.
sircode_v2 = '
data {
int<lower=0> N;
int<lower=0> T;
int<lower = 0,upper=1> S[T,N];
int<lower = 0,upper=1> I[T,N];
int<lower = 0,upper=1> R[T,N];
matrix[N,N] fmat;
int<lower=0,upper=1> adult[N];
real<lower=0> lambda_max;
}
transformed data {
int<lower=0, upper = 1> z_si[T-1,N];
int<lower=0, upper = 1> z_ir[T-1,N];
matrix[T,N] I_mat;
matrix[T,N] fam_I;
matrix[T,N] com_I;
I_mat = to_matrix(I);
fam_I = I_mat*fmat;
for(t in 1:(T-1)){
for(i in 1:N){
z_si[t,i] = S[t,i] - S[t+1,i];
z_ir[t,i] = R[t+1,i] - R[t,i];
com_I[t,i] = sum(I_mat[t]) - I_mat[t]*fmat[,i];
}
}
}
parameters {
real<lower=1> gamma_adult;
real<lower=1> gamma_child;
real<lower=0, upper = lambda_max> lambda_family;
real<lower=0, upper = lambda_max> lambda_community;
}
model {
real prob_infection;
for(i in 1:N){
for(t in 1:(T-1)){
if( sum(I_mat[t]) > 0 && S[t,i]){
prob_infection = 1- ( pow( 1-lambda_family, fam_I[t,i])*pow(1-lambda_community,com_I[t,i]));
z_si[t,i] ~ bernoulli(prob_infection);
}
if(I[t,i]){
if(adult[i]) z_ir[t,i] ~ bernoulli(1/gamma_adult);
if(1-adult[i]) z_ir[t,i] ~ bernoulli(1/gamma_child);
}
}
}
}
'
#### fit model ############
mod_v2 = stan_model(model_name = 'sir_v2',model_code = sircode_v2)## In file included from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,
## from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file1d982857109b.cpp:8:
## C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"stan_data = list(
N = n,
T = t+1,
S = as.numeric(wdata[,-1] == 'S') %>% matrix(t+1,N),
I = as.numeric(wdata[,-1] == 'I') %>% matrix(t+1,N),
R = as.numeric(wdata[,-1] == 'R') %>% matrix(t+1,N),
fmat = as.matrix(fmat),
adult = adult,
lambda_max = 1,
dt = rep(1,t+1)
)
fit = sampling(mod_v2,stan_data, iter = 1000, chains = 1)##
## SAMPLING FOR MODEL 'sir_v2' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%] (Sampling)
## Elapsed Time: 73.751 seconds (Warm-up)
## 58.441 seconds (Sampling)
## 132.192 seconds (Total)stan_dens(fit)
print(fit)## Inference for Stan model: sir_v2.
## 1 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=500.
##
## mean se_mean sd 2.5% 25% 50% 75%
## gamma_adult 3.02 0.01 0.21 2.62 2.89 3.00 3.14
## gamma_child 10.38 0.04 0.88 8.63 9.80 10.39 10.92
## lambda_family 0.11 0.00 0.01 0.09 0.10 0.11 0.11
## lambda_community 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## lp__ -1480.46 0.10 1.50 -1484.07 -1481.21 -1480.09 -1479.33
## 97.5% n_eff Rhat
## gamma_adult 3.52 500 1
## gamma_child 12.18 500 1
## lambda_family 0.13 500 1
## lambda_community 0.00 500 1
## lp__ -1478.65 207 1
##
## Samples were drawn using NUTS(diag_e) at Tue May 09 12:48:16 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Again, we seem to recover the original parameters of the model. Here, I didn’t calculate \(R_0\), but it is possible with slight modification of the simple SIR model. In defining \(R_0\), note that we would expect different numbers of secondary infections depending on whether the index case is a child or an adult.
SIR Model with variable time points:
In the current work, samples are not taken uniformly over time, but instead are sampled at 10 unevenly spaced time points. Note that \(\lambda\) is a probability per-unit of time. We can replace \(\lambda\) by \(\lambda \cdot \Delta t\) everywhere in the model, where \(\Delta t\) is the number of days between observations.
This is straightforward for the transmission acquisition component of the model. For the transmission duration, we have to replace \(1/\gamma\) by \(1-(1-(1/\gamma))^{\Delta t}\).
Another wrinkle introduced is censoring: with multiple days between observations, we do not know precisely when people start and stop shedding. For now I’ll ignore this, and assess what bias might be introduced through this simulation. As we will see, in the current simulation, it appears that the bias is small relative to statistical uncertainty in the parameters.
sircode_v3 = '
data {
int<lower=0> N;
int<lower=0> T;
int<lower = 0,upper=1> S[T,N];
int<lower = 0,upper=1> I[T,N];
int<lower = 0,upper=1> R[T,N];
matrix[N,N] fmat;
int<lower=0,upper=1> adult[N];
real<lower=0> lambda_max;
real<lower=0> dt[T];
}
transformed data {
int<lower=0, upper = 1> z_si[T-1,N];
int<lower=0, upper = 1> z_ir[T-1,N];
matrix[T,N] I_mat;
matrix[T,N] fam_I;
matrix[T,N] com_I;
I_mat = to_matrix(I);
fam_I = I_mat*fmat;
for(t in 1:(T-1)){
for(i in 1:N){
z_si[t,i] = S[t,i] - S[t+1,i];
z_ir[t,i] = R[t+1,i] - R[t,i];
com_I[t,i] = sum(I_mat[t]) - I_mat[t]*fmat[,i];
}
}
}
parameters {
real<lower=1> gamma_adult;
real<lower=1> gamma_child;
real<lower=0, upper = lambda_max> lambda_family;
real<lower=0, upper = lambda_max> lambda_community;
}
model {
real prob_infection;
for(i in 1:N){
for(t in 1:(T-1)){
if( sum(I_mat[t]) > 0 && S[t,i]){
prob_infection = 1- ( pow( 1-lambda_family*dt[t], fam_I[t,i])*pow(1-lambda_community*dt[t],com_I[t,i]));
z_si[t,i] ~ bernoulli(prob_infection);
}
if(I[t,i]){
if(adult[i]) z_ir[t,i] ~ bernoulli( 1-(1-(1/gamma_adult))^dt[t] );
if(1-adult[i]) z_ir[t,i] ~ bernoulli( 1- (1-(1/gamma_child))^dt[t] );
}
}
}
}
'
#### fit model ############
mod_v3 = stan_model(model_name = 'sir_v3',model_code = sircode_v3)## In file included from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,
## from C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/arendv/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file2da0f396cbb.cpp:8:
## C:/Users/arendv/Documents/R/win-library/3.3/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"time.pts = c(1,3, 6, 10)
stan_data = list(
N = n,
T = length(time.pts),
S = as.numeric(wdata[time.pts,-1] == 'S') %>% matrix(length(time.pts),N),
I = as.numeric(wdata[time.pts,-1] == 'I') %>% matrix(length(time.pts),N),
R = as.numeric(wdata[time.pts,-1] == 'R') %>% matrix(length(time.pts),N),
fmat = as.matrix(fmat),
adult = adult,
lambda_max = 1,
dt = c(diff(time.pts),0)
)
fit = sampling(mod_v3,stan_data, iter = 1000, chains = 1, seed =3)##
## SAMPLING FOR MODEL 'sir_v3' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1, Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1, Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1, Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1, Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1, Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1, Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1, Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1, Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1, Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1, Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1, Iteration: 1000 / 1000 [100%] (Sampling)
## Elapsed Time: 1.628 seconds (Warm-up)
## 1.14 seconds (Sampling)
## 2.768 seconds (Total)
##
## [1] "The following numerical problems occurred the indicated number of times on chain 1"
## count
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1, but must be in the interval 20
## Exception thrown at line 43: bernoulli_log: Probability parameter is -0.122553, but must be in the i 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -1.22294e+006, but must be in t 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -157778, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -2.75714e+011, but must be in t 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -252399, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -491.457, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -6.88957e+010, but must be in t 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -72.3147, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -74616.2, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -870.203, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -9.85909, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is -95535.9, but must be in the in 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00001, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00004, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00014, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00024, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00044, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00074, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.00101, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.72974, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.85247e+012, but must be in th 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.99783, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 1.99832, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 2.92791e+007, but must be in th 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 2239.5, but must be in the inte 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 2371.46, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 293552, but must be in the inte 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 35969.3, but must be in the int 1
## Exception thrown at line 43: bernoulli_log: Probability parameter is 7.38859e+006, but must be in th 1
## [1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
## [1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"
## [1] "If the number in the 'count' column is small, there is no need to ask about this message on stan-users."stan_dens(fit)
print(fit)## Inference for Stan model: sir_v3.
## 1 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=500.
##
## mean se_mean sd 2.5% 25% 50% 75%
## gamma_adult 2.44 0.02 0.53 1.65 2.12 2.39 2.69
## gamma_child 10.46 0.10 2.16 7.01 8.91 10.23 11.65
## lambda_family 0.10 0.00 0.01 0.08 0.09 0.10 0.10
## lambda_community 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## lp__ -311.47 0.11 1.57 -315.17 -312.25 -311.08 -310.32
## 97.5% n_eff Rhat
## gamma_adult 3.62 500 1.00
## gamma_child 15.45 500 1.00
## lambda_family 0.12 500 1.00
## lambda_community 0.00 500 1.01
## lp__ -309.64 220 1.02
##
## Samples were drawn using NUTS(diag_e) at Tue May 09 13:05:34 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Again, the posterior distribution provides good support for the true parameters.
It is interesting that with only 4 points of observation, the posterior variability is larger, but of a comparable magnitude, to that based on the estimates with 200 points of observation.