Nothing
#' Full-likelihood Bayes factor
#'
#' Computes the Bayes factor for co-segregation, as described by Thompson et al.
#' (2003).
#'
#' @param x A [pedtools::ped()] object.
#' @param carriers A character vector (or coercible to such), containing the ID
#' labels of pedigree members known to carry the variant in question.
#' @param noncarriers A character vector (or coercible to such), containing the
#' ID labels of pedigree members known *not* to carry the variant in question.
#' @param freq A single number strictly between 0 and 1: the population
#' frequency of the observed allele.
#' @param affected The affected pedigree members.
#' @param unknown Pedigree members with unknown affection status.
#' @param proband The ID label of the proband. This person must also be in both
#' `carriers` and `affected`.
#' @param penetrances Either a numeric vector of length 3, corresponding to
#' `(f0, f1, f2)` or a matrix or data frame with 3 columns. Each row contains
#' the penetrance values of a liability class.
#' @param liability A vector of length `pedsize(x)`, containing for each
#' pedigree member the row number of `penetrances` which should be used for
#' that individual. (If `penetrances` is just a vector, it will be used for
#' all classes.) If `liability` is NULL (the default), it is set to `1` for
#' all individuals.
#' @param details A logical, indicating if detailed output should be returned
#' (for debugging purposes).
#' @param plot A logical.
#' @param ... Optional plot parameters passed on to [pedtools::plot.ped()].
#'
#' @references Thompson D, Easton DF, Goldgar DE. *A full-likelihood method for
#' the evaluation of causality of sequence variants from family data.* Am J
#' Hum Genet, 2003. \doi{10.1086/378100}.
#' @return A positive number. If `details = TRUE`, a list of intermediate
#' results is returned.
#'
#' @examples
#'
#' x = nuclearPed(2)
#'
#' FLB(x, carriers = 3:4, aff = 3:4, unknown = 1:2,
#' freq = 0.0001, penetrances = c(0, 1, 1), proband = 3)
#'
#' @export
FLB = function(x, carriers, noncarriers = NULL, freq,
affected, unknown = NULL, proband,
penetrances, liability = NULL,
details = FALSE, plot = FALSE, ...) {
### Note to self: Don't mess with the order of input checks.
if(!is.ped(x))
stop2("The first argument must be a `ped` object")
labs = labels(x)
if(length(proband) == 0 || proband == "")
stop2("A proband must be specified")
if(!proband %in% labs)
stop2("Unknown proband: ", proband)
if(!proband %in% affected)
stop2("The proband must be affected")
if(!proband %in% carriers)
stop2("The proband must be a carrier")
for(ids in list(proband, affected, unknown, carriers, noncarriers)) {
validID = ids %in% labs
if(!all(validID))
stop2("Unknown ID label: ", ids[!validID])
}
if(length(err1 <- intersect(affected, unknown)))
stop2("Individual specified as both affected and unknown: ", err1)
if(length(err2 <- intersect(carriers, noncarriers)))
stop2("Individual specified as both a carrier and a non-carrier: ", err2)
if(is.null(freq) || is.na(freq))
stop2("An allele frequency must be specified")
if(!is.numeric(freq) || length(freq) != 1 || freq <= 0 || freq >= 1)
stop2("The allele frequency must be a single number strictly between 0 and 1")
if(plot) {
plotSegregation(x, affected, unknown, proband, carriers, noncarriers, ...)
}
# Affection status vector, sorted along labs
aff = logical(pedsize(x))
aff[internalID(x, affected)] = TRUE
aff[internalID(x, unknown)] = NA
# Empty marker (= disease locus)
dis = m = mProband = marker(x, afreq = c(a = 1 - freq, b = freq))
# Full marker
genotype(m, carriers) = c("a", "b")
genotype(m, noncarriers) = c("a", "a")
# Proband marker
genotype(mProband, proband) = c("a", "b")
# Penetrances and liability classes
penetMat = fixPenetrances(penetrances)
liability = liability %||% rep_len(1, pedsize(x))
if(!all(liability %in% 1:nrow(penetMat)))
stop2("Illegal liability class: ", setdiff(liability, 1:nrow(penetrances)))
# Setup for likelihood under causal hypothesis
peelOrder = peelingOrder(x)
peeler = pedprobr:::.peel_M_AUT
peelingProcess = pedprobr:::peelingProcess
startCausal = function(x, locus)
startdata_causative(x, locus, aff = aff, penetMat = penetMat, liability = liability)
likelihoodCausal = function(x, locus) {
peelingProcess(x, locus, startdata = startCausal,
peeler = peeler, peelOrder = peelOrder)
}
# Main Bayes factor
numer1 = likelihoodCausal(x, m)
likDis = likelihoodCausal(x, dis)
likM = likelihood(x, m)
denom1 = likDis * likM
BF1 = numer1/denom1
# Correction factor
likMproband = likelihood(x, mProband)
numer2 = likDis * likMproband
denom2 = likelihoodCausal(x, mProband)
BF2 = numer2/denom2
FLB = BF1 * BF2
if(details)
return(list(c(FLB = FLB, BF1 = BF1, BF2 = BF2),
c(numer1 = numer1, denom1 = denom1, numer2 = numer2, denom2 = denom2),
c(likM = likM, likMproband = likMproband, likDis = likDis)))
FLB
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.