# R/functions.R In covidprobability: Estimate the Unit-Wide Probability of COVID-19

#' The remaining individuals who would not be expected to show symptoms yet
#'
#' Every day, a certain number of people are expected to show symptoms, based
#' on the incubation period. This would typically lead to further investigation
#' and ongoing suspicion of an outbreak. This function calculates the proportion
#' of individuals on a given day that would not be expected to have developed
#' symptoms yet. So if no one has developed symptoms, this proportion of people
#' could still have undetected COVID-19.
#'
#' @param t day
#' @param asympt The proportion of positive patients who would be expected not
#'               to ever develop symptoms (true asymptomatic patients).
#' @param mu The mean of a lognormal distribution that approximates the
#'           incubation period for COVID-19. E.g. 1.63 (see reference).
#' @param sigma The standard deviation of a lognormal distribution that
#'              approximates the incubation period for COVID-19. E.g. 0.5
#'              (see reference).
#'
#' @return Proportion who would not be expected to show symptoms yet
#' @references See McAloon et al. <https://bmjopen.bmj.com/content/10/8/e039652/>
#' @export
#'
prop_remaining <- function(t, asympt, mu = 1.63, sigma = 0.5) {
potential_sympt <- 1 - asympt
prop_sympt_by <- stats::plnorm(t, mu, sigma)
prop <- (potential_sympt - (potential_sympt * prop_sympt_by)) + asympt
return(prop)
}

#' Calculate pretest probability change over time
#'
#' Calculates the pretest probability over time, assuming the individual does
#' not develop symptoms, by taking into account the distribution of incubation
#' periods (defined as the time from exposure to symptom onset).
#'
#' @param pre0 Initial pretest probability (on day of exposure)
#' @param asympt The proportion of positive patients who would be expected not
#'               to ever develop symptoms (true asymptomatic patients).
#' @param days Days since exposure for calculation range
#' @param mu The mean of a lognormal distribution that approximates the
#'           incubation period for COVID-19. E.g. 1.63 (see reference).
#' @param sigma The standard deviation of a lognormal distribution that
#'              approximates the incubation period for COVID-19. E.g. 0.5
#'              (see reference).
#'
#' @return pretest probability by day (time series)
#' @references See McAloon et al. <https://bmjopen.bmj.com/content/10/8/e039652/>
#' @export
#'
adjust_pretest <- function(pre0, asympt, days = 14,
mu = 1.63, sigma = 0.5) {
return(pre0 * sapply(1:days, prop_remaining, asympt, mu, sigma))
}

#' Calculate posttest probability from pretest probability and test
#' characteristics
#'
#'
#'
#'
#' @param pretest_prob Pretest probability
#' @param sens Test sensitivity
#' @param spec Test specificity
#'
#' @return posttest probability
#' @export
#'
calc_postest_prob <- function(pretest_prob, sens, spec) {
pretest_odds <- pretest_prob / (1 - pretest_prob)
LR_neg <- (1 - sens) / spec
posttest_odds <- pretest_odds * LR_neg
postest_prob <- posttest_odds / (posttest_odds + 1)
return(postest_prob)
}

