-
Notifications
You must be signed in to change notification settings - Fork 2
/
gendist_sensspec_cutoff.R
83 lines (71 loc) · 3.51 KB
/
gendist_sensspec_cutoff.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
##' Calculate sensitivity and specificity of a genetic distance cutoff
##'
##' Function to calculate the sensitivity and specificity of a genetic distance cutoff
##' given an underlying mutation rate and mean number of generations between cases
##'
##' @param cutoff the maximum genetic distance at which to consider cases linked
##' @param mut_rate mean number of mutations per generation, assumed to be Poisson distributed
##' @param mean_gens_pdf the density distribution of the mean number of generations between cases;
##' the index of this vector is assumed to be the discrete distance between cases
##' @param max_link_gens the maximum generations of separation for linked pairs
##' @param max_gens the maximum number of generations to consider, if `NULL` (default) value set to the highest
##' number of generations in mean_gens_pdf with a non-zero probability
##' @param max_dist the maximum distance to calculate, if `NULL` (default) value set to max_gens * 99.9th percentile
##' of mut_rate Poisson distribution
##'
##' @return a data frame with the sensitivity and specificity for a particular genetic distance cutoff
##'
##' @author Shirlee Wohl and Justin Lessler
##'
##' @examples
##' # calculate the sensitivity and specificity for a specific genetic distance threshold of 2 mutations
##' gendist_sensspec_cutoff(cutoff=2,
##' mut_rate=1,
##' mean_gens_pdf=c(0.02,0.08,0.15,0.75),
##' max_link_gens=1)
##'
##' # calculate the sensitivity and specificity for a a range of genetic distance thresholds
##' gendist_sensspec_cutoff(cutoff=1:10,
##' mut_rate=1,
##' mean_gens_pdf=c(0.02,0.08,0.15,0.75),
##' max_link_gens=1)
##'
##' @family genetic distance functions
##'
##' @export
##'
gendist_sensspec_cutoff <- function(cutoff, mut_rate, mean_gens_pdf, max_link_gens = 1,
max_gens = NULL, max_dist = NULL) {
if (is.null(max_gens))
max_gens <- which(mean_gens_pdf != 0)[length(which(mean_gens_pdf != 0))]
if (is.null(max_dist))
max_dist <- suppressWarnings(max_gens * stats::qpois(0.999, mut_rate))
# check that we have used a sensible cutoff the mutation rate should be
# high enough such that the cutoff used is less than the max possible
# distance
if (max_dist < max(cutoff + 1)) {
warning("Nonsensical cutoff given distances considered")
}
# get linked and unlinked probability distributions
gendist <- gendist_distribution(mut_rate = mut_rate, mean_gens_pdf = mean_gens_pdf,
max_link_gens = max_link_gens, max_gens = max_gens, max_dist = max_dist)
linked_pdf <- gendist[, "linked_prob"]
unlinked_pdf <- gendist[, "unlinked_prob"]
linked_cdf <- cumsum(linked_pdf)/sum(linked_pdf)
unlinked_cdf <- cumsum(unlinked_pdf)/sum(unlinked_pdf)
# Utility function to allow for multiple sensitivity and specificity
# cutoffs
get_sens_spec <- function(cutoff) {
# remember that cdf[cutoff] represents the probability of cutoff-1 but
# this is what we want because the threshold is <cutoff (not <=)
# decrease specificity by probability mass of unlinked
spec <- 1 - unlinked_cdf[cutoff]
# decrease sensitivity by probability mass of linked
sens <- 1 - (1 - linked_cdf[cutoff])
return(c(sens, spec))
}
rc <- t(sapply(cutoff, get_sens_spec))
rc <- cbind(cutoff, rc)
colnames(rc) <- c("cutoff", "sensitivity", "specificity")
return(rc)
}