Disease Burden: Part I

In the previous two labs you were able to measure water and air quality in a laboratory setting. We will now use some of the data collected from those experiments to make inferences about the burden of disease associated with varying levels of water and air quality.

For the first part of this lab, you will attempt to answer a simple question: how would the burden of disease be different for the United States if drinking water quality was similar to what we tested in the Willamette River? To answer this question, you will need to break it down into several sub-questions, including:

1. What is the current burden of disease associated with drinking water illnesses in the United States?

2. How would the burden of disease change based on the quality of water that we tested for the water quality lab?

   1. How many people would be affected by the reduced drinking water quality?

   2. How much water do these people drink on average each day?

   3. How many bacteria would they be ingesting, and how many of them would get sick as a result of ingesting more pathogens?

   4. What would be the total health impact, expressed in disability adjusted life years (DALYs)?

In other words, you will be conducting a quantitative microbial risk assessment (QMRA) to estimate the increased burden of disease associated with reduced drinking-water quality. We recommend that you use the following resources to aid in your investigation:

1. Havelaar, A.H., & Melse, J.M. (2003) Quantifying public health risk in the WHO Guidelines for drinking-water quality: A burden of disease approach. National Institute for Public Health and the Environment of the Netherlands
   • A discussion paper on the concepts and methodology of disability-adjusted life years (DALYs) as a common public health metric and its usefulness for drinking-water quality. http://www.who.int/water_sanitation_health/publications/quantifyinghealthrisks/en/
2. FAO & WHO. (2003) Hazard characterization for pathogens in food and water: Guidelines.
   • A practical framework and structured approach for the characterization of microbial hazards in food and water, to assist governmental and research scientists. http://www.who.int/water_sanitation_health/publications/hazard-characterization-forpathogens/ en/
3. Brown, J., & Clasen, T. (2012). High Adherence Is Necessary to Realize Health Gains from Water Quality Interventions. PLoS ONE, 7(5), e36735. http://doi.org/10.1371/journal.pone.0036735
   • A description of a QRMA model used to estimate the level of adherence necessary to significantly increase DALYs using different levels of contamination and water treatment efficacies.

1. Current Burden of Disease

To determine the current burden of disease from unsafe water, you’ll want to check out the Global Burden of Disease (GBD) Visualization at the Institute of Health Metrics and Evaluation (IHME): https://vizhub.healthdata.org/gbd-compare/.

Select the second icon on the left-hand side (“Risks by Cause”), select level four on the slider, and select “United States” for location. Then select “Diarrhea/LRI/other” on the right-hand side. Find “Unsafe water” in the rankings. What are the total number of DALYs, the percent of total DALYs, and the DALY rate (number of DALYs per 100,000)?

Now navigate to the arrow diagram showing the risks for each cause and make sure that level four is still selected on the slider and United States is selected for location. Download this data using the Download button on the top right and determine where diarrheal diseases ranked for the United States in 2016.

Now do the same thing (determine the total DALYs, the percent of total DALYs, the DALY rate, and the cause ranking) for Sub-Saharan Africa in 2016.

Screenshot of IHME Visualization Dashboard

2a. Population

What was the population of the United States in 2016?

library(MonteCarlo)
library(data.table)

population = 0

2b. Water Consumption

How much water do people drink on average per day? Should the amount of water consumed be averaged or simulated? What are the benefits of conducting a simulation versus a single estimate based on average values?

If you’re not familiar with generating random numbers or writing a function in R, here’s an example to simulate the approximate number of liters that a person could drink each day in a year. See if you can generate a data table that randomly samples the amount of water 100,000 people would drink in a year. Hint: you might try the function replicate.

# function to randomly selected between 1, 2, and 5 liters of water per person per day
generate_water = function(n=1) {
  sample(c(1,2,5),n,replace = TRUE)
}
generate_water(365)

  [1] 1 1 2 2 2 2 1 1 5 5 1 5 2 5 1 2 1 1 5 1 1 1 5 1 1 5 2 1 1 1 1 5 5 1 2 2 5 2 1 5 2 5 2 2 1
 [46] 2 2 2 5 5 1 5 5 5 1 5 5 5 1 1 1 2 5 5 1 1 1 1 5 2 2 2 1 2 1 1 1 1 5 2 1 2 5 2 5 5 1 1 2 1
 [91] 5 2 1 1 5 1 2 2 1 2 1 2 5 2 1 1 1 2 1 1 2 2 1 2 2 1 5 2 5 5 1 1 5 1 1 2 1 5 2 1 1 1 1 1 2
