From d973473b5cb54fcbc70d5a7de0e672a5f8dbe3ba Mon Sep 17 00:00:00 2001 From: Peter Winskill Date: Tue, 26 Jul 2022 11:38:22 +0100 Subject: [PATCH] Remove tidyr dependency from demography vignette --- DESCRIPTION | 1 - R/variables.R | 3 +- tests/testthat/test-demography.R | 4 +-- vignettes/Demography.Rmd | 60 ++++++++++++++++++++------------ 4 files changed, 42 insertions(+), 26 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index bf6e19eb..e60e26c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -64,7 +64,6 @@ Suggests: DiagrammeR, cowplot, ggplot2, - tidyr, covr, mgcv RoxygenNote: 7.2.1 diff --git a/R/variables.R b/R/variables.R index 51197dfb..b12d1d8d 100644 --- a/R/variables.R +++ b/R/variables.R @@ -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) diff --git a/tests/testthat/test-demography.R b/tests/testthat/test-demography.R index 572706c2..3b08dd37 100644 --- a/tests/testthat/test-demography.R +++ b/tests/testthat/test-demography.R @@ -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( @@ -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)) diff --git a/vignettes/Demography.Rmd b/vignettes/Demography.Rmd index 20c0dac9..cfa68a00 100644 --- a/vignettes/Demography.Rmd +++ b/vignettes/Demography.Rmd @@ -16,7 +16,6 @@ knitr::opts_chunk$set( ```{r setup} suppressPackageStartupMessages(library(ggplot2)) -library(tidyr) library(malariasimulation) ``` @@ -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 ) ) @@ -68,8 +70,7 @@ 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 @@ -77,19 +78,36 @@ 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) @@ -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 @@ -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)