seantalts/covid_party_risk.Rmd

## covid_party_risk.Rmd
---
title: "Covid Party Risk - 10,000 Parties Edition"
author: "Sean Talts"
date: "8/1/2020"
output:  html_document
runtime: shiny
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

I'm not an epidemiologist! I attempt here to source all of the necessary epidemiological information from **you** via these inputs and, using very basic probability, see what that would imply about large gatherings with an isolation period beforehand. Please do not take these default values as meaningful; it's a very coarse approximation I wrote to deal with anxiety over an upcoming trip a group of friends we've been planning since the before-times.

Rather than assuming all of our inputs are correct, we're imagining throwing 10,000 parties with some uncertainty over all of our inputs. That means each of these graphs is a histogram of the various calculations over all of the 10,000 parties. Scroll down for more about the method and calculations. Source code is [here](https://gist.github.com/seantalts/f371ed854662ef185aff33275badc0d4). Please review!

## Base rate in your city
First we'll take a back-of-the-envelope stab at modeling the number of people in your city who are out spreading the disease. Take a look at adjusting the parameters according to your best guesses for your city based on recent news and then adjust the standard deviation sliders until the graphs for simulated values look like they reasonably capture your uncertainty about those numbers.
```{r echo=FALSE}
inputPanel(
  numericInput("days_infectious", label = "Number of days we think covid is contagious if you are infected, X", 10),
  numericInput("new_cases", label = "Average number of new cases per day for the past X days", 320),
  numericInput("city_pop", label = "Number of people who are living in your city", 7.5e6),
  sliderInput("asymptomatic_rate", label="Percentage of covid cases that are asymptomatic (and thus not tested)", 75,
              min=0, max=98),
  sliderInput("days_stdev", label="Days infected standard deviation (sd)", 1,
               min=0, max=6),
  sliderInput("cases_stdev", label="New cases stdev %", 20,
               min=2, max=100),
  sliderInput("city_pop_stdev", label="City population stdev %", 2,
               min=0, max=20),
  sliderInput("asym_rate_stdev", label="Asymptomatic rate stdev", 15,
               min=0, max=40)
)

bound_pct = function(draws) {
  n = length(draws)
  draws = draws[draws > 1 & draws < 99]
  missing = n - length(draws)
  bootstrapped = sample(draws, missing, replace=T)
  c(draws, bootstrapped) / 100
}

trunc_at_one = function(a) {
  a[a > 1] = 1
  a
}

get_draws = reactive({
  n = 10000
  draws_sdpct = function(mean, sd) {
    abs(rnorm(n, mean, sd / 100 * mean ))
  }
  draws = function(mean, sd) {
    abs(rnorm(n, mean, sd))
  }

  df = data.frame(
    days_infectious = draws(input$days_infectious, input$days_stdev),
    cases = draws_sdpct(input$new_cases, input$cases_stdev),
    city_pop = draws_sdpct(input$city_pop, input$city_pop_stdev),
    asym_rate = bound_pct(draws(input$asymptomatic_rate, input$asym_rate_stdev))
  )
  df$num_walking_infected = df$days_infectious * (df$cases / (1 - df$asym_rate) - df$cases)
  df$chance_random_person_infected = trunc_at_one(df$num_walking_infected / df$city_pop)
  df
})
# TODO: Change asymptomatic rate to "ascertainment error?" or break it out

renderText({
  df = get_draws3()
  paste(sep=" ", "Your average chance of being infected across", nrow(df),
        "simulations given all the data on this page is",
        round(mean(df$chance_you_contract) * 100, 3), "%")
})

fluidRow(
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$days_infectious,
                main="Histogram of # days infectious used in these simulations",
                xlab="Days",
                ylab="", yaxt='n')
         })),
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$cases,
                main="Histogram of number of new cases yesterday",
                xlab="Cases",
                ylab="", yaxt='n')
         })))

fluidRow(
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$city_pop,
                main="Histogram of city population",
                xlab="People",
                ylab="", yaxt='n')
         })),
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$asym_rate, breaks=100,
                main="Distribution of asymptomatic infection rates",
                xlab="Infection rate (of total infections)",
                ylab="", yaxt='n')
         })))


fluidRow(
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$num_walking_infected, breaks=100,
                main="Histogram of # asymptomatic infections at large",
                xlab="People infected",
                ylab="", yaxt='n')
         })),
  column(6,
         renderPlot({
           df = get_draws()
           hist(df$chance_random_person_infected, breaks=100,
                main="Histogram of the chance a random person is infected",
                xlab="Chance",
                ylab="", yaxt='n')
         })))
```

