relationLR: Relationship Likelihood Ratio

Description

Computes likelihood for two pedigrees and their ratio, the likelihood ratio (LR).

Usage

relationLR(
  ped_numerator,
  ped_denominator,
  ids,
  alleles,
  afreq = NULL,
  known_genotypes = list(),
  loop_breakers = NULL,
  Xchrom = FALSE,
  plot = TRUE,
  title1 = "",
  title2 = ""
)

Arguments

ped_numerator

a linkdat object, or a list of several linkdat and/or singleton objects, describing the relationship corresponding to the hypothesis H1 (numerator). If a list, the sets of ID labels must be disjoint, that is, all ID labels must be unique.

ped_denominator

a linkdat object, or a list of several linkdat and/or singleton objects, describing the relationship corresponding to the hypothesis H2 (denominator). ID labels must be consistent with ped_claim.

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 (a, b, c), indicating that individual a has genotype b/c. Missing value is 0.

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.

Value

lik.numerator

likelihood of data given ped_numerator

lik.denominator

likelihood of data given ped_denominator

likelihood ratio lik.numerator/lik.denominator

Details

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.

Examples

Run this code

# NOT RUN {
############################################
# 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)

# }