ELRH: Expected LR for pairwise inbred relationship

In thoree/inbred: Different functions related to LR, bootstrap etc

Description

The expected value of LR = P(data|HP)/P(data|HD=unrelated) is calculated assuming that the data is generated by T. There are some examples verifying formulae for the inbred (Jacquard) case in HKB et al. (2019)

Usage

ELRH(p, deltaHP, deltaHT = deltaHP)

Arguments

p	Double vector, allele frequencies
deltaHP	Jacquard cofficients under HP, row vector, length 9
deltaHT	Jacquard cofficients under T, row vector, length 9

Details

The approach of HKB et al. (2019) is implemented

Value

The expected value of LR

Examples

require(ribd)
require(pedprobr)

# Example 1. From inbred.pdf
# HP: 1, whose parents are first cousins, i.e., f = 1/16, is father of 3
# HD: 1 and 3 are unrelated
# HP = T

# First expected value from Eq (1.2) in inbred.pdf
p = c(0.2, 0.8)
# p = c(1, 2, 3, 4)/10 # Remove comment to try other p
L = length(p)
f = 1/16
s = sum(1/p)
mu1 = (0.5*s-(L+1)/4)*f^2+(L-1)/2*f+(L+3)/4

# Next using inbreed function, based on matrix formulation
x = nuclearPed(1)
founderInbreeding(x,1) = f
deltaHP = matrix(ribd::condensedIdentity(x, c(1,3)), nrow=1)
mu2 = ELRH(p, deltaHP, deltaHT = deltaHP)
mu1 == mu2 #OK

# Next verifying using oneMarkerDistribution
x = nuclearPed(1)
founderInbreeding(x,1) = f
m = marker(x, alleles = 1:length(p), afreq = p)
pHP = oneMarkerDistribution(x, ids = c(1,3), partialmarker = m, verbose = F)
founderInbreeding(x,1) = 0
pHD = oneMarkerDistribution(x, ids = c(1,2), partialmarker = m, verbose = F)
pHT = pHP
mu3 = sum(pHP/pHD*pHT)
abs(mu3 - mu1) < 10^(-15) #OK

# Example 2 (Section 'The general pairwise inbred case', HKB et al, 2019)
# HP: 1 parent of 3 (no inbreeding)
# HD: 1 and 3 are unrelated
# T: 1, whose parents are related,  is father of 3
# Need to consider Jacquard, kappa does not fully
# specify relationship between 1 and 3.
#
# delta-s below can be calcuated in by ribd::condensedIdentity
f = 1/16
deltaHP = rbind(c(0, 0, 0,0, 0, 0, 0,   1, 0))
deltaHT = rbind(c(0, 0, f,0, 0, 0, 0, 1-f, 0))
mu4 = ELRH(p, deltaHP, deltaHT)

# Check numerically above expectation
x = nuclearPed(1)
m = marker(x, alleles = 1:length(p), afreq = p)
pHP = oneMarkerDistribution(x, ids = c(1,3), partialmarker = m, verbose = F)
pHD = oneMarkerDistribution(x, ids = c(1,2), partialmarker = m, verbose = F)
founderInbreeding(x,1) = f
pHT = oneMarkerDistribution(x, ids = c(1,3), partialmarker = m, verbose = F)
mu5 = sum(pHP/pHD*pHT)
abs(mu4 - mu5) < 10^(-15) #OK

# Check against matrix equation
mu6 = (L+1)/2*deltaHT[3]+(L+3)/4*deltaHT[8]
abs(mu4 - mu6) < 10^(-15) #OK

# Variance numerically
mu2 = sum((pHP/pHD)^2*pHT)
var.1 = mu2 - mu5^2

# Next, check using Eq (1.5) of note
s1 = sum(1/p)
h883 = (3*L+s1)/4
h888 = (5*L+3)/8 + (s1-L)/16
var.2 = deltaHT[3]*h883 + deltaHT[8]*h888 -mu5^2
abs(var.1-var.2) < 10^(-15)

thoree/inbred documentation built on March 28, 2021, 7:42 p.m.