In grantbrown/ABSEIR: Posterior Approximation to Fully Bayesian Spatial SEIR Models Using Multithreaded C++
Note: This tutorial assumes that you have successfully installed ABSEIR, and are
at least passingly familiar with compartmental models. Some introductory information is available
in this vignette.
Step 1: Setup
The data set we will consider here records daily case counts of a 1995 outbreak of Ebola in the city of Kikwit, in the Democratic Republic of the Congo. A classic analysis of this data is
available in the work of Lekone and Finkenstädt, 2006. A copy is provided along with the ABSEIR library:
# load the ABSEIR library
library(ABSEIR)
# read in the data set
data(Kikwit1995)
Step 2: Exploration
We can now take a look at the raw case counts:
# make barplot of new case counts
barplot(t(Kikwit1995$Count), main = "New Cases", xlab = "Day", ylab = "Cases")
# add an X axis to give a reference epidemic time
axis(side = 1, at = seq(0, nrow(Kikwit1995), 50))
Step 3: Build some plausible models
We're now ready to start specifying models. Unfortunately, there are quite a
few moving parts to an epidemic model, so we'll address them one at a time.
Data Model
A data model describes the way in which the observed epidemic outcome relates
to the underlying model. This single location analysis is defined over time
points ${t_i : i = 1,...,T}$, and divides individuals in the population
into four categories:
- Susceptible: individuals capable of contracting an illness
- Exposed: individuals who have contracted the pathogen, but are not yet infectious
- Infectious: individuals who are capable of spreading the pathogen
- Removed: individuals who are no longer infectious, either through mortality or recovery with immunity.
Mathematically, we keep track of these counts in the $T\times 1$ column vectors:
$\bf{S}$, $\bf{E}$, $\bf{I}$, and $\bf{R}$. In spatial analyses, these become $T \times n$ matrices,
where $n$ is equal to the number of spatial locations.
Our model needs to keep track of transitions between these categories, which we denote
with asterisks. This gives rise to the following transition equations:
$$
S_{i+1} = S_{i} - E^{i} \
E{i+1} = E_{i} + E^{i} - I^{i} \
I{i+1} = I_{i} + I^{i} - R^{i} \
R{i+1} = R_{i} + R^_{i}
$$
The reason for all this exposition and discussion of notation is that
it explains the way we've allowed users to configure data models. In particular,
we implement models for "I_star", or new infectious counts, "I", or current infectious
size, and "R_star", or the newly removed individuals. In the case of the
Ebola data, we believe that it most closely resembles new infectious counts, and
therefore choose the "I_star" compartment type.
In addition, we're going to choose the simplest type of data model, denoted
"identity". The identity data model assumes that we observe, without error,
the values of the chosen epidemic compartment. We also allow users to model
this quantity under overdispersion (i.e., with noise).
Finally, users can specify whether a compartment is reported on a cumulative scale.
In this case, counts are not cumulative. Additional details are available
in the help page of the DataModel function.
# Create a model to relate observe data to epidemic process
data_model = DataModel(Kikwit1995$Count,
type = "identity",
compartment="I_star",
cumulative=FALSE)
Exposure Model
The next model component we need to specify is the "exposure model", which descscribes
factors relating to changes in epidemic intensity. This model component specifically
relates to the probability of exposure, or moving from the susceptible to exposed
category. Specification of the exposure model requires careful thought, and can be a
little tricky. In particular, we need to specify the "design matrix" of covariates.
In the single location case, this is considerably simplified. The single location
design matrix is a $T \times p$ matrix, where $p$ is the number of covariates and
$T$ is the number of time points. Here we will consider two exposure models. The
first assumes a simple intercept model, and the second introduces a covariate
which captures the effect of a public health intervention which was launched
around May 9th.
In addition to the need to determine the temporal (and potentially spatial)
structure of the epidemic intensity, we need to specify prior distributions
for the associated linear parameters. These are modeled as normal random
variables with mean vector "betaPriorMean" and precision vector "betaPriorPrecision".
# Create a model to describe relationship of any covariates to epidemic intensity
exposure_model_1 = ExposureModel(matrix(1, nrow = nrow(Kikwit1995)),
nTpt = nrow(Kikwit1995),
nLoc = 1,
betaPriorPrecision = 0.5,
betaPriorMean = 0)
# Create a (more realistic) Model to describe relationship of any covariates to epidemic intensity
# In this analysis we know on which date interventions began, so we may include a linear term
# starting on that date.
intervention_term = cumsum(Kikwit1995$Date > as.Date("05-09-1995", "%m-%d-%Y"))
exposure_model_2 = ExposureModel(cbind(1,intervention_term),
nTpt = nrow(Kikwit1995),
nLoc = 1,
betaPriorPrecision = 0.5,
betaPriorMean = 0)
Reinfection Model
Reinfection refers to the process of losing immunity, and returning to the susceptible population.
This isn't generally considered a concern for Ebola, so we'll just employ a SEIR model, rather
than an SEIRS model.
# There's no reinfection in this case, so we just use a "SEIR" model.
reinfection_model = ReinfectionModel("SEIR")
Distance Model
A distance model is the mechanism by which spatial dependence is introduced.
In this analysis, we have only one spatial location, so the "distance model"
is just a placeholder. We'll initialize it with an empty matrix.
distance_model = DistanceModel(list(matrix(0)))
Initial Values
The transition equations described above are helpful for understanding how
the epidemic models fit by ABSEIR progress, but they don't describe
how the models are initialized. Models require the specification of an initial
state: the number of susceptible, exposed, infectious, and removed individuals
present at the beginning of the anlysis.
# Set initial population sizes
initial_value_container = InitialValueContainer(S0=5.36e6,
E0=2,
I0=2,
R0=0)
Other Transitions
Thus far, we have said little of the transition processes
between the E, I and R compartments. These processes are determined
by specifying a "transition model". Transition models are implemented
in several forms. The simplest form is provided by the
"ExponentialTransitionPriors" function. The "exponential"
in "ExponentialTransitionPriors" refers to the assumption that the time
individuals spend in the latent and infectious states is exponentially
distributed. This is a common assumption in compartmental models, and
provides computational benefits. Neverthe less, it may be unrealistic
for some pathogens. We'll stick to the exponential case for the purposes
of this tutorial, but ABSEIR also allows fully parameterized Weibull
membership times and arbitrarily flexible fixed membership distributions
using path specific compartmental modeling techniques.
Exponential transition models are specified by including an E to I
and I to R transition probability, and an associated effective sample
size for each.
transition_priors = ExponentialTransitionPriors(p_ei = 1-exp(-1/5),
p_ir= 1-exp(-1/7),
p_ei_ess = 10,
p_ir_ess = 10)
Sampling Control
Finally, we need to specify and tune the algorithm ABSEIR should use to perform
inference. The basic ABC rejection algorithm simulates epidemics from the prior
distribution, and selects parameter values producing epidemics most closely
resembling the data provided. We prefer to use the sequential Monte Carlo
algorithm of Beaumont 2009, though both algorithms are implemented by ABSEIR.
Other algorithms may be introduced if a need is identified.
The sequential Monte Carlo approach requires the specification of several tuning
parameters:
- n_cores: Number of CPU cores to use for simulation.
- batch_size: Number of epidemics to simulate in parallel
before reevaluating convergence.
- epochs: Maximum number of complete updates to perform before declaring convergence
- max_batches: Maximum number of incomplete sample batches to perform before declaring convergence
- shrinkage: Multiplicative decrease in maximum distance between simulated and observed epidemics between iterations
- multivariate_perturbation: Logical indicator determining whether new proposed parameter values are made
from a univariate or multivariate normal distribution.
# Set algorithm configuration
sampling_control = SamplingControl(seed = 123123,
n_cores = 8,
algorithm="Beaumont2009",
list(batch_size = 2000,
epochs = 1e6,
max_batches = 100,
shrinkage = 0.99,
multivariate_perturbation=FALSE,
keep_compartments = TRUE
)
)
Run the Models
We can now combine these components and fit the models. In an interactive session,
users can request verbose output to get an idea of the convergence status of the model.
# Underspecified intensity
runtime1 = system.time(result1 <- SpatialSEIRModel(data_model,
exposure_model_1,
reinfection_model,
distance_model,
transition_priors,
initial_value_container,
sampling_control,
samples = 100,
verbose = 0))
# user system elapsed
# 2258.391 17.449 294.584
# Reasonable intensity
runtime2 = system.time(result2 <- SpatialSEIRModel(data_model,
exposure_model_2,
reinfection_model,
distance_model,
transition_priors,
initial_value_container,
sampling_control,
samples = 100,
verbose = 0))
Let's take a look at the required runtimes:
timeMatrix = rbind(runtime1,runtime2)
rownames(timeMatrix) = paste("model", 1:2)
print(timeMatrix[,1:3])
We can also get a summary of parameter estimates:
summary(result1)
summary(result2)
Model Comparison
Models can be formally compared using Bayes Factors,
and by producing posterior predictive plots. We'll take a look at both techniques.
compareModels(list(result1, result2), n_samples = 1000,
batch_size = 2000)
We observe strong evidence for model 2 over model 1 in this case, which is unsurprising given the presence of a public health intervention.
We can also consider model fit by looking at posterior predictive distribution of epidemic values.
# Simulate new epidemics based on the accepted parameters for model 1
simulations1 <- epidemic.simulations(result1, replicates = 50)
simulations2 <- epidemic.simulations(result2, replicates = 50)
# Add reproductive number estimates to simulation results
system.time(simulations1 <- ComputeR0(simulations1))
system.time(simulations2 <- ComputeR0(simulations2))
plotPosteriorPredictive = function(simulations, main)
{
allSimulatedI_star = sapply(simulations$simulationResults, function(x){x$I_star})
lowerQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.025))
posteriorMean = apply(allSimulatedI_star, 1, mean)
upperQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.975))
plot(Kikwit1995$Count, ylim = c(0, max(Kikwit1995$Count)*2),
xlab = "Epidemic Day", ylab = "New Cases", main = main)
lines(upperQuantile, lty = 2, col = "blue")
lines(lowerQuantile, lty = 2, col = "blue")
lines(posteriorMean, lty = 1, col = "blue")
legend(x = 100, y = 12, legend = c("Mean", "95% CI", "Observed"), lty = c(1,2,0),
pch = c(NA,NA,1), col = c("blue", "blue", "black"), cex = 1)
}
plotPosteriorPredictive(simulations1, "Model 1: Posterior Predictive Distribution")
plotPosteriorPredictive(result1, "Model 1: Posterior Distribution")
plotPosteriorPredictive(simulations2, "Model 2: Posterior Predictive Distribution")
plotPosteriorPredictive(result2, "Model 2: Posterior Distribution")
Reproductive Numbers
Reproductive numbers are an important tool to understand the behavior of epidemic processes, and
capture the number of secondary infectious caused by a single infectious individual. Beyond that
simple description, reproductive numbers differ in many details. We provide estimates of the time
varying reproductive number, the effective reproductive number, and the empirically adjusted
reproductive number. A more complete discussion of these quantities is available here.
simulations1.R0 <- ComputeR0(simulations1, cores = 8)
simulations2.R0 <- ComputeR0(simulations2, cores = 8)
plotR0 = function(simulations, main)
{
allSimulatedEA_R0 = sapply(simulations$simulationResults, function(x){x$R_EA})
plot(apply(allSimulatedEA_R0, 1, mean), type = "l", ylim = c(0, 3), lwd =2,
ylab = "Reproductive Number", main = main)
lines(apply(allSimulatedEA_R0, 1, mean), lwd = 2, lty = 2, col = "blue")
lines(apply(allSimulatedEA_R0, 1, quantile, probs = c(0.1)), lwd = 2, lty = 2, col = "blue")
lines(apply(allSimulatedEA_R0, 1, quantile, probs = c(0.9)), lwd = 2, lty = 2, col = "blue")
}
plotR0(simulations1.R0, "Model 1: EA-R(t)")
plotR0(simulations2.R0, "Model 2: EA-R(t)")
Here, we also observe that model 2 provides much more reasonable reproductive number curves. Of course, numerous other modeling
choices are possible. In particular, we did not explore different potential transition models. This can be easily accomplished,
in much the same way that we examined two potential exposure models in this example.
introduction to working with ABSEIR.
We can also calculate the basic reproductive number in each of these models, though this requires us to be clear about
how that should be defined, and what condition should be used to calculate it.
p1 <- as.data.frame(result1$param.samples)
p2 <- as.data.frame(result2$param.samples)
hist(exp(p1$Beta_SE_1)/p1$gamma_IR, main= "Model 1 R0 - Posterior Samples")
hist(exp(p2$Beta_SE_1)/p2$gamma_IR, main= "Model 2 R0 Posterior Samples")
Questions?
Questions and suggestions can be directed to the issue page of the ABSEIR repository.
grantbrown/ABSEIR documentation built on Oct. 14, 2021, 2:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Note: This tutorial assumes that you have successfully installed ABSEIR, and are at least passingly familiar with compartmental models. Some introductory information is available in this vignette.
Step 1: Setup
The data set we will consider here records daily case counts of a 1995 outbreak of Ebola in the city of Kikwit, in the Democratic Republic of the Congo. A classic analysis of this data is available in the work of Lekone and Finkenstädt, 2006. A copy is provided along with the ABSEIR library:
# load the ABSEIR library library(ABSEIR) # read in the data set data(Kikwit1995)
Step 2: Exploration
We can now take a look at the raw case counts:
# make barplot of new case counts barplot(t(Kikwit1995$Count), main = "New Cases", xlab = "Day", ylab = "Cases") # add an X axis to give a reference epidemic time axis(side = 1, at = seq(0, nrow(Kikwit1995), 50))
Step 3: Build some plausible models
We're now ready to start specifying models. Unfortunately, there are quite a few moving parts to an epidemic model, so we'll address them one at a time.
Data Model
A data model describes the way in which the observed epidemic outcome relates to the underlying model. This single location analysis is defined over time points ${t_i : i = 1,...,T}$, and divides individuals in the population into four categories:
- Susceptible: individuals capable of contracting an illness
- Exposed: individuals who have contracted the pathogen, but are not yet infectious
- Infectious: individuals who are capable of spreading the pathogen
- Removed: individuals who are no longer infectious, either through mortality or recovery with immunity.
Mathematically, we keep track of these counts in the $T\times 1$ column vectors: $\bf{S}$, $\bf{E}$, $\bf{I}$, and $\bf{R}$. In spatial analyses, these become $T \times n$ matrices, where $n$ is equal to the number of spatial locations.
Our model needs to keep track of transitions between these categories, which we denote with asterisks. This gives rise to the following transition equations:
$$ S_{i+1} = S_{i} - E^{i} \ E{i+1} = E_{i} + E^{i} - I^{i} \ I{i+1} = I_{i} + I^{i} - R^{i} \ R{i+1} = R_{i} + R^_{i} $$
The reason for all this exposition and discussion of notation is that it explains the way we've allowed users to configure data models. In particular, we implement models for "I_star", or new infectious counts, "I", or current infectious size, and "R_star", or the newly removed individuals. In the case of the Ebola data, we believe that it most closely resembles new infectious counts, and therefore choose the "I_star" compartment type.
In addition, we're going to choose the simplest type of data model, denoted "identity". The identity data model assumes that we observe, without error, the values of the chosen epidemic compartment. We also allow users to model this quantity under overdispersion (i.e., with noise).
Finally, users can specify whether a compartment is reported on a cumulative scale. In this case, counts are not cumulative. Additional details are available in the help page of the DataModel function.
# Create a model to relate observe data to epidemic process data_model = DataModel(Kikwit1995$Count, type = "identity", compartment="I_star", cumulative=FALSE)
Exposure Model
The next model component we need to specify is the "exposure model", which descscribes factors relating to changes in epidemic intensity. This model component specifically relates to the probability of exposure, or moving from the susceptible to exposed category. Specification of the exposure model requires careful thought, and can be a little tricky. In particular, we need to specify the "design matrix" of covariates. In the single location case, this is considerably simplified. The single location design matrix is a $T \times p$ matrix, where $p$ is the number of covariates and $T$ is the number of time points. Here we will consider two exposure models. The first assumes a simple intercept model, and the second introduces a covariate which captures the effect of a public health intervention which was launched around May 9th.
In addition to the need to determine the temporal (and potentially spatial) structure of the epidemic intensity, we need to specify prior distributions for the associated linear parameters. These are modeled as normal random variables with mean vector "betaPriorMean" and precision vector "betaPriorPrecision".
# Create a model to describe relationship of any covariates to epidemic intensity exposure_model_1 = ExposureModel(matrix(1, nrow = nrow(Kikwit1995)), nTpt = nrow(Kikwit1995), nLoc = 1, betaPriorPrecision = 0.5, betaPriorMean = 0) # Create a (more realistic) Model to describe relationship of any covariates to epidemic intensity # In this analysis we know on which date interventions began, so we may include a linear term # starting on that date. intervention_term = cumsum(Kikwit1995$Date > as.Date("05-09-1995", "%m-%d-%Y")) exposure_model_2 = ExposureModel(cbind(1,intervention_term), nTpt = nrow(Kikwit1995), nLoc = 1, betaPriorPrecision = 0.5, betaPriorMean = 0)
Reinfection Model
Reinfection refers to the process of losing immunity, and returning to the susceptible population. This isn't generally considered a concern for Ebola, so we'll just employ a SEIR model, rather than an SEIRS model.
# There's no reinfection in this case, so we just use a "SEIR" model. reinfection_model = ReinfectionModel("SEIR")
Distance Model
A distance model is the mechanism by which spatial dependence is introduced. In this analysis, we have only one spatial location, so the "distance model" is just a placeholder. We'll initialize it with an empty matrix.
distance_model = DistanceModel(list(matrix(0)))
Initial Values
The transition equations described above are helpful for understanding how the epidemic models fit by ABSEIR progress, but they don't describe how the models are initialized. Models require the specification of an initial state: the number of susceptible, exposed, infectious, and removed individuals present at the beginning of the anlysis.
# Set initial population sizes initial_value_container = InitialValueContainer(S0=5.36e6, E0=2, I0=2, R0=0)
Other Transitions
Thus far, we have said little of the transition processes between the E, I and R compartments. These processes are determined by specifying a "transition model". Transition models are implemented in several forms. The simplest form is provided by the "ExponentialTransitionPriors" function. The "exponential" in "ExponentialTransitionPriors" refers to the assumption that the time individuals spend in the latent and infectious states is exponentially distributed. This is a common assumption in compartmental models, and provides computational benefits. Neverthe less, it may be unrealistic for some pathogens. We'll stick to the exponential case for the purposes of this tutorial, but ABSEIR also allows fully parameterized Weibull membership times and arbitrarily flexible fixed membership distributions using path specific compartmental modeling techniques.
Exponential transition models are specified by including an E to I and I to R transition probability, and an associated effective sample size for each.
transition_priors = ExponentialTransitionPriors(p_ei = 1-exp(-1/5), p_ir= 1-exp(-1/7), p_ei_ess = 10, p_ir_ess = 10)
Sampling Control
Finally, we need to specify and tune the algorithm ABSEIR should use to perform inference. The basic ABC rejection algorithm simulates epidemics from the prior distribution, and selects parameter values producing epidemics most closely resembling the data provided. We prefer to use the sequential Monte Carlo algorithm of Beaumont 2009, though both algorithms are implemented by ABSEIR. Other algorithms may be introduced if a need is identified.
The sequential Monte Carlo approach requires the specification of several tuning parameters:
- n_cores: Number of CPU cores to use for simulation.
- batch_size: Number of epidemics to simulate in parallel before reevaluating convergence.
- epochs: Maximum number of complete updates to perform before declaring convergence
- max_batches: Maximum number of incomplete sample batches to perform before declaring convergence
- shrinkage: Multiplicative decrease in maximum distance between simulated and observed epidemics between iterations
- multivariate_perturbation: Logical indicator determining whether new proposed parameter values are made from a univariate or multivariate normal distribution.
# Set algorithm configuration sampling_control = SamplingControl(seed = 123123, n_cores = 8, algorithm="Beaumont2009", list(batch_size = 2000, epochs = 1e6, max_batches = 100, shrinkage = 0.99, multivariate_perturbation=FALSE, keep_compartments = TRUE ) )
Run the Models
We can now combine these components and fit the models. In an interactive session, users can request verbose output to get an idea of the convergence status of the model.
# Underspecified intensity runtime1 = system.time(result1 <- SpatialSEIRModel(data_model, exposure_model_1, reinfection_model, distance_model, transition_priors, initial_value_container, sampling_control, samples = 100, verbose = 0)) # user system elapsed # 2258.391 17.449 294.584 # Reasonable intensity runtime2 = system.time(result2 <- SpatialSEIRModel(data_model, exposure_model_2, reinfection_model, distance_model, transition_priors, initial_value_container, sampling_control, samples = 100, verbose = 0))
Let's take a look at the required runtimes:
timeMatrix = rbind(runtime1,runtime2) rownames(timeMatrix) = paste("model", 1:2) print(timeMatrix[,1:3])
We can also get a summary of parameter estimates:
summary(result1) summary(result2)
Model Comparison
Models can be formally compared using Bayes Factors, and by producing posterior predictive plots. We'll take a look at both techniques.
compareModels(list(result1, result2), n_samples = 1000, batch_size = 2000)
We observe strong evidence for model 2 over model 1 in this case, which is unsurprising given the presence of a public health intervention. We can also consider model fit by looking at posterior predictive distribution of epidemic values.
# Simulate new epidemics based on the accepted parameters for model 1 simulations1 <- epidemic.simulations(result1, replicates = 50) simulations2 <- epidemic.simulations(result2, replicates = 50) # Add reproductive number estimates to simulation results system.time(simulations1 <- ComputeR0(simulations1)) system.time(simulations2 <- ComputeR0(simulations2)) plotPosteriorPredictive = function(simulations, main) { allSimulatedI_star = sapply(simulations$simulationResults, function(x){x$I_star}) lowerQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.025)) posteriorMean = apply(allSimulatedI_star, 1, mean) upperQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.975)) plot(Kikwit1995$Count, ylim = c(0, max(Kikwit1995$Count)*2), xlab = "Epidemic Day", ylab = "New Cases", main = main) lines(upperQuantile, lty = 2, col = "blue") lines(lowerQuantile, lty = 2, col = "blue") lines(posteriorMean, lty = 1, col = "blue") legend(x = 100, y = 12, legend = c("Mean", "95% CI", "Observed"), lty = c(1,2,0), pch = c(NA,NA,1), col = c("blue", "blue", "black"), cex = 1) } plotPosteriorPredictive(simulations1, "Model 1: Posterior Predictive Distribution") plotPosteriorPredictive(result1, "Model 1: Posterior Distribution") plotPosteriorPredictive(simulations2, "Model 2: Posterior Predictive Distribution") plotPosteriorPredictive(result2, "Model 2: Posterior Distribution")
Reproductive Numbers
Reproductive numbers are an important tool to understand the behavior of epidemic processes, and capture the number of secondary infectious caused by a single infectious individual. Beyond that simple description, reproductive numbers differ in many details. We provide estimates of the time varying reproductive number, the effective reproductive number, and the empirically adjusted reproductive number. A more complete discussion of these quantities is available here.
simulations1.R0 <- ComputeR0(simulations1, cores = 8) simulations2.R0 <- ComputeR0(simulations2, cores = 8) plotR0 = function(simulations, main) { allSimulatedEA_R0 = sapply(simulations$simulationResults, function(x){x$R_EA}) plot(apply(allSimulatedEA_R0, 1, mean), type = "l", ylim = c(0, 3), lwd =2, ylab = "Reproductive Number", main = main) lines(apply(allSimulatedEA_R0, 1, mean), lwd = 2, lty = 2, col = "blue") lines(apply(allSimulatedEA_R0, 1, quantile, probs = c(0.1)), lwd = 2, lty = 2, col = "blue") lines(apply(allSimulatedEA_R0, 1, quantile, probs = c(0.9)), lwd = 2, lty = 2, col = "blue") } plotR0(simulations1.R0, "Model 1: EA-R(t)") plotR0(simulations2.R0, "Model 2: EA-R(t)")
Here, we also observe that model 2 provides much more reasonable reproductive number curves. Of course, numerous other modeling choices are possible. In particular, we did not explore different potential transition models. This can be easily accomplished, in much the same way that we examined two potential exposure models in this example. introduction to working with ABSEIR.
We can also calculate the basic reproductive number in each of these models, though this requires us to be clear about how that should be defined, and what condition should be used to calculate it.
p1 <- as.data.frame(result1$param.samples) p2 <- as.data.frame(result2$param.samples) hist(exp(p1$Beta_SE_1)/p1$gamma_IR, main= "Model 1 R0 - Posterior Samples") hist(exp(p2$Beta_SE_1)/p2$gamma_IR, main= "Model 2 R0 Posterior Samples")
Questions?
Questions and suggestions can be directed to the issue page of the ABSEIR repository.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.