Comments updated and added

R/TZ_p_multiple.bites.R

@@ -1,16 +1,5 @@
-#' @title Calculate proportion of bite victims that recieved multiple bites
-#' @description This function produces a dataframe, giving the proportion of
-#' individuals that recieved mutliple bites during the bite incident, the numerator
-#' and denominator for the proportion, as well as 95% confidence intervals.
-#' @param dataframe the dataframe from which the proportion is calculated
-#' @param head_col (default="Head") the column stating a bite wound to the head
-#' @param arm_col (default="Arm") the column stating a bite wound to the arm
-#' @param hand_col (default="Hand") the column stating a bite wound to the hand
-#' @param trunk_col (default="Trunk") the column stating a bite wound to the trunk
-#' @param leg_col (default="Leg") the column stating a bite wound to the leg
-#' @param foot_col (default="Foot") the column stating a bite wound to the foot
-#' @return a dataframe providing the results of the calculations
-#' @export
+# This function gives information on individuals with multiple bite wounds from 
+# the same bite incident.
 TZ_p_multiple.bites <- function(dataframe, head_col="Head", arm_col="Arm", hand_col="Hand",
                                 trunk_col="Trunk", leg_col="Leg", foot_col="Foot"){
 

R/mixture_model.R

@@ -1,6 +1,8 @@
 # Load in data
 transmission <- read.csv("output/transmission.csv", stringsAsFactors=FALSE)
 
+# This function uses information on number of bite incidents, expoosures, no-PEP,
+# deaths and location on the body of bites.
 mixture_model<-function(n){ # Cleaveland 2002 mixture model
 
   N = transmission$n_rabid.bites # N BITES
@@ -52,6 +54,7 @@ mixture_model<-function(n){ # Cleaveland 2002 mixture model
                     p_death_leg=prob_death_given_leg*prob_leg))
 }
 
