include isolation of set 3

Author: avallecam

instructors/03-practical-tutors.qmd

@@ -243,87 +243,10 @@ Write your answers to the questions above:
 
 ##### Set 1 (sample)
 
-```{r}
-#| warning: false
-#| eval: false
+```{r, file = "fig/03-practical-instructor-1.R", eval = FALSE}
 
-# Load packages -----------------------------------------------------------
-library(epicontacts)
-library(fitdistrplus)
-library(tidyverse)
-
-
-# Read linelist and contacts ----------------------------------------------
-dat_contacts <- readr::read_rds(
-  "https://epiverse-trace.github.io/tutorials-middle/data/set-01-contacts.rds"  #<DIFFERENT PER GROUP>
-)
-
-dat_linelist <- readr::read_rds(
-  "https://epiverse-trace.github.io/tutorials-middle/data/set-01-linelist.rds"  #<DIFFERENT PER GROUP>
-)
-
-
-# Create an epicontacts object -------------------------------------------
-epi_contacts <-
-  epicontacts::make_epicontacts(
-    linelist = dat_linelist,
-    contacts = dat_contacts,
-    directed = TRUE
-  )
-
-epi_contacts
-
-# visualize the contact network
-contact_network <- epicontacts::vis_epicontacts(epi_contacts)
-
-contact_network
-
-
-# Count secondary cases per subject in contacts and linelist --------
-secondary_cases <- epicontacts::get_degree(
-  x = epi_contacts,
-  type = "out",
-  only_linelist = TRUE
-)
-
-# plot the histogram of secondary cases
-individual_reproduction_num <- secondary_cases %>%
-  enframe() %>%
-  ggplot(aes(value)) +
-  geom_histogram(binwidth = 1) +
-  labs(
-    x = "Number of secondary cases",
-    y = "Frequency"
-  )
-
-individual_reproduction_num
-
-
-# Fit a negative binomial distribution -----------------------------------
-offspring_fit <- secondary_cases %>%
-  fitdistrplus::fitdist(distr = "nbinom")
-
-offspring_fit
-
-
-# Estimate proportion of new cases from a cluster of secondary cases -----
-
-# Set seed for random number generator
-set.seed(33)
-
-# Estimate the proportion of new cases originating from 
-# a transmission cluster of at least 5, 10, or 25 cases
-proportion_cases_by_cluster_size <-
-  superspreading::proportion_cluster_size(
-    R = offspring_fit$estimate["mu"],
-    k = offspring_fit$estimate["size"],
-    cluster_size = c(5, 10, 25)
-  )
-
-proportion_cases_by_cluster_size
 ```
 
-
 #### Outputs
 
 Group 1
@@ -556,100 +479,10 @@ Write your answers to the questions above:
 
 ##### Set 1 (sample)
 
-```{r}
-#| warning: false
-#| eval: false
-
-# Load packages -----------------------------------------------------------
-library(epiparameter)
-library(epichains)
-library(tidyverse)
-
-
-# Set input parameters ---------------------------------------------------
-known_basic_reproduction_number <- 0.8 #<DIFFERENT PER GROUP>
-known_dispersion <- 0.01 #<DIFFERENT PER GROUP>
-chain_to_observe <- 957 #<DIFFERENT PER GROUP>
-
+```{r, file = "fig/03-practical-instructor-2.R", eval = FALSE}
 
-# Set iteration parameters -----------------------------------------------
-
-# Create generation time as <epiparameter> object
-generation_time <- epiparameter::epiparameter(
-  disease = "disease x",
-  epi_name = "generation time",
-  prob_distribution = "gamma",
-  summary_stats = list(mean = 3, sd = 1)
-)
-
-
-# Simulate multiple chains -----------------------------------------------
-# run set.seed() and epichains::simulate_chains() together, in the same run
-
-# Set seed for random number generator
-set.seed(33)
-
-multiple_chains <- epichains::simulate_chains(
-  # simulation controls
-  n_chains = 1000, # number of chains to simulate
-  statistic = "size",
-  stat_threshold = 500, # stopping criteria
-  # offspring
-  offspring_dist = rnbinom,
-  mu = known_basic_reproduction_number,
-  size = known_dispersion,
-  # generation
-  generation_time = function(x) generate(x = generation_time, times = x)
-)
-
-multiple_chains
-
-
-# Explore suggested chain ------------------------------------------------
-multiple_chains %>%
-  # use data.frame output from <epichains> object
-  as_tibble() %>%
-  filter(chain == chain_to_observe) %>%
-  print(n = Inf)
-
-
-# visualize ---------------------------------------------------------------
-
-# daily aggregate of cases
-aggregate_chains <- multiple_chains %>%
-  as_tibble() %>%
-  # count the daily number of cases in each chain
-  mutate(day = ceiling(time)) %>%
-  count(chain, day, name = "cases") %>%
-  # calculate the cumulative number of cases for each chain
-  group_by(chain) %>%
-  mutate(cumulative_cases = cumsum(cases)) %>%
-  ungroup()
-
-# Visualize transmission chains by cumulative cases
-aggregate_chains %>%
-  # create grouped chain trajectories
-  ggplot(aes(x = day, y = cumulative_cases, group = chain)) +
-  geom_line(color = "black", alpha = 0.25, show.legend = FALSE) +
-  # define a 100-case threshold
-  geom_hline(aes(yintercept = 100), lty = 2) +
-  labs(x = "Day", y = "Cumulative cases")
-
-# count chains over 100 cases
-aggregate_chains %>%
-  filter(cumulative_cases >= 100) %>%
-  count(chain)
-# distribution of size of chains
-aggregate_chains %>%
-  filter(cumulative_cases >= 100) %>% 
-  skimr::skim(cumulative_cases)
-# distribution of lenght of chains
-aggregate_chains %>%
-  filter(cumulative_cases >= 100) %>% 
-  skimr::skim(day)
 ```
 
