Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add epicontacts + fitdistrplus + superspreading #55

Open
avallecam opened this issue May 4, 2024 · 0 comments
Open

add epicontacts + fitdistrplus + superspreading #55

avallecam opened this issue May 4, 2024 · 0 comments
Labels
enhancement New feature or request

Comments

@avallecam
Copy link
Member

# activity 1 --------------------------------------------------------------

library(epicontacts)
library(tidyverse)

contacts <- read_rds("https://epiverse-trace.github.io/tutorials-middle/data/set-01-contacts.rds")
linelist <- read_rds("https://epiverse-trace.github.io/tutorials-middle/data/set-01-linelist.rds")

contacts
#> # A tibble: 137 × 3
#>     from    to  time
#>    <dbl> <dbl> <dbl>
#>  1     1     2  1.66
#>  2     1     3  3.24
#>  3     1     4  2.46
#>  4     1     5  2.98
#>  5     2     6  5.87
#>  6     2     7  6.19
#>  7     2     8  2.97
#>  8     2     9  6.64
#>  9     2    10  4.83
#> 10     2    11  4.53
#> # ℹ 127 more rows
linelist
#> # A tibble: 138 × 1
#>       id
#>    <dbl>
#>  1     1
#>  2     2
#>  3     3
#>  4     4
#>  5     5
#>  6     6
#>  7     7
#>  8     8
#>  9     9
#> 10    10
#> # ℹ 128 more rows

epi_contacts <- epicontacts::make_epicontacts(linelist = linelist,contacts = contacts)

# epicontacts::vis_epicontacts(epi_contacts)


# activity 2 --------------------------------------------------------------

# count secondary cases per infector
infector_secondary <- epi_contacts %>%
  pluck("contacts") %>%
  count(from, name = "secondary_cases") # --------------------- emphasis

# infector_secondary

all_secondary <- epi_contacts %>%
  # extract ids in contact *and* linelist using "which" argument
  epicontacts::get_id(which = "all") %>% # --------------------- emphasis
  # transform vector to dataframe to use left_join()
  enframe(name = NULL, value = "from") %>%
  # join count secondary cases per infectee
  left_join(infector_secondary) %>% # --------------------- emphasis
  # infectee with missing secondary cases are replaced with zero
  replace_na(
    replace = list(secondary_cases = 0)
  )
#> Joining with `by = join_by(from)`

## plot the distribution
all_secondary %>%
  ggplot(aes(secondary_cases)) +
  geom_histogram(binwidth = 1) +
  labs(
    x = "Number of secondary cases",
    y = "Frequency"
  )

# activity 3 --------------------------------------------------------------

library(fitdistrplus)
#> Loading required package: MASS
#> 
#> Attaching package: 'MASS'
#> 
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> 
#> Loading required package: survival

## fit distribution
offspring_fit <- all_secondary %>%
  pull(secondary_cases) %>%
  fitdist(distr = "nbinom") # --------------------- emphasis

offspring_fit
#> Fitting of the distribution ' nbinom ' by maximum likelihood 
#> Parameters:
#>       estimate  Std. Error
#> size 0.0246670 0.008351478
#> mu   0.9929577 0.544939632


# activity 4 --------------------------------------------------------------

library(superspreading)
#> Warning: package 'superspreading' was built under R version 4.3.3

proportion_cluster_size(
  R = 0.8,
  k = offspring_fit$estimate["size"],
  cluster_size = c(5,10,25) # --------------------- emphasis
)
#>     R        k prop_5 prop_10 prop_25
#> 1 0.8 0.024667  89.5%   77.2%   49.3%

Created on 2024-05-04 with reprex v2.1.0

@avallecam avallecam added the enhancement New feature or request label May 4, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
Status: No status
Development

No branches or pull requests

1 participant