+# This function calculates quantiles for confidence of bite location and deaths
 get_quantiles <- function(model_results){
 
   # Create a vector of probability of bite site for each bite location (+ CIs)

ms.TZ_figures.R

@@ -1,7 +1,5 @@
-################################################################################
-#                            PRODUCE FIGURES/TABLES                            #
-################################################################################
-
+# This script will produce all figures and tables given in the manuscript, saved
+# in the figs/ and output/paper/ folders respectively.
 rm(list=ls())
 
 # Load in libraries
@@ -21,42 +19,45 @@ fig_palette <- c("#0072B2", "#E69F00")
 # 1. Load data #
 ################
 
+# Load data from output folder
 CT_data <- read.csv("output/CT_data.csv", stringsAsFactors=FALSE)
 CT_data_rabid <- read.csv("output/CT_data_rabid.csv", stringsAsFactors=FALSE)
 HC_bite_summary <- read.csv("output/HC_bite_summary_patients.csv", stringsAsFactors=FALSE)
 infect_results <- read.csv("output/p_infect.csv", stringsAsFactors=FALSE)
 model_results <- read.csv("output/p_infect_model.csv", stringsAsFactors=FALSE)
 HC_data <- read.csv("output/HC_data.csv", stringsAsFactors=FALSE)
 
-# Limit CT to Serengeti & Ngorongoro
+# Limit contact tracing data to Serengeti & Ngorongoro
 ser_ngor <- CT_data[which(CT_data$District %in% c("Serengeti", "Ngorongoro")),]
 ser_ngor_rabid <- ser_ngor[which(ser_ngor$Rabid=="Yes"),]
 
 ###################################################################
 # Figure 1: Variation in annual incidence of patient presentation #
 ###################################################################
 
-# Set the df order and District factor by dog population
+# Set the order and District factor by dog population
 boxplot_data <- HC_bite_summary[order(HC_bite_summary$Dogs),]
 row.names(boxplot_data) <- NULL
 levsord <- unique(boxplot_data$District)
 boxplot_data$District <- factor(boxplot_data$District, levels = levsord, ordered = TRUE)
 boxplot_data <- boxplot_data[with(boxplot_data, order(HDR)), ]
 boxplot_data$District <- factor(boxplot_data$District, levels = unique(boxplot_data$District))
 
-# Produce Figure 1 - Option 1 (District on x-axis)
-pdf("figs/Bite_incidence_Boxplot_1.pdf", width=10, height=8)
+# Produce Figure 1
+pdf("figs/Bite_incidence_Boxplot.pdf", width=10, height=8)
 ggplot(boxplot_data, aes(x=District, y=bites_per_pop_per_100000, fill=Setting)) +
   theme_classic()+
   geom_boxplot(position=position_dodge(0.9)) +
   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
   geom_point(shape=21, size=2, position=position_dodge(0.9)) +
-  stat_summary(fun.y=mean, geom="point", size = 2, colour="black",position=position_dodge(0.9)) +
+  stat_summary(fun.y=mean, geom="point", size = 2, colour="black", position=position_dodge(0.9)) +
   labs(x = "District", y="Incidence of bite patients per 100,000") +
   scale_fill_manual(name="Setting", values=fig_palette) + ylim(0,165) +
-  geom_segment(aes(x=1, y=160, xend=length(unique(boxplot_data$District)), yend=160), arrow = arrow(ends="both", length=unit(0.5, "cm"))) +
+  geom_segment(aes(x=1, y=160, xend=length(unique(boxplot_data$District)), yend=160), 
+               arrow = arrow(ends="both", length=unit(0.5, "cm"))) +
   annotate("text", x=3.5, y=155, label=paste0("Low HDR (", round(min(boxplot_data$HDR), digits=1), ")")) +
-  annotate("text", x=length(unique(boxplot_data$District))-2.5, y=155, label=paste0("High HDR (", round(max(boxplot_data$HDR), digits=1), ")")) #+
+  annotate("text", x=length(unique(boxplot_data$District))-2.5, y=155, 
+           label=paste0("High HDR (", round(max(boxplot_data$HDR), digits=1), ")"))
 dev.off()
 
 ###############################################################
@@ -75,15 +76,23 @@ final_infect_df <- data.frame(Bite_loc = infect_results$Bite.loc,
                               P_bite = quantiles_list$p_bite_site,
                               N_rabid_bites = infect_results$n_bite.from.rabid.animal,
                               stringsAsFactors = FALSE)
+
+# Use function to calculate quantiles
 quantiles_list <- get_quantiles(model_results)
+
+# Create new row for totals
 new_row <- data.frame(Bite_loc = "", Order = "", P_death = quantiles_list$p_infect_overall,
                       N_death = sum(final_infect_df$N_death),
                       N_no_pep = sum(final_infect_df$N_no_pep), P_bite = "",
                       N_rabid_bites = sum(final_infect_df$N_rabid_bites),
                       stringsAsFactors = FALSE)
 final_infect_df <- rbind(final_infect_df, new_row)
+
+# Re-order dataframe 
 final_infect_df <- final_infect_df[order(final_infect_df$Order, decreasing=TRUE),]
 final_infect_df$Order <- NULL
+
+# Update column names
 colnames(final_infect_df) <- c("Bite location",
                                "Probability of death (95% CI)", "Number of deaths",
                                "Number of bite victims that did not receive PEP",
@@ -96,16 +105,23 @@ write.csv(final_infect_df, "output/paper/Probablity_of_death.csv", row.names=FAL
 # Table 2: Details of human deaths where late/incomplete PEP was received #
 ###########################################################################
 
-# Create df of details on which bite victims died following late/incomplete pep
+# Subset data for individuals that died of Rabies, with either incomplete or late PEP
 deaths_details <- CT_data_rabid[which(CT_data_rabid$Patient.outcome=="Died"),]
 deaths_details <- deaths_details[which(deaths_details$complete=="incomplete" |  deaths_details$late.PEP=="late"),]
 
-deaths_details <- data.frame(PEP_Failure=NA, Region=deaths_details$Region, Age=deaths_details$Age..in.years.,
-                             Sex=deaths_details$Sex, PEP=deaths_details$PEP.administration,
-                             n_pep_received=deaths_details$n.PEP.recieved, delay_to_PEP.1=deaths_details$delay.pep.1,
-                             n_bites=deaths_details$n.bite.locs, Head=deaths_details$Head, Arm=deaths_details$Arm,
-                             Hand=deaths_details$Hand, Trunk=deaths_details$Trunk, Leg=deaths_details$Leg,
+# Create table
+deaths_details <- data.frame(PEP_Failure=NA, Region=deaths_details$Region, 
+                             Age=deaths_details$Age..in.years., Sex=deaths_details$Sex, 
+                             PEP=deaths_details$PEP.administration,
+                             n_pep_received=deaths_details$n.PEP.recieved, 
+                             delay_to_PEP.1=deaths_details$delay.pep.1,
+                             n_bites=deaths_details$n.bite.locs, Head=deaths_details$Head, 
+                             Arm=deaths_details$Arm, Hand=deaths_details$Hand, 
+                             Trunk=deaths_details$Trunk, Leg=deaths_details$Leg,
                              Foot=deaths_details$Foot, Bite_loc=NA)
+
+# Run loop to find where on the body each individual was bitten, saving results
+# into a vector and then pasting into single string
 for(i in 1:nrow(deaths_details)){
   bite_loc_string <- c()
   if(deaths_details$Head[i]==1){
@@ -129,7 +145,7 @@ for(i in 1:nrow(deaths_details)){
   deaths_details$Bite_loc[i] <- toString(bite_loc_string)
 }
 
-# Set PEP Failure & PEP route
+# Run loop to give detail on PEP Failure & PEP route (ID/IM)
 deaths_details$PEP <- as.character(deaths_details$PEP)
 for(i in 1:nrow(deaths_details)){
   if(deaths_details$n_pep_received[i]>2 & deaths_details$delay_to_PEP.1[i]>0){
@@ -147,6 +163,7 @@ for(i in 1:nrow(deaths_details)){
   }
 }
 
+# Select only required rows
 deaths_details <- dplyr::select(deaths_details, PEP_Failure, Age, Sex, n_pep_received,
                                 PEP, delay_to_PEP.1, n_bites, Bite_loc)
 
@@ -160,25 +177,26 @@ write.csv(deaths_details, "output/paper/Death_details.csv", row.names=FALSE)
 # Figure 2A: Delay between bite date & PEP 1 Hist #
 ##################################################
 
-# Remove any remaining delays outside of range
+# Remove any delays outside of range (i.e. outliers)
 summary(CT_data_rabid$delay.pep.1)
 barplot_df <- CT_data_rabid[which(CT_data_rabid$delay.pep.1>=0 & CT_data_rabid$delay.pep.1<91),]
 
-# Read any "Other" patient outcomes as "Fine"
+# Change "Other" patient outcomes to "Fine"
 barplot_df$Patient.outcome[which(barplot_df$Patient.outcome=="Other")] <- "Fine"
 
 # Seperate off NTz and STz
 NTz_barplot <- barplot_df[which(barplot_df$Study.Area=="NTz"),]; nrow(NTz_barplot)
 STz_barplot <- barplot_df[which(barplot_df$Study.Area %in% c("STz", "Pemba")),]; nrow(STz_barplot)
 
-# NTZ: Create dataframe and transform to matrix
+# NTZ: Find counts per day delay, create dataframe and transform to matrix
 NTz_delays <- hist(NTz_barplot$delay.pep.1, breaks=c(seq(-1, 14, by=1), 21, 35, 64, 84))
 NTz_barplot.df <- data.frame(Delay=NTz_delays$breaks[2:20],
                              N_Fine=NTz_delays$counts,
                              Location=rep("NTz", 19))
 NTz_barplot.df$Prop = NTz_barplot.df$N_Fine/sum(NTz_barplot.df$N_Fine)
 NTz_barplot.df$Delay <- factor(NTz_barplot.df$Delay, levels=sort(unique(NTz_barplot.df$Delay)))
 
+# STZ: Find counts per day delay, create dataframe and transform to matrix
 STz_delays <- hist(STz_barplot$delay.pep.1, breaks=c(seq(-1, 14, by=1), 21, 35, 64, 84))
 STz_barplot.df <- data.frame(Delay=STz_delays$breaks[2:20],
                              N_Fine=STz_delays$counts,
@@ -189,12 +207,14 @@ STz_barplot.df$Delay <- factor(STz_barplot.df$Delay, levels=sort(unique(STz_barp
 # Bind together
 barplot.df <- rbind(NTz_barplot.df, STz_barplot.df)
 
+# Create plot
 pdf("figs/delay_in_pep_V3.pdf", width=10, height=6)
 ggplot(data=barplot.df, aes(x=Delay, y=Prop, fill=Location)) +
   geom_col(position=position_dodge(), color="black") + theme_classic() + ylim(0,1) +
   labs(x="Delay between exposure and initiation of PEP",
        y="Proportion that presented at a health facility") +
-  scale_fill_manual(name="", breaks=c("NTz", "STz"), labels=c("Serengeti & Ngorongoro", "Southern Tanzania"),
+  scale_fill_manual(name="", breaks=c("NTz", "STz"), 
+                    labels=c("Serengeti & Ngorongoro", "Southern Tanzania"),
                     values=fig_palette) +
   theme(legend.position="top",
         text = element_text(size=12),
@@ -208,32 +228,40 @@ dev.off()
 # Figure 2B - Barplot of proportion doses receieved #
 #####################################################
 
-# Barplots of patients that received PEP dose 1 (as a proportion, 95%; include
-# clinical signs as no PEP), then dose 2, then 3, then 4, then 5 - for Mobile
-# surveillance (panel 1) and for contact tracing rabies exposures (SD/ND where
-# people pay), panel 2.
+# Exclude 'patients' that arrived at hospital with Rabies symptoms
 HC_data_sub <- HC_data[which(HC_data$visit_status>0),]
+
+# Find number of patients where PEP was available upon arrival for treatment
 hc_visits <- HC_data_sub %>%
   group_by(visit_status) %>%
   summarise(n_bite_victims = length(which(is_antirabies_available=="yes")))
+
+# Calculate proportion of individuals that received PEP
 hc_visits$p_receive <- hc_visits$n_bite_victims/length(which(HC_data_sub$visit_status==1))
 hc_visits$PEP_available <- "yes"
 hc_visits$data <- "MS"
 
+# Subset for individuals in CT data that sought healthcare
 CT_barplot <- ser_ngor_rabid[which(ser_ngor_rabid$Sought==1),]
+
+# Find number of individuals that received each number of PEP doses
 ct_doses <- data.frame(visit_status=1:5,
                        n_bite_victims=c(length(which(CT_barplot$n.PEP.recieved>0)),
                                         length(which(CT_barplot$n.PEP.recieved>1 & CT_barplot$PEP.1==1)),
                                         length(which(CT_barplot$n.PEP.recieved>2 & CT_barplot$PEP.1==1)),
                                         length(which(CT_barplot$n.PEP.recieved>3 & CT_barplot$PEP.1==1)),
                                         length(which(CT_barplot$n.PEP.recieved>4 & CT_barplot$PEP.1==1))))
+
+# Calculate proportion of individuals that received PEP
 ct_doses$p_receive <- ct_doses$n_bite_victims/nrow(CT_barplot)
 ct_doses$PEP_available <- "yes"
 ct_doses$data <- "CT"
 
+# Bind dataframes together & set factor levels
 barplot_data <- rbind(hc_visits, ct_doses)
 barplot_data$data <- factor(barplot_data$data, levels = c("CT", "MS"))
 
+# Create plot
 pdf("figs/Proportion_PEP_completion.pdf", width=10, height=6)
 ggplot(data=barplot_data, aes(x=visit_status, y=p_receive, fill=data)) +
   geom_bar(stat="identity", position=position_dodge(), colour="black") +

ms.TZ_paper_stats.R

@@ -289,21 +289,28 @@ paste0("Lowest ideal average travelling distance ", round(mean(HC_distance_data$
 # Caption for figure 2 #
 ########################
 
-# Remove any remaining delays outside of range
-summary(traced_CT_data$delay.pep.1)
-hist_cap <- traced_CT_data[which(traced_CT_data$delay.pep.1>=0 & traced_CT_data$delay.pep.1<91),]
-
-paste0("Figure 2: Of ", length(which(hist_cap$late.PEP=="late")),
-       " patients with delayed PEP, ", length(which(hist_cap$late.PEP=="late" & hist_cap$Patient.outcome=="Died")),
-       " deaths occured.")
-
-########################
-# Caption for figure 3 #
-########################
+# Remove delays outside of range
+summary(CT_data_rabid$delay.pep.1)
+hist_cap <- CT_data_rabid[which(CT_data_rabid$delay.pep.1>=0 & CT_data_rabid$delay.pep.1<91),]
+hist_cap_traced <- traced_CT_data[which(traced_CT_data$delay.pep.1>=0 & traced_CT_data$delay.pep.1<91),]
 
+# Subset HC for individuals that received PEP
 HC_barplot_cap <- HC_data[which(HC_data$visit_status>0 & HC_data$is_antirabies_available=="yes"),]
+
+# Subset CT for uspect exposure individuals that sought PEP and did not die of Rabies
 CT_barplot_cap <- ser_ngor_rabid[which(ser_ngor_rabid$Sought==1 & ser_ngor_rabid$Patient.outcome!="Died"),]
 
-paste0("Dark grey indicates patients recorded in mobile phone-based surveillance data from Southern Tanzania (",
-       format(nrow(HC_barplot_cap), big.mark=","), ") and light grey indicates rabies exposed patients from Serengeti and Ngorongoro districts (",
-       format(nrow(CT_barplot_cap), big.mark=","), ")")
+paste0("Panel A shows contact tracing data on delays between exposure and initiation of PEP 
+       for rabies exposed persons ( ", length(which(hist_cap$Study.Area=="NTz")), 
+       " exposures from Serengeti and Ngorongoro districts, and ", 
+       length(which(hist_cap$Study.Area %in% c("STz", "Pemba"))), 
+       "exposures from 11 districts in southern Tanzania). Out of ", 
+       length(which(hist_cap_traced$late.PEP=="late")),
+       " patients identified through contact tracing who had delayed PEP (more than 1 day late), ",
+       length(which(hist_cap_traced$late.PEP=="late" & hist_cap_traced$Patient.outcome=="Died")),
+       " deaths occurred (Table 2). ")
+
+paste0("Panel B shows mobile phone-based surveillance records from Southern Tanzania (yellow, n = ",
+       format(nrow(HC_barplot_cap), big.mark=","),
+       ") of PEP completion and contact tracing data on rabies exposed patients from Serengeti and 
+       gorongoro districts (n = ", format(nrow(CT_barplot_cap), big.mark=","), "). ")