# NOT RUN {
#### Format data using Y, X, and Prev functions #############################
## Input data must be in long format
y <- Y( # Cases
data = sim_SA$cases,
y = "Human",
type = "Type",
time = "Time",
location = "Location"
)
x <- X( # Sources
data = sim_SA$sources,
x = "Count",
type = "Type",
time = "Time",
source = "Source"
)
k <- Prev( # Prevalences
data = sim_SA$prev,
prev = "Value",
time = "Time",
source = "Source"
)
#### Create Dirichlet(1) priors #############################################
## Create alpha prior data frame
prior_alpha_long <- expand.grid(
Source = unique(sim_SA$sources$Source),
Time = unique(sim_SA$sources$Time),
Location = unique(sim_SA$cases$Location),
Alpha = 1
)
# Use the Alpha() constructor to specify alpha prior
prior_alpha <- Alpha(
data = prior_alpha_long,
alpha = 'Alpha',
source = 'Source',
time = 'Time',
location = 'Location'
)
## Create r prior data frame
prior_r_long <- expand.grid(
Type = unique(sim_SA$sources$Type),
Source = unique(sim_SA$sources$Source),
Time = unique(sim_SA$sources$Time),
Value = 0.1
)
# Use X() constructor to specify r prior
prior_r <- X(
data = prior_r_long,
x = 'Value',
type = 'Type',
time = 'Time',
source = 'Source'
)
## Pack all priors into a list
priors <- list(
a_theta = 0.01,
b_theta = 0.00001,
a_alpha = prior_alpha,
a_r = prior_r
)
## If all prior values are the same, they can be specified in shorthand
## Equivalent result to the longform priors specified above
priors <- list(
a_theta = 0.01,
b_theta = 0.00001,
a_alpha = 1,
a_r = 0.1
)
#### Set initial values (optional) ##########################################
types <- unique(sim_SA$cases$Type)
q_long <- data.frame(q=rep(15, length(types)), Type=types)
init_q <- Q(q_long, q = 'q', type = 'Type')
inits <- list(q = init_q) # Pack starting values into a list
#### Construct model ########################################################
my_model <- HaldDP(y = y, x = x, k = k, priors = priors, inits = inits, a_q = 0.1)
#### Set mcmc parameters ####################################################
my_model$mcmc_params(n_iter = 2, burn_in = 2, thin = 1)
#### Update model ###########################################################
my_model$update()
## Add an additional 10 iterations
my_model$update(n_iter = 2, append = TRUE)
#### Extract posterior ######################################################
## returns the posterior for the r, alpha, q, c,
## lambda_i, xi and xi_prop parameters,
## for all times, locations, sources and types
## the posterior is returned as a list or arrays
# }
# NOT RUN {
str(my_model$extract())
# }
# NOT RUN {
## returns the posterior for the r and alpha parameters,
## for time 1, location B, sources Source3, and Source4,
## types 5, 25, and 50, and iterations 200:300
## the posterior is returned as a list of dataframes
# }
# NOT RUN {
str(my_model$extract(params = c("r", "alpha"),
times = "1", location = "B",
sources = c("Source3", "Source4"),
types = c("5", "25", "50"),
iters = 5:15,
flatten = TRUE))
# }
# NOT RUN {
#### Calculate medians and credible intervals ###############################
# }
# NOT RUN {
my_model$summary(alpha = 0.05, CI_type = "chen-shao")
# }
# NOT RUN {
## subsetting is done in the same way as extract()
# }
# NOT RUN {
my_model$summary(alpha = 0.05, CI_type = "chen-shao",
params = c("r", "alpha"),
times = "1", location = "B",
sources = c("Source3", "Source4"),
types = c("5", "25", "50"),
iters = 5:15,
flatten = TRUE)
# }
# NOT RUN {
#### Plot heatmap and dendrogram of the type effect grouping ################
my_model$plot_heatmap()
#### Extract data, initial values, prior values, acceptance
## rates for the mcmc algorithm and mcmc parameters
my_model$get_data()
my_model$get_inits()
my_model$get_priors()
my_model$get_acceptance()
my_model$get_mcmc_params()
# }
Run the code above in your browser using DataLab