## Adding an isolation period
Here we'll try to model your friends interactions over the past week to see how likely each of them will be to be infected when they arrive. This lets us examine the effects of varying levels of isolation for the X days prior. These uncertainty parameters are pure standard deviations in the units of their associated parameters; transmission rate sd is in % chance of transission and friends of friends sd's unit is people.
```{r echo=FALSE}
inputPanel(
  sliderInput("transmission_rate", "Chance of transmission from hanging out with an infected person", 65,
              min=0, max=100),
  sliderInput("trans_rate_stdev", label="Chance of transmission uncertainty (sd)", 25,
               min=1, max=50),
  numericInput("transitive_guests", label="Number of people each guest has seen indoors in the past X days", 3),
  sliderInput("friends_of_friends_stdev", label="Friends-of-friends uncertainty (sd)", 3,
               min=1, max=50))

bootstrap_positive = function(draws) {
  n = length(draws)
  draws = draws[draws >= 0]
  missing = n - length(draws)
  bootstrapped = sample(draws, missing, replace=T)
  c(draws, bootstrapped)
}

get_draws2 = reactive({
  df = get_draws()
  n = nrow(df)
  df$trans_rate = bound_pct(rnorm(n, input$transmission_rate, input$trans_rate_stdev))
  df$friends_of_friends = bootstrap_positive(rnorm(n, input$transitive_guests, input$friends_of_friends_stdev))
  df$chance_friend_contracts = 1 - (1 - df$trans_rate * df$chance_random_person_infected)^df$friends_of_friends
  df
})

renderText({
  df = get_draws3()
  paste(sep=" ", "Your average chance of being infected across", nrow(df),
        "simulations given all the data on this page is",
        round(mean(df$chance_you_contract) * 100, 3), "%")
})

renderPlot({
  df = get_draws2()
  hist(df$friends_of_friends, breaks=-0.5:ceiling(max(df$friends_of_friends + 1)),
       main="Number of friends each guests hangs out with before your party",
       xlab="People",
       ylab="", yaxt='n')
})

fluidRow(
  column(6,
         renderPlot({
           df = get_draws2()
           hist(df$trans_rate, breaks=100,
                main="Distribution of transmission rates",
                xlab="Probability",
                ylab="", yaxt='n')
         })),
  column(6,
         renderPlot({
           df = get_draws2()
           hist(df$chance_friend_contracts, breaks=100,
                main="Histogram of the % chance any friend gets covid elsewhere",
                xlab="Probability",
                ylab="", yaxt='n')
         })))
```

