Demonstration of package

An example of running the model with ST level character data

## Load the library
library(kilde)

## Loading required package: R2OpenBUGS

## Get some data
df <- sample_data()
## This dataset contains 4 different countries, we'll pick Canada:
df <- df[df$country == "Canada",]
## We need to format the data to the form accepted by the model function.
ob <- dataformatting_ST(DATA = df, UM = 2)
## Initialize the mcmc data objects. Note: 100 iterations is too few.
result <- initialize_mcmc_ST(ob$ns, ob$STu, MCMC = 100, ob$Nisolates)
## Run the MCMC for the ST model in R
mcmc_ob <- runmcmc_ST(result = result,
           ob = ob,
           h = 0,
           FULL = 0)
## plot the history of the MCMC
plot_history(mcmc_ob, 50)

## plot the modelfit
plot_modelfit(mcmc_ob, 50)

## A summary of the model fit
summary_kilde(mcmc_ob, 50)

## $ES
## [1] 0.0008236704
## 
## $EShum
## [1] 0.02115501

## The sample attribution
plot_sample_attribution(mcmc_ob, 50)

## a plot of the population attribution
plot_population_attribution(mcmc_ob, 50)

##
##
## We can run the same ST model in BUGS
##
## Then the same model in BUGS
## Initialize the MCMC objects
initial_result <- initialize_bugs_ST(ob)
## Run the model, this time with 1000 iterations, because bugs will
## throw an error if we try to run too few.
result_bugs <- run_bugs(result = initial_result,
                        ob = ob,
                        MCMC = 1000,
                        n.burnin = 100,
                        FULL = 0,
                        model = "SA_ST_model.jag",
                        n.chains = 1)
plot_history(result_bugs, 100)

plot_modelfit(result_bugs, 100)

summary_kilde(result_bugs, 100)

## $ES
## [1] 0.0001580709
## 
## $EShum
## [1] 0.02108002

plot_sample_attribution(result_bugs, 100)

plot_population_attribution(result_bugs, 100)

### The model estimated by using the Allele level results

library(kilde)
## Read in a format data
df <- sample_data()
df <- df[df$country == "Canada",]
ob <- dataformatting(DATA = df,
                     UM = 2)
######################################
## Initialize and then run mcmc in R
######################################
result <- initialize_mcmc(ns = ob$inits$ns,
                          nat = ob$inits$nat,
                          MCMC = 100,
                          Nisolates = ob$inits$Nisolates)
mcmc_ob <- runmcmc(result, ob, MCMC = 100, h = 0, FULL = 0)
##  Plot the results of this model.
plot_history(mcmc_ob, 50)

plot_modelfit(mcmc_ob, 50)

summary_kilde(mcmc_ob, 50)

## $ES
## [1] 0.005014784
## 
## $EShum
## [1] 0.0899753

plot_sample_attribution(mcmc_ob, 50)

plot_population_attribution(mcmc_ob, 50)

##
################################################
## Initialize and then run the model in OpenBugs
################################################
##
## Below, BUGS model cannot handle a large number of MCMC iterations
## for all parameters.  Therefore, it is advisable to try with smaller
## number of iterations to start, perhaps 1000.
##
initial_result <- initialize_bugs(ob)
result_bugs <- run_bugs(result = initial_result,
                        ob = ob,
                        MCMC = 1000,
                        n.burnin = 100,
                        FULL = 0)
plot_history(result_bugs, 100)

plot_modelfit(result_bugs, 100)

summary_kilde(result_bugs, 100)

## $ES
## [1] 0.002545353
## 
## $EShum
## [1] 0.09235274

plot_sample_attribution(result_bugs, 100)

plot_population_attribution(result_bugs, 100)