[136] 2 2 1 1 1 5 5 5 1 2 5 5 1 1 1 2 1 2 5 2 5 1 1 2 5 2 2 5 2 1 5 1 1 5 5 5 1 1 1 2 2 5 2 1 1
[181] 5 2 2 2 2 2 5 2 2 2 2 1 5 2 2 1 2 5 2 5 1 2 2 1 5 1 5 2 1 1 5 1 2 5 2 5 5 2 1 1 2 1 2 1 1
[226] 2 2 1 2 1 2 1 5 5 5 2 5 5 2 2 2 5 5 2 2 2 2 2 2 1 5 5 2 2 2 2 1 2 5 1 2 5 5 1 2 5 2 2 1 5
[271] 5 2 2 2 5 5 1 2 5 2 1 1 2 5 1 2 5 1 1 1 1 2 1 2 2 5 5 2 5 1 5 1 2 1 5 2 5 1 5 2 5 2 1 5 2
[316] 2 5 5 1 1 1 1 1 2 1 5 1 2 5 1 5 2 1 5 1 2 5 1 2 1 1 2 5 2 1 1 2 5 5 1 2 2 2 1 1 5 1 5 5 2
[361] 2 5 5 2 1

# generate a data table (or data frame or matrix) with 365 columns and 100,000 rows with randomly sample amounts of water for each entry

2c. Dose-Response

How many organisms would people be ingesting on average per liter of water? What is the probability of infection for the number of microbes ingested? What is the probability of illness given infection?

Here are two different ways to randomly select the number of organisms present in a liter of water: one samples from a log-normal distribution and the other samples from an empirical distribution based on water quality testing. Make sure that you can justify which one that you use, as well as the parameters used (e.g., the mean for the log-normal distribution). How do your two data tables compare?

# generate number of organisms consumed empirically (based on log-normal distribution)
organisms_lnorm <- function(n = 1,mean = 1) {
  
  rlnorm(n, meanlog = mean, sdlog = mean)
}
# generate number of organisms consumed empirically (based on E. coli water testing)
organisms_empirical <- function(n = 1) {
 
  sample(c(14,75,9,18,14,22,8,14,6,12,26,6,1,2,8,4),n,replace = TRUE) * 10
}
# now create two new data tables with the number of organisms per liter for 100,000 people for both distributions

You now have two data tables that represent the number of organisms 100,000 people would ingest in each liter of water for a year. Pick the one that you think is more accurate. Now you need to determine each person’s daily dose of organisms, their probability of infection, and their probability of illness given infection.

# number of microbes ingested per day
dose = dt_water * dt_lnorm
# probabiliy density function for campylobacter infection
beta_poisson = function(dose, alpha = 0.145, n50 = 896) {
  
  1 - ( 1 + (dose / n50) * ( 2^(1/alpha) - 1 ) )^(-alpha)
}
# generate probability of daily infection (beta-poisson function)
p_inf <- beta_poisson(dose)
  
# generate probability of illness given infection (0.3 for campylobacter)
p_illness_inf = 0.3
p_illness <- p_inf * p_illness_inf

2d. Health Impact

You should have a large matrix with the risk of daily infection due to ingesting different amounts of campylobacter bacteria over the course of one year for 100,000 individuals. Now we will use the disease burden associated with campylobacteriosis (try saying that five times really fast) to calculate the total number of disability adjusted life years associated with increased exposure to campylobacter bacteria. To do this you have to make assumptions (that can be justified) about the proportion of the population that is susceptible to infection and the disease burden associated with campylobacteriosis.

# percent of the population susceptible to infection (campylobacter)
susceptible = 1.00

# disease burden per case (campylobacter)
db_case = 4.6e-3

# generate DALYs per year
DALYs_year = p_illness * db_case * susceptible
  
# generate total DALYs per 100,000
DALYs_total <- sum(DALYs_year)

DALYs_total

Simulation of Disease Burden

Now you have the total disease burden associated with increased exposure to campylobacter. However, there are two problems. First, we assumed that a person could have campylobacteriosis every day, whereas in reality infections and illnesses have average durations that are baked into the disease burden per case estimates. Second, we made assumptions about the number of campylobacter bacteria present in each liter of water when in reality we’re not sure how many bacteria might be present. Also, for some illnesses the disease burden per case isn’t a fixed number. As a result, we need to re-run the experiment with multiple repetitions to generate a confidence interval for our DALY estimates.