## The party
```{r echo=FALSE}
inputPanel(
  numericInput("num_guests", label="Number of guests in your party", 17),
  sliderInput("stdev3", label="Standard deviation in guests", 5,
               min=2, max=50))

get_draws3 = reactive({
  df = get_draws2()
  draws = function(mean) {
    abs(rnorm(nrow(df), mean, input$stdev3 / 100 * mean ))
  }
  df$guests = draws(input$num_guests)
  df$num_friends_infected = rbinom(nrow(df), floor(df$guests), df$chance_friend_contracts)
  df$chance_you_contract = 1 - (1 - df$trans_rate)^df$num_friends_infected
  df
})

renderText({
  df = get_draws3()
  paste(sep=" ", "Your average chance of being infected across", nrow(df),
        "simulations given all the data on this page is",
        round(mean(df$chance_you_contract) * 100, 3), "%")
})

fluidRow(
  column(6,
         renderPlot({
           df = get_draws3()
           hist(df$guests, breaks=-0.5:ceiling(max(df$guests + 1)),
                main="Histogram of the number of guests in our simulations",
                xlab="People",
                ylab="", yaxt='n'
           )})),
  column(6,
         renderPlot({
           df = get_draws3()
           hist(df$num_friends_infected, breaks=-0.5:ceiling(max(df$num_friends_infected) + 1),
                main="Number of friends infected",
                xlab="People",
                ylab="", yaxt='n'
           )
         })
  ))
```

#### Detailed stat calculation simulation summaries
```{r echo=F}
renderTable(t(do.call(cbind, lapply(get_draws3(), summary))), rownames = T)
```

## Method
### Simulations
I wanted to do simulations that capture some uncertainty rather than using single estimated numbers for everything and plugging those into the formulas, so every number input is used as the mean of a normal distribution with the percent standard deviation input times the mean as the distribution's standard deviation.

In maths, $ x_{draws} \sim \mathcal{N}(x, \sigma x) $. For the number of days we think someone is infectious, with days_infectious = 10 and stdev = 10%, this looks like e.g. $ days.infectious \sim \mathcal{N}(10, 1) $ and generates this histogram,
```{r}
days_infectious = 10
stdev = 0.1 * days_infectious
hist(rnorm(10000, days_infectious, stdev))
```

### Formulæ
Here's the code, but let's walk through it. For non-R folks just ignore the `df$` bits.
```{r eval=F}
  df$num_walking_infected = df$days_infectious * (df$cases / (1 - df$asym_rate) - df$cases)
  df$chance_random_person_infected = trunc_at_one(df$num_walking_infected / df$city_pop)
  df$chance_friend_contracts = 1 - (1 - df$trans_rate * df$chance_random_person_infected)^df$friends_of_friends
  df$num_friends_infected = rbinom(n, floor(df$guests), df$chance_friend_contracts)
  df$chance_you_contract = 1 - (1 - df$trans_rate)^df$num_friends_infected
```

#### Number infected walking around your city
First we estimate the number of people walking around who are infected somewhat conservatively by saying that only those who are asymptomatic are out and about. If we think `cases = total * rate_of_symptoms` then we can find the total infected by dividing cases by rate of having symptoms (which is 1 - asymptomatic rate). So that gives us `total = cases / (1 - asym_rate)`, to find the number of people walking around spreading we can then just subtract the number of cases.

#### Number of people at your party infected
A binomial distribution gives us the number of "successes" in N trials with p probability of success. If we think about success as the chance of infection and each guest as an independent trial, we can use R's `rbinom(n, df$guests, df$chance_friend_contracts)` to model the number of guests who will have it. We have to convert the number of guests to an integer and also pass in n because we're doing everything as n simulations, so we'll end up with 10000 simulations each of which says how many people at your party are infected.

#### Given X infected people, what is your chance of infection?
The formula for at least one success in n trials with a probability p of success in each trial is `1 - (1 -p)^n`.  We can use the binomial distribution again to find the chance of infection given the number of people at your party who are infected and the transmission rate between people. that gives us `1 - (1 - df$trans_rate)^df$num_friends_infected`. We also use this to estimate the chance that each friend has contracted the virus based on their number of friends, the transmission rate, and the chance a random person is infected.