-
 #### Outputs
 
 Group 1

instructors/fig/03-practical-instructor-1.R

@@ -0,0 +1,85 @@
+# nolint start
+
+# Practical 3
+# Activity 1
+
+# Load packages -----------------------------------------------------------
+library(epicontacts)
+library(fitdistrplus)
+library(tidyverse)
+
+
+# Read linelist and contacts ----------------------------------------------
+dat_contacts <- readr::read_rds(
+  "https://epiverse-trace.github.io/tutorials-middle/data/set-01-contacts.rds"  #<DIFFERENT PER GROUP>
+)
+
+dat_linelist <- readr::read_rds(
+  "https://epiverse-trace.github.io/tutorials-middle/data/set-01-linelist.rds"  #<DIFFERENT PER GROUP>
+)
+
+
+# Create an epicontacts object -------------------------------------------
+epi_contacts <- epicontacts::make_epicontacts(
+  linelist = dat_linelist,
+  contacts = dat_contacts,
+  directed = TRUE
+)
+
+# Print output
+epi_contacts
+
+# Visualize the contact network
+contact_network <- epicontacts::vis_epicontacts(epi_contacts)
+
+# Print output
+contact_network
+
+
+# Count secondary cases per subject in contacts and linelist --------------
+secondary_cases <- epicontacts::get_degree(
+  x = epi_contacts,
+  type = "out",
+  only_linelist = TRUE
+)
+
+# Plot the histogram of secondary cases
+individual_reproduction_num <- secondary_cases %>%
+  enframe() %>% 
+  ggplot(aes(value)) +
+  geom_histogram(binwidth = 1) +
+  labs(
+    x = "Number of secondary cases",
+    y = "Frequency"
+  )
+
+# Print output
+individual_reproduction_num
+
+
+# Fit a negative binomial distribution -----------------------------------
+offspring_fit <- secondary_cases %>%
+  fitdistrplus::fitdist(distr = "nbinom")
+
+# Print output
+offspring_fit
+
+
+# Estimate proportion of new cases from a cluster of secondary cases ------
+
+# Set seed for random number generator
+set.seed(33)
+
+# Estimate the proportion of new cases originating from 
+# a transmission cluster of at least 5, 10, or 25 cases
+proportion_cases_by_cluster_size <- 
+  superspreading::proportion_cluster_size(
+    R = offspring_fit$estimate["mu"],
+    k = offspring_fit$estimate["size"],
+    cluster_size = c(5, 10, 25)
+  )
+
+# Print output
+proportion_cases_by_cluster_size
+
+# nolint end

instructors/fig/03-practical-instructor-2.R

@@ -0,0 +1,94 @@
+# nolint start
+
+# Practical 3
+# Activity 2
+
+# Load packages -----------------------------------------------------------
+library(epiparameter)
+library(epichains)
+library(tidyverse)
+
+
+# Set input parameters ---------------------------------------------------
+known_basic_reproduction_number <- 0.8
+known_dispersion <- 0.01
+chain_to_observe <- 957
+
+
+# Set iteration parameters -----------------------------------------------
+
+# Create generation time as an <epiparameter> object
+generation_time <- epiparameter::epiparameter(
+  disease = "disease x",
+  epi_name = "generation time",
+  prob_distribution = "gamma",
+  summary_stats = list(mean = 3, sd = 1)
+)
+
+
+# Simulate multiple chains -----------------------------------------------
+# Run set.seed() and epichains::simulate_chains() together, in the same run
+
+# Set seed for random number generator
+set.seed(33)
+
+multiple_chains <- epichains::simulate_chains(
+  # Simulation controls
+  n_chains = 1000, # Number of chains to simulate
+  statistic = "size",
+  stat_threshold = 500, # Stopping criteria
+  # Offspring
+  offspring_dist = rnbinom,
+  mu = known_basic_reproduction_number,
+  size = known_dispersion,
+  # Generation
+  generation_time = function(x) generate(x = generation_time, times = x)
+)
+
+multiple_chains
+
+
+# Explore suggested chain ------------------------------------------------
+multiple_chains %>%
+  # Use data.frame output from <epichains> object
+  as_tibble() %>%
+  filter(chain == chain_to_observe) %>%
+  print(n = Inf)
+
+
+# Visualize --------------------------------------------------------------
+
+# Daily aggregate of cases
+aggregate_chains <- multiple_chains %>%
+  as_tibble() %>%
+  # Count the daily number of cases in each chain
+  mutate(day = ceiling(time)) %>%
+  count(chain, day, name = "cases") %>%
+  # Calculate the cumulative number of cases for each chain
+  group_by(chain) %>%
+  mutate(cumulative_cases = cumsum(cases)) %>%
+  ungroup()
+
+# Visualize transmission chains by cumulative cases
+aggregate_chains %>%
+  # Create grouped chain trajectories
+  ggplot(aes(x = day, y = cumulative_cases, group = chain)) +
+  geom_line(color = "black", alpha = 0.25, show.legend = FALSE) +
+  # Define a 100-case threshold
+  geom_hline(aes(yintercept = 100), lty = 2) +
+  labs(x = "Day", y = "Cumulative cases")
+
+# count chains over 100 cases
+aggregate_chains %>%
+  filter(cumulative_cases >= 100) %>%
+  count(chain)
+# distribution of size of chains
+aggregate_chains %>%
+  filter(cumulative_cases >= 100) %>% 
+  skimr::skim(cumulative_cases)
+# distribution of lenght of chains
+aggregate_chains %>%
+  filter(cumulative_cases >= 100) %>% 
+  skimr::skim(day)
+
+# nolint end