Preparing CRAN submission

DESCRIPTION

@@ -1,13 +1,13 @@
 Package: epiworldR
 Type: Package
 Title: Fast Agent-Based Epi Models
-Version: 0.6-0.0
+Version: 0.6.0.0
 Authors@R: c(
-  person(given="George", family="Vega Yon", role=c("aut","cre"),
+  person(given="George", family="Vega Yon", role=c("aut"),
   email="[email protected]", comment = c(ORCID = "0000-0002-3171-0844")),
   person(given="Derek", family="Meyer", role=c("aut"),
   email="[email protected]", comment = c(ORCID = "0009-0005-1350-6988")),
-  person(given="Andrew", family="Pulsipher", role=c("aut"),
+  person(given="Andrew", family="Pulsipher", role=c("aut", "cre"),
   email="[email protected]", comment = c(ORCID = "0000-0002-0773-3210")),
   person(given="Susan", family="Holmes", role = "rev", comment =
          c(what = "JOSS reviewer", ORCID="0000-0002-2208-8168")),

NEWS.md

@@ -1,4 +1,4 @@
-# epiworldR 0.4-3 (development version)
+# epiworldR 0.6.0.0
 
 ## New features
 
@@ -14,6 +14,12 @@
 
 * The function `today()` returns the current day (step) of the
   simulation.
+
+## Misc
+
+* We changed the versioning system. To allow the R package to increase
+  version number while preserving epiworld (C++) versioning, we added a fourth
+  number that indicates R-only patches (similar to RcppArmadillo).
 
 
 # epiworldR 0.3-2

README.md

@@ -123,7 +123,7 @@ sir
 #> Susceptible-Infected-Recovered (SIR)
 #> It features 100000 agents, 1 virus(es), and 0 tool(s).
 #> The model has 3 states.
-#> The final distribution is: 822 Susceptible, 415 Infected, and 98763 Recovered.
+#> The final distribution is: 1209 Susceptible, 499 Infected, and 98292 Recovered.
 ```
 
 Visualizing the outputs
@@ -133,37 +133,37 @@ summary(sir)
 #> ________________________________________________________________________________
 #> ________________________________________________________________________________
 #> SIMULATION STUDY
-#>
+#> 
 #> Name of the model   : Susceptible-Infected-Recovered (SIR)
 #> Population size     : 100000
 #> Agents' data        : (none)
 #> Number of entities  : 0
 #> Days (duration)     : 50 (of 50)
 #> Number of viruses   : 1
-#> Last run elapsed t  : 62.00ms
-#> Last run speed      : 79.47 million agents x day / second
+#> Last run elapsed t  : 66.00ms
+#> Last run speed      : 74.96 million agents x day / second
 #> Rewiring            : off
-#>
+#> 
 #> Global events:
 #>  (none)
-#>
+#> 
 #> Virus(es):
 #>  - COVID-19
-#>
+#> 
 #> Tool(s):
 #>  (none)
-#>
+#> 
 #> Model parameters:
 #>  - Recovery rate     : 0.3000
 #>  - Transmission rate : 0.7000
-#>
+#> 
 #> Distribution of the population at time 50:
-#>   - (0) Susceptible :  99000 -> 822
-#>   - (1) Infected    :   1000 -> 415
-#>   - (2) Recovered   :      0 -> 98763
-#>
+#>   - (0) Susceptible :  99000 -> 1209
+#>   - (1) Infected    :   1000 -> 499
+#>   - (2) Recovered   :      0 -> 98292
+#> 
 #> Transition Probabilities:
-#>  - Susceptible  0.91  0.09  0.00
+#>  - Susceptible  0.92  0.08  0.00
 #>  - Infected     0.00  0.70  0.30
 #>  - Recovered    0.00  0.00  1.00
 ```
@@ -221,39 +221,39 @@ summary(model_seirconn)
 #> ________________________________________________________________________________
 #> ________________________________________________________________________________
 #> SIMULATION STUDY
-#>
+#> 
 #> Name of the model   : Susceptible-Exposed-Infected-Removed (SEIR) (connected)
 #> Population size     : 10000
 #> Agents' data        : (none)
 #> Number of entities  : 0
 #> Days (duration)     : 100 (of 100)
 #> Number of viruses   : 2
 #> Last run elapsed t  : 15.00ms
-#> Last run speed      : 64.01 million agents x day / second
+#> Last run speed      : 65.69 million agents x day / second
 #> Rewiring            : off
-#>
+#> 
 #> Global events:
 #>  - Update infected individuals (runs daily)
-#>
+#> 
 #> Virus(es):
 #>  - COVID-19
 #>  - COVID-19
-#>
+#> 
 #> Tool(s):
 #>  (none)
-#>
+#> 
 #> Model parameters:
 #>  - Avg. Incubation days : 7.0000
 #>  - Contact rate         : 10.0000
 #>  - Prob. Recovery       : 0.1429
 #>  - Prob. Transmission   : 0.1000
-#>
+#> 
 #> Distribution of the population at time 100:
 #>   - (0) Susceptible :  9800 -> 11
 #>   - (1) Exposed     :   200 -> 0
 #>   - (2) Infected    :     0 -> 3
 #>   - (3) Recovered   :     0 -> 9986
-#>
+#> 
 #> Transition Probabilities:
 #>  - Susceptible  0.94  0.06  0.00  0.00
 #>  - Exposed      0.00  0.85  0.15  0.00
@@ -367,25 +367,25 @@ rn <- get_reproductive_number(model_logit)
   X[, "Female"],
   (1:n %in% rn$source)
 ) |> prop.table())[, 2]
-#>       0       1
-#> 0.12984 0.14201
+#>       0       1 
+#> 0.13466 0.14878
 ```
 
 ``` r
 
 # Looking into the agents
 get_agents(model_logit)
 #> Agents from the model "Susceptible-Infected-Removed (SIR) (logit)":
-#> Agent: 0, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 1, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 2, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 3, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 4, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 5, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 6, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 7, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 0, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 1, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 2, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 3, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 4, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 5, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 6, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 7, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
 #> Agent: 8, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
-#> Agent: 9, state: Recovered (2), Has virus: no, NTools: 0i NNeigh: 8
+#> Agent: 9, state: Susceptible (0), Has virus: no, NTools: 0i NNeigh: 8
 #> ... 99990 more agents ...
 ```
 
@@ -477,20 +477,20 @@ head(ans$total_hist)
 #> 1       1    0        1 Susceptible    990
 #> 2       1    0        1    Infected     10
 #> 3       1    0        1   Recovered      0
-#> 4       1    1        1 Susceptible    971
-#> 5       1    1        1    Infected     29
+#> 4       1    1        1 Susceptible    974
+#> 5       1    1        1    Infected     26
 #> 6       1    1        1   Recovered      0
 ```
 
 ``` r
 head(ans$reproductive)
 #>   sim_num virus_id    virus source source_exposure_date rt
-#> 1       1        0 COVID-19    943                    9  0
-#> 2       1        0 COVID-19     40                    9  0
-#> 3       1        0 COVID-19      6                    9  0
-#> 4       1        0 COVID-19    973                    8  0
-#> 5       1        0 COVID-19    495                    9  0
-#> 6       1        0 COVID-19    480                    8  0
+#> 1       1        0 COVID-19     77                   10  0
+#> 2       1        0 COVID-19    970                    8  0
+#> 3       1        0 COVID-19    895                    8  0
+#> 4       1        0 COVID-19    866                    8  0
+#> 5       1        0 COVID-19    811                    8  0
+#> 6       1        0 COVID-19    752                    8  0
 ```
 
 ``` r
@@ -515,16 +515,17 @@ If you use `epiworldR` in your research, please cite it as follows:
 ``` r
 citation("epiworldR")
 #> To cite epiworldR in publications use:
-#>
+#> 
 #>   Meyer, Derek and Vega Yon, George (2023). epiworldR: Fast Agent-Based
 #>   Epi Models. Journal of Open Source Software, 8(90), 5781,
 #>   https://doi.org/10.21105/joss.05781
-#>
+#> 
 #> And the actual R package:
-#>
-#>   Meyer D, Vega Yon G (????). _epiworldR: Fast Agent-Based Epi Models_.
-#>   R package version 0.3-2, <https://github.com/UofUEpiBio/epiworldR>.
-#>
+#> 
+#>   Meyer D, Pulsipher A, Vega Yon G (2024). _epiworldR: Fast Agent-Based
+#>   Epi Models_. R package version 0.6.0.0,
+#>   <https://github.com/UofUEpiBio/epiworldR>.
+#> 
 #> To see these entries in BibTeX format, use 'print(<citation>,
 #> bibtex=TRUE)', 'toBibtex(.)', or set
 #> 'options(citation.bibtex.max=999)'.

epiworldR.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: b00320e3-afb5-47c8-8486-cc8d2715c18b
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

inst/CITATION

@@ -1,4 +1,4 @@
-year <- sub("-.*", "", meta$Date)
+year <- 2024
 note <- sprintf("R package version %s", meta$Version)
 
 bibentry(
@@ -28,6 +28,7 @@ bibentry(bibtype = "Manual",
          title = "{{epiworldR: Fast Agent-Based Epi Models}}",
          author = c(
           person("Derek", "Meyer", comment = c(ORCID = "0009-0005-1350-6988")),
+          person(given="Andrew", family="Pulsipher", comment = c(ORCID = "0000-0002-0773-3210")),
           person("George", "Vega Yon", comment = c(ORCID = "0000-0002-3171-0844"))
           ),
          year = year,

vignettes/likelihood-free-mcmc.Rmd

@@ -29,10 +29,10 @@ In this example, we use LFMCMC to recover the parameters of an SIR model.
 Our SIR model will have the following characteristics:
 
 * **Virus Name:** COVID-19
-* **Initial Virus Prevalence:** 0.1
-* **Recovery Rate:** 0.3
-* **Transmission Rate:** 0.3
-* **Number of Agents:** 1,000
+* **Initial Virus Prevalence:** 0.01
+* **Recovery Rate:** 1/7 (0.14)
+* **Transmission Rate:** 0.1
+* **Number of Agents:** 2,000
 
 We use the `ModelSIR` and `agents_smallworld` functions to construct the model in epiworldR.
 
@@ -43,14 +43,14 @@ model_seed <- 122
 
 model_sir <- ModelSIR(
   name = "COVID-19",
-  prevalence = .1,
-  transmission_rate = .3,
-  recovery_rate = .3
+  prevalence = .01,
+  transmission_rate = .1,
+  recovery_rate = 1/7
 )
 
 agents_smallworld(
   model_sir,
-  n = 1000,
+  n = 2000,
   k = 5,
   d = FALSE,
   p = 0.01
@@ -114,16 +114,21 @@ The **proposal function** returns a new set of parameters, which it is "proposin
 
 ```{r propfun}
 proposal_fun <- function(old_params, lfmcmc_obj) {
-  res <- plogis(qlogis(old_params) + rnorm(length(old_params)))
+  res <- plogis(qlogis(old_params) + rnorm(length(old_params), sd = .1))
   return(res)
 }
 ```
 
 The **kernel function** effectively scores the results of the latest simulation run against the observed data, by comparing the summary statistics from `summary_fun` for each. LFMCMC uses the kernel score and the Hastings Ratio to determine whether or not to accept the parameters for that run. In our example, since `summary_fun` simply passes the data through, `simulated_stats` and `observed_stats` are our simulated and observed data respectively.
 
 ```{r kernfun}
-kernel_fun <- function(simulated_stats, observed_stats, epsilon, lfmcmc_obj) {
-  dnorm(sqrt(sum((simulated_stats - observed_stats)^2)))
+kernel_fun <- function(
+    simulated_stats, observed_stats, epsilon, lfmcmc_obj
+    ) {
+  
+  diff <- ((simulated_stats - observed_stats)^2)^epsilon
+  dnorm(sqrt(sum(diff)))
+  
 }
 ```
 
@@ -140,16 +145,14 @@ lfmcmc_model <- LFMCMC(model_sir) |>
 
 # Run LFMCMC Simulation
 
-To run LFMCMC, we need to set the initial model parameters. For our example, we use an initial Recovery rate of 0.1 and an initial Transmission rate of 0.5. We set the kernel epsilon to 1.0 and run the simulation for 2,000 samples (iterations) using the `run_lfmcmc` function.
+To run LFMCMC, we need to set the initial model parameters. For our example, we use an initial Recovery rate of 0.3 and an initial Transmission rate of 0.3. We set the kernel epsilon to 1.0 and run the simulation for 2,000 samples (iterations) using the `run_lfmcmc` function.
 
 ```{r lfmcmc-run}
-initial_params <- c(0.1, 0.5)
+initial_params <- c(0.3, 0.3)
 epsilon <- 1.0
 n_samples <- 2000
 
 # Run the LFMCMC simulation
-verbose_off(lfmcmc_model)
-
 run_lfmcmc(
   lfmcmc = lfmcmc_model,
   params_init = initial_params,
@@ -160,13 +163,40 @@ run_lfmcmc(
 ```
 
 # Results
-To make the printed results easier to read, we use the `set_params_names` and `set_stats_names` functions before calling `print`. We also use a burn-in period of 500 samples.
+To make the printed results easier to read, we use the `set_params_names` and `set_stats_names` functions before calling `print`. We also use a burn-in period of 1,500 samples.
 
 ```{r lfmcmc-print}
 set_params_names(lfmcmc_model, c("Recovery rate", "Transmission rate"))
 set_stats_names(lfmcmc_model, get_states(model_sir))
 
-print(lfmcmc_model, burnin = 500)
+print(lfmcmc_model, burnin = 1500)
+```
+
+We can also look at the trace of the parameters:
+
+```{r post-dist}
+# Extracting the accepted parameters
+accepted <- get_all_accepted_params(lfmcmc_model)
+
+# Plotting the trace
+plot(
+  accepted[, 1], type = "l", ylim = c(0, 1),
+  main = "Trace of the parameters",
+  lwd = 2,
+  col = "tomato",
+  xlab = "Step",
+  ylab = "Parameter value"
+  )
+
+lines(accepted[, 2], type = "l", lwd = 2, col = "steelblue")
+
+legend(
+  "topright",
+  bty = "n",
+  legend = c("Recovery rate", "Transmission rate"),
+  pch    = 20,
+  col    = c("tomato", "steelblue")
+)
 ```
 
 Recall that the observed data came from a model with a Recovery rate of 0.3 and a Transmission rate of 0.3. As the above output shows, LFMCMC made a close approximation of the parameters, which resulted in a close approximation of the observed population distribution. This example highlights the effectiveness of using LFMCMC for highly complex models.