This model encodes a 2 step process - as if each friend has an "event" with infection risk from everyone they come into contact with in the past X days before they come to your party or new isolation pod. I did this because I wanted to model the effects of isolation prior to the event.
	---
	title: "Covid Party Risk - 10,000 Parties Edition"
	author: "Sean Talts"
	date: "8/1/2020"
	output: html_document
	runtime: shiny
	---

	```{r setup, include=FALSE}
	knitr::opts_chunk$set(echo = TRUE)
	```

	I'm not an epidemiologist! I attempt here to source all of the necessary epidemiological information from you via these inputs and, using very basic probability, see what that would imply about large gatherings with an isolation period beforehand. Please do not take these default values as meaningful; it's a very coarse approximation I wrote to deal with anxiety over an upcoming trip a group of friends we've been planning since the before-times.

	Rather than assuming all of our inputs are correct, we're imagining throwing 10,000 parties with some uncertainty over all of our inputs. That means each of these graphs is a histogram of the various calculations over all of the 10,000 parties. Scroll down for more about the method and calculations. Source code is [here](https://gist.github.com/seantalts/f371ed854662ef185aff33275badc0d4). Please review!

	## Base rate in your city
	First we'll take a back-of-the-envelope stab at modeling the number of people in your city who are out spreading the disease. Take a look at adjusting the parameters according to your best guesses for your city based on recent news and then adjust the standard deviation sliders until the graphs for simulated values look like they reasonably capture your uncertainty about those numbers.
	```{r echo=FALSE}
	inputPanel(
	numericInput("days_infectious", label = "Number of days we think covid is contagious if you are infected, X", 10),
	numericInput("new_cases", label = "Average number of new cases per day for the past X days", 320),
	numericInput("city_pop", label = "Number of people who are living in your city", 7.5e6),
	sliderInput("asymptomatic_rate", label="Percentage of covid cases that are asymptomatic (and thus not tested)", 75,
	min=0, max=98),
	sliderInput("days_stdev", label="Days infected standard deviation (sd)", 1,
	min=0, max=6),
	sliderInput("cases_stdev", label="New cases stdev %", 20,
	min=2, max=100),
	sliderInput("city_pop_stdev", label="City population stdev %", 2,
	min=0, max=20),
	sliderInput("asym_rate_stdev", label="Asymptomatic rate stdev", 15,
	min=0, max=40)
	)

	bound_pct = function(draws) {
	n = length(draws)
	draws = draws[draws > 1 & draws < 99]
	missing = n - length(draws)
	bootstrapped = sample(draws, missing, replace=T)
	c(draws, bootstrapped) / 100
	}

	trunc_at_one = function(a) {
	a[a > 1] = 1
	a
	}

	get_draws = reactive({
	n = 10000
	draws_sdpct = function(mean, sd) {
	abs(rnorm(n, mean, sd / 100 * mean ))
	}
	draws = function(mean, sd) {
	abs(rnorm(n, mean, sd))
	}

	df = data.frame(
	days_infectious = draws(input$days_infectious, input$days_stdev),
	cases = draws_sdpct(input$new_cases, input$cases_stdev),
	city_pop = draws_sdpct(input$city_pop, input$city_pop_stdev),
	asym_rate = bound_pct(draws(input$asymptomatic_rate, input$asym_rate_stdev))
	)
	df$num_walking_infected = df$days_infectious * (df$cases / (1 - df$asym_rate) - df$cases)
	df$chance_random_person_infected = trunc_at_one(df$num_walking_infected / df$city_pop)
	df
	})
	# TODO: Change asymptomatic rate to "ascertainment error?" or break it out

	renderText({
	df = get_draws3()
	paste(sep=" ", "Your average chance of being infected across", nrow(df),
	"simulations given all the data on this page is",
	round(mean(df$chance_you_contract) * 100, 3), "%")
	})

	fluidRow(
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$days_infectious,
	main="Histogram of # days infectious used in these simulations",
	xlab="Days",
	ylab="", yaxt='n')
	})),
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$cases,
	main="Histogram of number of new cases yesterday",
	xlab="Cases",
	ylab="", yaxt='n')
	})))

	fluidRow(
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$city_pop,
	main="Histogram of city population",
	xlab="People",
	ylab="", yaxt='n')
	})),
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$asym_rate, breaks=100,
	main="Distribution of asymptomatic infection rates",
	xlab="Infection rate (of total infections)",
	ylab="", yaxt='n')
	})))


	fluidRow(
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$num_walking_infected, breaks=100,
	main="Histogram of # asymptomatic infections at large",
	xlab="People infected",
	ylab="", yaxt='n')
	})),
	column(6,
	renderPlot({
	df = get_draws()
	hist(df$chance_random_person_infected, breaks=100,
	main="Histogram of the chance a random person is infected",
	xlab="Chance",
	ylab="", yaxt='n')
	})))
	```

	## Adding an isolation period
	Here we'll try to model your friends interactions over the past week to see how likely each of them will be to be infected when they arrive. This lets us examine the effects of varying levels of isolation for the X days prior. These uncertainty parameters are pure standard deviations in the units of their associated parameters; transmission rate sd is in % chance of transission and friends of friends sd's unit is people.
	```{r echo=FALSE}
	inputPanel(
	sliderInput("transmission_rate", "Chance of transmission from hanging out with an infected person", 65,
	min=0, max=100),
	sliderInput("trans_rate_stdev", label="Chance of transmission uncertainty (sd)", 25,
	min=1, max=50),
	numericInput("transitive_guests", label="Number of people each guest has seen indoors in the past X days", 3),
	sliderInput("friends_of_friends_stdev", label="Friends-of-friends uncertainty (sd)", 3,
	min=1, max=50))

	bootstrap_positive = function(draws) {
	n = length(draws)
	draws = draws[draws >= 0]
	missing = n - length(draws)
	bootstrapped = sample(draws, missing, replace=T)
	c(draws, bootstrapped)
	}

	get_draws2 = reactive({
	df = get_draws()
	n = nrow(df)
	df$trans_rate = bound_pct(rnorm(n, input$transmission_rate, input$trans_rate_stdev))
	df$friends_of_friends = bootstrap_positive(rnorm(n, input$transitive_guests, input$friends_of_friends_stdev))
	df$chance_friend_contracts = 1 - (1 - df$trans_rate * df$chance_random_person_infected)^df$friends_of_friends
	df
	})

	renderText({
	df = get_draws3()
	paste(sep=" ", "Your average chance of being infected across", nrow(df),
	"simulations given all the data on this page is",
	round(mean(df$chance_you_contract) * 100, 3), "%")
	})

	renderPlot({
	df = get_draws2()
	hist(df$friends_of_friends, breaks=-0.5:ceiling(max(df$friends_of_friends + 1)),
	main="Number of friends each guests hangs out with before your party",
	xlab="People",
	ylab="", yaxt='n')
	})

	fluidRow(
	column(6,
	renderPlot({
	df = get_draws2()
	hist(df$trans_rate, breaks=100,
	main="Distribution of transmission rates",
	xlab="Probability",
	ylab="", yaxt='n')
	})),
	column(6,
	renderPlot({
	df = get_draws2()
	hist(df$chance_friend_contracts, breaks=100,
	main="Histogram of the % chance any friend gets covid elsewhere",
	xlab="Probability",
	ylab="", yaxt='n')
	})))
	```

	## The party
	```{r echo=FALSE}
	inputPanel(
	numericInput("num_guests", label="Number of guests in your party", 17),
	sliderInput("stdev3", label="Standard deviation in guests", 5,
	min=2, max=50))

	get_draws3 = reactive({
	df = get_draws2()
	draws = function(mean) {
	abs(rnorm(nrow(df), mean, input$stdev3 / 100 * mean ))
	}
	df$guests = draws(input$num_guests)
	df$num_friends_infected = rbinom(nrow(df), floor(df$guests), df$chance_friend_contracts)
	df$chance_you_contract = 1 - (1 - df$trans_rate)^df$num_friends_infected
	df
	})

	renderText({
	df = get_draws3()
	paste(sep=" ", "Your average chance of being infected across", nrow(df),
	"simulations given all the data on this page is",
	round(mean(df$chance_you_contract) * 100, 3), "%")
	})

	fluidRow(
	column(6,
	renderPlot({
	df = get_draws3()
	hist(df$guests, breaks=-0.5:ceiling(max(df$guests + 1)),
	main="Histogram of the number of guests in our simulations",
	xlab="People",
	ylab="", yaxt='n'
	)})),
	column(6,
	renderPlot({
	df = get_draws3()
	hist(df$num_friends_infected, breaks=-0.5:ceiling(max(df$num_friends_infected) + 1),
	main="Number of friends infected",
	xlab="People",
	ylab="", yaxt='n'
	)
	})
	))
	```

	#### Detailed stat calculation simulation summaries
	```{r echo=F}
	renderTable(t(do.call(cbind, lapply(get_draws3(), summary))), rownames = T)
	```

	## Method
	### Simulations
	I wanted to do simulations that capture some uncertainty rather than using single estimated numbers for everything and plugging those into the formulas, so every number input is used as the mean of a normal distribution with the percent standard deviation input times the mean as the distribution's standard deviation.

	In maths, $ x_{draws} \sim \mathcal{N}(x, \sigma x) $. For the number of days we think someone is infectious, with days_infectious = 10 and stdev = 10%, this looks like e.g. $ days.infectious \sim \mathcal{N}(10, 1) $ and generates this histogram,
	```{r}
	days_infectious = 10
	stdev = 0.1 * days_infectious
	hist(rnorm(10000, days_infectious, stdev))
	```

	### Formulæ
	Here's the code, but let's walk through it. For non-R folks just ignore the `df$` bits.
	```{r eval=F}
	df$num_walking_infected = df$days_infectious * (df$cases / (1 - df$asym_rate) - df$cases)
	df$chance_random_person_infected = trunc_at_one(df$num_walking_infected / df$city_pop)
	df$chance_friend_contracts = 1 - (1 - df$trans_rate * df$chance_random_person_infected)^df$friends_of_friends
	df$num_friends_infected = rbinom(n, floor(df$guests), df$chance_friend_contracts)
	df$chance_you_contract = 1 - (1 - df$trans_rate)^df$num_friends_infected
	```

	#### Number infected walking around your city
	First we estimate the number of people walking around who are infected somewhat conservatively by saying that only those who are asymptomatic are out and about. If we think `cases = total * rate_of_symptoms` then we can find the total infected by dividing cases by rate of having symptoms (which is 1 - asymptomatic rate). So that gives us `total = cases / (1 - asym_rate)`, to find the number of people walking around spreading we can then just subtract the number of cases.

	#### Number of people at your party infected
	A binomial distribution gives us the number of "successes" in N trials with p probability of success. If we think about success as the chance of infection and each guest as an independent trial, we can use R's `rbinom(n, df$guests, df$chance_friend_contracts)` to model the number of guests who will have it. We have to convert the number of guests to an integer and also pass in n because we're doing everything as n simulations, so we'll end up with 10000 simulations each of which says how many people at your party are infected.

	#### Given X infected people, what is your chance of infection?
	The formula for at least one success in n trials with a probability p of success in each trial is `1 - (1 -p)^n`. We can use the binomial distribution again to find the chance of infection given the number of people at your party who are infected and the transmission rate between people. that gives us `1 - (1 - df$trans_rate)^df$num_friends_infected`. We also use this to estimate the chance that each friend has contracted the virus based on their number of friends, the transmission rate, and the chance a random person is infected.

	This model encodes a 2 step process - as if each friend has an "event" with infection risk from everyone they come into contact with in the past X days before they come to your party or new isolation pod. I did this because I wanted to model the effects of isolation prior to the event.