Skip to content

Commit

Permalink
Remove tidyr dependency from demography vignette
Browse files Browse the repository at this point in the history
  • Loading branch information
pwinskill committed Jul 26, 2022
1 parent 73022ed commit d973473
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 26 deletions.
1 change: 0 additions & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,6 @@ Suggests:
DiagrammeR,
cowplot,
ggplot2,
tidyr,
covr,
mgcv
RoxygenNote: 7.2.1
Expand Down
3 changes: 2 additions & 1 deletion R/variables.R
Original file line number Diff line number Diff line change
Expand Up @@ -406,11 +406,12 @@ calculate_initial_ages <- function(parameters) {
)))
}

deathrates <- parameters$deathrates[1, , drop = FALSE]
age_high <- parameters$deathrate_agegroups
age_width <- diff(c(0, age_high))
age_low <- age_high - age_width
n_age <- length(age_high)
birthrate <- find_birthrates(parameters$human_population, agegroups, deathrates[1,])
birthrate <- find_birthrates(parameters$human_population, age_high, deathrates)
deathrates <- parameters$deathrates[1,]

eq_pop <- get_equilibrium_population(age_high, birthrate, deathrates)
Expand Down
4 changes: 2 additions & 2 deletions tests/testthat/test-demography.R
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ test_that('find_birthrates is consistent with get_equilibrium_population', {


test_that('calculate_initial_ages calculates truncated exp custom demographic', {
parameters <- get_parameters()
parameters <- get_parameters(list(human_population = 4))
ages <- c(50, 100) * 365
deathrates <- c(.5, .75)
parameters <- set_demography(
Expand All @@ -50,7 +50,7 @@ test_that('calculate_initial_ages calculates truncated exp custom demographic',
mockery::mock(c(3, 1))
)
ages <- calculate_initial_ages(parameters)
mockery::expect_args(mock_groups, 1, 2, 100, replace = TRUE, prob = c(3, 1))
mockery::expect_args(mock_groups, 1, 2, 4, replace = TRUE, prob = c(3, 1))
mockery::expect_args(mock_rtexp, 1, 2, .5, 50 * 365)
mockery::expect_args(mock_rtexp, 2, 2, .75, 50 * 365)
expect_setequal(ages, c(25 * 365, 75 * 365, 30 * 365, 80 * 365))
Expand Down
60 changes: 38 additions & 22 deletions vignettes/Demography.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ knitr::opts_chunk$set(

```{r setup}
suppressPackageStartupMessages(library(ggplot2))
library(tidyr)
library(malariasimulation)
```

Expand All @@ -31,11 +30,14 @@ sim_length <- 5 * year
human_population <- 1000
starting_EIR <- 5
age_min <- seq(0, 80, 5) * 365
age_max <- seq(5, 85, 5) * 365
simparams <- get_parameters(
list(
human_population = human_population,
age_group_rendering_min_ages = seq(0, 80, 5) * 365,
age_group_rendering_max_ages = seq(5, 85, 5) * 365
age_group_rendering_min_ages = age_min,
age_group_rendering_max_ages = age_max
)
)
Expand Down Expand Up @@ -68,28 +70,44 @@ demography_params <- set_demography(
)
```

Let's run the simulations and compare the age distributions in the populations
in year 5:
Let's run the simulations

```{r}
# run and combine the outputs
exp_output <- run_simulation(sim_length, simparams)
exp_output$run <- 'exponential'
custom_output <- run_simulation(sim_length, demography_params)
custom_output$run <- 'custom'
output <- rbind(exp_output, custom_output)
```

and now we can compare the age distributions in the populations
in year 5:

```{r}
# wrangle outputs
output_long <- output |>
dplyr::filter(timestep == 5 * 365) |>
dplyr::select(run, timestep, contains("n_age")) |>
tidyr::pivot_longer(cols = -c(run, timestep), names_to = "age_group", values_to = "n", names_prefix = "n_age_") |>
tidyr::separate(col = age_group, sep = "_", into = c("age_lower", "age_upper"), convert = TRUE) |>
dplyr::mutate(age_mid = (age_lower + (age_upper - age_lower) / 2) / 365,
run = factor(run, levels = c("exponential", "custom")))
output <- rbind(exp_output, custom_output)
output <- output[output$timestep == 5 * 365,]
# A function to extract the age variables and convert to long format
convert_to_long <- function(age_min, age_max, output){
output <- lapply(
seq_along(age_min),
function(i) {
data.frame(
age_lower = age_min[[i]],
age_upper = age_max[[i]],
n = output[,paste0('n_age_', age_min[[i]], '_',age_max[[i]])],
age_mid = (age_min[[i]] + (age_min[[i]] - age_max[[i]]) / 2) / 365,
run = output$run,
timestep = output$timestep)
}
)
output <- do.call("rbind", output)
}
output <- convert_to_long(age_min, age_max, output)
# Plot the age distributions
ggplot(output_long, aes(x = age_mid, y = n)) +
ggplot(output, aes(x = age_mid, y = n)) +
geom_bar(stat = "identity") +
theme_bw() +
facet_wrap(~ run)
Expand All @@ -98,7 +116,7 @@ ggplot(output_long, aes(x = age_mid, y = n)) +

We can also specify time-varying death rates to capture a dynamic demography

```{r}
```{r, fig.width=7}
dynamic_demography_params <- simparams
# Set a flat demography
Expand All @@ -125,15 +143,13 @@ dynamic_demography_params <- set_demography(
dynamic_demography_output <- run_simulation(sim_length, dynamic_demography_params)
# wrangle outputs
dynamic_output_long <- dynamic_demography_output |>
dplyr::filter(timestep %in% (1:5 * 365)) |>
dplyr::select(timestep, contains("n_age")) |>
tidyr::pivot_longer(cols = -timestep, names_to = "age_group", values_to = "n", names_prefix = "n_age_") |>
tidyr::separate(col = age_group, sep = "_", into = c("age_lower", "age_upper"), convert = TRUE) |>
dplyr::mutate(age_mid = (age_lower + (age_upper - age_lower) / 2) / 365)
dynamic_demography_output <- dynamic_demography_output[dynamic_demography_output$timestep %in% (1:5 * 365),]
dynamic_demography_output$run <- "dynamic"
dynamic_demography_output <- convert_to_long(age_min, age_max, dynamic_demography_output)
# Plot the age distributions each year
ggplot(dynamic_output_long, aes(x = age_mid, y = n)) +
ggplot(dynamic_demography_output, aes(x = age_mid, y = n)) +
geom_bar(stat = "identity") +
theme_bw() +
facet_wrap(~ timestep / 365, ncol = 5)
Expand Down

0 comments on commit d973473

Please sign in to comment.