Relationship Likelihood Ratio
Computes likelihood for two pedigrees and their ratio, the likelihood ratio (LR).
relationLR(ped_numerator, ped_denominator, ids, alleles, afreq = NULL, known_genotypes = list(), loop_breakers = NULL, Xchrom = FALSE, plot = TRUE, title1 = "", title2 = "")
ped_numerator |
a |
ped_denominator |
a |
ids |
genotyped individuals. |
alleles |
a numeric or character vector containing marker alleles names |
afreq |
a numerical vector with allele frequencies. An error is given if they don't sum to 1 (rounded to 3 decimals). |
known_genotypes |
list of triplets |
loop_breakers |
Not yet implemented, only default value NULL currently handled |
Xchrom |
a logical: Is the marker on the X chromosome? |
plot |
either a logical or the character 'plot_only', controlling if a plot should be produced. If 'plot_only', a plot is drawn, but no further computations are done. |
title1 |
a character, title of leftmost plot. |
title2 |
a character, title of rightmost plot. |
This function computes the likelihood of two pedigrees (each corresponding
to a hypothesis describing a family relationship). The likelihood ratio is
also reported. Unlike other implementations we are aware of, partial DNA
profiles are allowed here. For instance, if the genotype of a person is
reported as 1/0 (0 is 'missing') for a triallelic marker with uniform allele
frequencies, the possible ordered genotypes (1,1), (2,1), (1,2), (1,3) and
(3,1) are treated as equally likely. (For general allele frequencies,
genotype probabilities are obtained by assuming Hardy-Weinberg equilibrium.)
A reasonable future extension would be to allow the user to weigh these
genotypes; typically (1,1) may be more likely than the others. If
plot='plot_only', the function returns NULL after producing the plot.
lik.numerator |
likelihood of data given ped_numerator |
lik.denominator |
likelihood of data given ped_denominator |
LR |
likelihood ratio lik.numerator/lik.denominator |
Thore Egeland, Magnus Dehli Vigeland
############################################
# A partial DNA profile is obtained from the person
# denoted 4 in the figure produced below
# There are two possibilities:
# H1: 4 is the missing relative of 3 and 6 (as shown to the left) or
# H2: 4 is unrelated to 3 and 6.
############################################
p = c(0.2, 0.8)
alleles = 1:length(p)
g3 = c(1,1); g4 = c(1,0); g6 = c(2,2)
x1 = nuclearPed(2)
x1 = addOffspring(x1, father = 4, sex = 1, noff = 1)
m = marker(x1, 3, g3, 4, g4, 6, g6, alleles = alleles, afreq = p)
x1 = addMarker(x1, m)
x2 = nuclearPed(2)
x2 = addOffspring(x2, father = 4, sex = 1, noff = 1)
m = marker(x2, 3, g3, 6, g6, alleles = alleles, afreq = p)
x2 = addMarker(x2, m)
missing = singleton(4, sex = 1)
m.miss = marker(missing, g4, alleles = alleles, afreq = p)
missing = addMarker(missing, m.miss)
x2 = relabel(x2, c(1:3, 99, 5:6), 1:6)
known = list(c(3, g3), c(4,g4), c(6, g6))
LR = relationLR(x1, list(x2, missing), ids = c(3,4,6),
alleles = alleles, afreq = p, known = known,
title1 = 'H1: Missing person 4 related',
title2 = 'H2:Missing person 4 unrelated')$LR
# Formula:
p = p[1]
LR.a = (1+p)/(2*p*(2-p))
stopifnot(abs(LR - LR.a) < 1e-10)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.