In order to simulate the disease burden associated with increased exposure to pathogens, we will conduct a Monte Carlo simulation. I’ve opted to use the MonteCarlo package (https://www.r-bloggers.com/introducing-the-montecarlo-package/), but you can feel free to use another method as you wish.

I am going to include the code for campylobacter below. Your responsibility is to understand what I’m doing below and then run simulations for rotavirus, shiga-toxin producing E. coli STEC 0157, and cryptosporidium based on parameters and methodologies you can find in the literature. Also, because we’re running multiple iterations, the simulation is modified to more closely approximate the methodology used in Brown & Clasen (cited above).

Following the instructions from the R Bloggers website, I will create the disease_burden function to calculate the total DALYs for 100,000 people (though slightly different from how we did it up above). I then specify the parameters I want to vary in the simulation. In this case I will vary the mean from the log-normal distribution that we’re using to generate random numbers of pathogens.

Campylobacter, Rotavirus, and E. Coli

# function to calculate disease burden
disease_burden = function(people, meanlog, alpha, n50, p_illness_inf, susceptible, db_case){
  
  # generate water consumed
  water_consumed <- sample(c(1,2,5),people,replace = TRUE)
  
  # generate number of organism consumed from log-normal distribution
  number_of_organisms <- rlnorm(people, meanlog = meanlog)
  
  # generate daily dose
  dose <- water_consumed * number_of_organisms
  
  # generate probability of daily infection (beta-poisson function)
  p_inf <- 1 - ( 1 + (dose / n50) * ( 2^(1/alpha) - 1 ) )^(-alpha)
  
  # generate probability of illness
  p_illness <- p_inf * p_illness_inf
  
  # generate yearly infection risk
  p_year = 1 - (1 - p_illness)^365
  
  # generate DALYs per year
  DALYs_year = p_year * db_case * susceptible
  
  # generate total DALYs per 100,000
  DALYs_total <- sum(DALYs_year)
  
  # returns value for MonteCarlo function
  return(list("DALYs"=DALYs_total))
  
}

# parameter list for campylobacter
param_campylobacter <- list("people" = 100000,
                            "meanlog" = c(1, 0.1, 0.01, 0.001, 0.0001),
                            "alpha" = 0.145,
                            "n50" = 896,
                            "p_illness_inf" = 0.3,
                            "susceptible" = 1.00,
                            "db_case" = 4.6e-3)

# function to run the simulation
campylobacter <- MonteCarlo(func=disease_burden, nrep=1000, param_list = param_campylobacter)

# summary of results
summary(campylobacter$results$DALYs)



# parameter list for rotavirus (fill this in)


# parameter list for E. coli (fill this in)




# function to run the simulation
rotavirus <- MonteCarlo(func=disease_burden, nrep=1000, param_list = param_rotavirus)

# summary of results
summary(rotavirus$results$DALYs)

# function to run the simulation
ecoli <- MonteCarlo(func=disease_burden, nrep=1000, param_list = param_ecoli)

# summary of results
summary(ecoli$results$DALYs)

Cryptosporidium

Why is cryptosporidium being modeled in its own section?

disease_burden_crypto = function(people, meanlog, infectivity, p_illness_inf, susceptible, db_case){

  # fill this in
  
}

# parameter list for cryptosporidium (fill this in)


# function to run the simulation
crypto <- MonteCarlo(func=disease_burden_crypto, nrep=1000, param_list = param_crypto)

# summary of results
summary(crypto$results$DALYs)

You should now have estimates of the total DALYs due to infections from campylobacter, rotavirus, shiga-toxin producing E. coli STEC 0157, and cryptosporidium. How do these estimates compare with the first DALY estimates from 2a. through 2d. above? What different assumptions were made, and how did that affect the outcome?

Lab Assignment

You will need to compose a formal lab report (that includes part one and part two) not exceeding five pages (excluding appendices). This can be submitted as an R Notebook (PDF version of the main report less than 5 pages) or as a PDF. You will need to include a brief introduction, a methods section, results, discussion, and a brief conclusion. See the syllabus and grading rubric for more a more detailed description.

In particular, you should address the following questions:

1. How does the newly calculated disease burden compare to the previous disease burden? How does it compare with other causes of disease burden in the United States (https://vizhub.healthdata.org/gbd-compare/)? How about in Low SDI countries?

2. Simulations require a lot of assumptions. Comment on which assumptions felt dubious (or that you weren’t able to justify from the literature) and which assumptions most significantly influenced your results. How much faith do you have in your results?

3. Why is it important to conduct a Monte Carlo simulation instead of just using average values?

4. Is disease burden per case the same for all locations? How would that assumption affect the outcome of this experiment?

The final report will be need to be submitted on D2L through Dropbox by 3pm on Tuesday, March 13th.