These functions compute the expected distributions of distances between pairs of cases given a case reporting probability 'pi'. Analytical results are used for some special cases, including:

dpaircase(x, type = c("temporal", "genetic", "spatial", "empiric"), gamma_shape, gamma_rate = 1, gamma_scale = 1/gamma_rate, poisson_rate, sd_spatial, p, pi, alpha = 0.001) dtemporal(x, shape, rate = 1, scale = 1/rate, pi, alpha = 0.001) dspatial(x, sd, pi, alpha = 0.001) dgenetic(x, gamma_shape, gamma_rate = 1, poisson_rate, pi, gamma_scale = 1/gamma_rate, alpha = 0.001) dempiric(p, pi, alpha = 0.001)

- x
- vector of quantiles.
- type
- type of distance to be considered (one of "temporal","genetic", "spatial" or "empiric").
- gamma_shape, gamma_scale
- shape and scale of the gamma distribution used for the serial interval
- gamma_rate
- an alternative way to specify the scale of the gamma distribution used for the serial interval
- poisson_rate
- rate (i.e. mean) of the poisson distribution used for the per time unit genetic mutation rate
- sd_spatial
- standard deviation of the Normal spatial kernel.
- p
- A
`numeric`

vector providing the probability mass function, or empirical frequencies, of pairwise distances. - pi
- The reporting probability, i.e. the proportion of cases of the outbreak that have been reported.
- alpha
- The probability threshold to be used to determine the maximum value of generations between two successive cases to consider; this value ('max_kappa') will be the smallest k so that p(k > max_kappa) < alpha. Defaults to 0.001.
- shape, scale
- shape and scale of the gamma distribution used for the serial interval
- rate
- an alternative way to specify the scale of the gamma distribution used for the serial interval
- sd
- standard deviation of the Normal spatial kernel.

- temporal distances: the serial interval is assumed to be gamma-distributed; the procedure returns a discretized, weighted convolution of gamma distributions.
- genetic distances: assumed ... TBC

## COMPARE DEMPIRIC AND DTEMPORAL ## Note in this comparison we are not expecting to get exactly the same ## results since dempiric does the convolution between discretised gamma ## distributions whilst dtemporal does the convolution between gamma ## distributions. ## compute empirical distribution correponding to exponential(mean 50) mean_exp <- 50 x <- 0:300 reporting_rate <- 0.5 p <- dgamma(x, shape = mean_exp, rate = 1) ## computes pdf of a gamma distr with shape mean_exp and scale = rate = 1 ## (i.e. an exponential distr with mean mean_exp) ## use this as an empirical distribution to feed into dempiric empiric_exp_distr_r50 <- dpaircase(x, type = "empiric", p = p, pi = reporting_rate) temporal_distr_r50 <- dpaircase(x, type = "temporal", gamma_shape = mean_exp, gamma_rate = 1, pi = reporting_rate) ## compare the two correlation <- cor(empiric_exp_distr_r50, temporal_distr_r50) ## graphical comparison plot(x, empiric_exp_distr_r50, xlab = "Time", ylab = "pdf", main = "Time between linked cases", cex.main = 1, pch = 3)mtext("SI ~ exp(mean=50), pi = 0.5", side = 3)lines(x, temporal_distr_r50, col = "red")legend("topright", c("dempiric","dtemporal"), pch = c(3, -1), lwd = c(-1, 1), col = c("black","red"), bty = "n")## COMPARE DEMPIRIC AND DGENETIC ## compute empirical distribution correponding to ## an Exponential(mean 50)-Poisson(mean 0.6) mixture mean_exp <- 50 mutation_rate <- 0.6 x <- 0:300 reporting_rate <- 0.5 prob <- 1-mutation_rate/(mutation_rate+1) ## pmf of a negative binomial distr with parameters size and prob p <- dnbinom(x, size = mean_exp,prob = prob) ## use this as an empirical distribution to feed into dempiric empiric_exp_distr_r50 <- dpaircase(x, type = "empiric", p = p, pi = reporting_rate) genetic_distr_r50 <- dpaircase(x, type = "genetic", gamma_shape = mean_exp, gamma_rate = 1, poisson_rate = mutation_rate, pi = reporting_rate) ## compare the two correlation <- cor(empiric_exp_distr_r50, genetic_distr_r50) ## graphical comparison plot(x, empiric_exp_distr_r50, xlab = "Number of mutations", ylab = "pmf", main = "Mutations between linked cases", cex.main = 1, pch = 3)mtext("SI ~ exp(mean=50), pi = 0.6", side = 3)lines(x, genetic_distr_r50, col = "red")legend("topright", c("dempiric", "dgenetic"), pch = c(3, -1), lwd = c(-1, 1), col = c("black","red"), bty = "n")