#' Calculate post-test probability if testing occurred on each day in a series
#'
#' Given an initial pretest probability, and assuming symptoms never arise, with
#' each passing day the pretest probability will be lower, given the person did
#' not experience symptoms. This returns a vector of posttest probabilities
#' which takes all of the above into account, assuming a negative test on each
#' day. Note this is not a time series, and does not reflect if serial testing
#' were done each day and assumes testing was only done once.
#'
#' @param pre0 The pretest probability on day 0 (at exposure)
#' @param asympt The proportion of infected patients expected to remain
#'               asymptomatic throughout the course of infection
#' @param days Days since exposure for calculation range
#' @param mu The mean of a lognormal distribution that approximates the
#'           incubation period for COVID-19. E.g. 1.63 (see reference).
#' @param sigma The standard deviation of a lognormal distribution that
#'              approximates the incubation period for COVID-19. E.g. 0.5
#'              (see reference).
#' @param sens A vector of sensitivities by day since exposure
#' @param spec The test specificity
#'
#' @return A vector of posttest probabilities
#' @export
#'
posttest_series <- function(pre0, asympt, days = 14, mu = 1.63, sigma = 0.5,
sens, spec) {
pretest <- adjust_pretest(pre0 = pre0, asympt = asympt,
days = days, mu = mu, sigma = sigma)

posttest <- data.frame(x=1:days, y=rep(NA, days))

if (days > length(sens)) { stop("Insufficient sensitivity data for
specified days") }

for (t in 1:days) {
posttest\$y[t] <- calc_postest_prob(pretest[t], sens[t], spec)
}
return(posttest)
}

#' Find the probability of any (at least one) event happening
#'
#' For an event that occurs with probability p, this function returns
#' the probability of an occurrence given n repetitions. p is numeric
#' and can be a vector.
#'
#' The probability that any event p occurs with n repetitions is equal to the
#' reciprocal of the probability that p never occurs. The probability that p
#' never occurs with n repetitions is `(1 - p) ^ n`. Thus, the probability
#' that any event p occurs after n repetitions is `1 - ( (1 - p) ^ n )`.
#'
#' @param n The number of times to repeat the event (independent)
#' @param p The individual probability of the event happening
#' @return The probability of an event with the specified probability, after n
#'         repetitions
#' @export
#' @examples
#' probability_any(1, 0.5)
#' probability_any(2, 0.5)
#' probability_any(2, c(0.5, 1/3, 0.25))
probability_any <- function(n, p) {
return(1 - ( (1 - p)^n))
}

#' Calculate a time series of probability for an individual following exposure
#'
#' The probability that an individual has COVID-19 will change over time as new
#' information is gleaned. The initial probability is the pretest probability
#' (pre0) associated with the nature of the interaction/exposure. This
#' probability will decrease with each passing day that the individual does not
#' develop symptoms. When a test is done, the probability is the posttest
#' probability; this reduces the probability based on the test characteristics
#' at the time of testing. Subsequently, the probability will continue to
#' decrease with each passing day that no symptoms develop. This function
#' returns a time series including those 3 phases.
#'
#' @param test_day Day of PCR test (days since exposure)
#' @param pre0 Pre-test probability of person on day of exposure
#' @param sens A vector of sensitivities by day since exposure
#' @param spec The specificity of the PCR test
#' @param asympt The proportion of infected patients expected to remain
#'               asymptomatic throughout the course of infection
#' @param days Days since exposure for calculation range
#' @param mu The mean of a lognormal distribution that approximates the
#'           incubation period for COVID-19. E.g. 1.63 (see reference).
#' @param sigma The standard deviation of a lognormal distribution that
#'              approximates the incubation period for COVID-19. E.g. 0.5
#'              (see reference).
#' @return A time series of probabilities
#' @export
individual_probability <- function(test_day, pre0, sens, spec,
asympt, days, mu, sigma) {

days, mu, sigma)[1:(test_day-1)]

posttest <- posttest_series(pre0, asympt, days, mu, sigma, sens\$point, spec)
posttest_upper <- posttest_series(pre0, asympt, days,
mu, sigma, sens\$upper, spec)
posttest_lower <- posttest_series(pre0, asympt, days,
mu, sigma, sens\$lower, spec)

posttest_day <- posttest\$y[test_day]
posttest_day_upper <- posttest_upper\$y[test_day]
posttest_day_lower <- posttest_lower\$y[test_day]

after_series <- rep(posttest_day, utils::tail(1:days, 1)-test_day)
after_series_upper <- rep(posttest_day_upper, utils::tail(1:days, 1)-test_day)
after_series_lower <- rep(posttest_day_lower, utils::tail(1:days, 1)-test_day)

additional_symptomatic <- c(sapply(1:days, prop_remaining, asympt, mu, sigma),
NA) - c(NA, sapply(1:days, prop_remaining, asympt,
mu, sigma))

after_series <- after_series -
after_series_upper <- after_series_upper -
after_series_lower <- after_series_lower -

probability_series <- c(pretest_series, posttest_day, after_series)
probability_series_upper <- c(rep(NA, test_day-1), posttest_day_upper,
after_series_upper)
probability_series_lower <- c(rep(NA, test_day-1), posttest_day_lower,
after_series_lower)

series <- data.frame(point = probability_series,
upper = probability_series_upper,
lower = probability_series_lower)
return(series)
}

#' Calculate a time series of unit-wide probability following exposure
#'
#' To calculate the probability that any asymptomatic person has COVID-19,
#' this function treats each person/exposure as independent events and
#' calculates the probability time series using the individuals time series from
#' `individual_probability()`.
#'
#' @param test_day Day of PCR test (days since exposure)
#' @param pre0 Pre-test probability of person on day of exposure
#' @param sens A vector of sensitivities by day since exposure
#' @param spec The specificity of the PCR test
#' @param asympt The proportion of infected patients expected to remain
#'               asymptomatic throughout the course of infection
#' @param days Days since exposure for calculation range
#' @param mu The mean of a lognormal distribution that approximates the
#'           incubation period for COVID-19. E.g. 1.63 (see reference).
#' @param sigma The standard deviation of a lognormal distribution that
#'              approximates the incubation period for COVID-19. E.g. 0.5
#'              (see reference).
#' @param n Number of exposed individuals
#' @return The probability of an event with the specified probability, after n
#'         repetitions
#' @export
unit_probability <- function(test_day, pre0, sens, spec, asympt, days, mu,
sigma, n) {

series <- individual_probability(test_day = test_day, pre0 = pre0,
sens = sens, spec = spec, asympt = asympt,
days = days, mu = mu, sigma = sigma)

series_n <- data.frame(point = probability_any(n, series\$point),
upper = probability_any(n, series\$upper),
lower = probability_any(n, series\$lower))

return(series_n)
}

## Try the covidprobability package in your browser

Any scripts or data that you put into this service are public.

covidprobability documentation built on Feb. 11, 2021, 5:08 p.m.