|
| 1 | +# nolint start |
| 2 | + |
| 3 | +# Practical 2 -miscellaneous on final size of heterogeneous population |
| 4 | +# Activity 1 |
| 5 | + |
| 6 | +room_number <- 1 |
| 7 | + |
| 8 | + # Load packages ---------------------------------------------------------- |
| 9 | +library(finalsize) |
| 10 | +library(tidyverse) |
| 11 | +library(socialmixr) |
| 12 | + |
| 13 | +# Declare the value of the given R_0 --------------------------------------- |
| 14 | +r0 <- 1.8 |
| 15 | + |
| 16 | +# load the polymod survey ---------------------------------------- |
| 17 | +polymod <- socialmixr::polymod |
| 18 | + |
| 19 | +# Extract data for the country and age groups assigned to your room-------- |
| 20 | + |
| 21 | +# get your country polymod data |
| 22 | +contact_data <- socialmixr::contact_matrix( |
| 23 | + polymod, |
| 24 | + countries = "Italy", |
| 25 | + age_limits = c(0,30,50,80), |
| 26 | + symmetric = TRUE |
| 27 | +) |
| 28 | + |
| 29 | +# view the elements of the contact data list |
| 30 | +# the contact matrix |
| 31 | +contact_data$matrix |
| 32 | + |
| 33 | +# the demography data |
| 34 | +contact_data$demography |
| 35 | + |
| 36 | +# get the contact matrix and demography data |
| 37 | +# contact_matrix <- contact_data$matrix |
| 38 | +demography_vector <- contact_data$demography$population |
| 39 | +demography_data <- contact_data$demography |
| 40 | + |
| 41 | + |
| 42 | +# Transpose and normalize contact matrix--------------------------------------- |
| 43 | + |
| 44 | +# scale the contact matrix so the largest eigenvalue is 1.0 |
| 45 | +# this is to ensure that the overall epidemic dynamics correctly reflect |
| 46 | +# the assumed value of R0 |
| 47 | +contact_matrix <- t(contact_data$matrix) |
| 48 | +scaling_factor <- 1 / max(eigen(contact_matrix)$values) |
| 49 | +normalised_matrix <- contact_matrix * scaling_factor |
| 50 | + |
| 51 | +normalised_matrix |
| 52 | + |
| 53 | +# divide each row of the contact matrix by the corresponding demography |
| 54 | +# this reflects the assumption that each individual in group {j} make contacts |
| 55 | +# at random with individuals in group {i} |
| 56 | +divided_matrix <- normalised_matrix/demography_vector |
| 57 | + |
| 58 | +# Create susceptibility, and demography-susceptibility distribution matrices--- |
| 59 | +n_demo_grps <- length(demography_vector) |
| 60 | + |
| 61 | +# Number of susceptible groups |
| 62 | +n_susc_groups <- 1 |
| 63 | +# susceptibility level |
| 64 | +susc_guess <- 1 |
| 65 | + |
| 66 | +# Declare susceptibility matrix |
| 67 | + |
| 68 | +susc_uniform <- matrix( |
| 69 | + data = susc_guess, |
| 70 | + nrow = n_demo_grps, |
| 71 | + ncol = n_susc_groups |
| 72 | +) |
| 73 | + |
| 74 | +# Declare demography-susceptibility distribution matrix |
| 75 | +p_susc_uniform <- matrix( |
| 76 | + data = 1, |
| 77 | + nrow = n_demo_grps, |
| 78 | + ncol = n_susc_groups |
| 79 | +) |
| 80 | + |
| 81 | +# Calculate the final size ----------------------------------------------------- |
| 82 | + |
| 83 | +final_size_data <- final_size( |
| 84 | + r0 = r0, |
| 85 | + contact_matrix = divided_matrix, |
| 86 | + demography_vector = demography_vector, |
| 87 | + susceptibility = susc_uniform, |
| 88 | + p_susceptibility = p_susc_uniform |
| 89 | +) |
| 90 | + |
| 91 | +# View the output data frame |
| 92 | +final_size_data |
| 93 | + |
| 94 | +# Visualize the proportion infected in each demographic group------------------- |
| 95 | + |
| 96 | +# order demographic groups as factors |
| 97 | +final_size_data$demo_grp <- factor( |
| 98 | + final_size_data$demo_grp, |
| 99 | + levels = demography_data$age.group |
| 100 | +) |
| 101 | +# plot data |
| 102 | +ggplot(final_size_data) + |
| 103 | + geom_col( |
| 104 | + aes( |
| 105 | + demo_grp, p_infected |
| 106 | + ), |
| 107 | + colour = "black", fill = "grey" |
| 108 | + ) + |
| 109 | + scale_y_continuous( |
| 110 | + labels = scales::percent, |
| 111 | + limits = c(0, 1) |
| 112 | + ) + |
| 113 | + expand_limits( |
| 114 | + x = c(0.5, nrow(final_size_data) + 0.5) |
| 115 | + ) + |
| 116 | + theme_classic() + |
| 117 | + coord_cartesian( |
| 118 | + expand = FALSE |
| 119 | + ) + |
| 120 | + labs( |
| 121 | + x = "Age group", |
| 122 | + y = "% Infected" |
| 123 | + ) |
| 124 | +# Visualize the total number infected in each demographic group---------------- |
| 125 | +# prepare demography data |
| 126 | +demography_data <- contact_data$demography |
| 127 | + |
| 128 | +# merge final size counts with demography vector |
| 129 | +final_size_data <- merge( |
| 130 | + final_size_data, |
| 131 | + demography_data, |
| 132 | + by.x = "demo_grp", |
| 133 | + by.y = "age.group" |
| 134 | +) |
| 135 | + |
| 136 | +# reset age group order |
| 137 | +final_size_data$demo_grp <- factor( |
| 138 | + final_size_data$demo_grp, |
| 139 | + levels = contact_data$demography$age.group |
| 140 | +) |
| 141 | + |
| 142 | +# multiply counts with proportion infected |
| 143 | +final_size_data$n_infected <- final_size_data$p_infected * |
| 144 | + final_size_data$population |
| 145 | + |
| 146 | +final_size_data |
| 147 | + |
| 148 | +# nolint end |
0 commit comments