R/ELR.R

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

#' Expected LR for pairwise inbred relationship
#'
#' 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)
#'
#' @param p Double vector, allele frequencies
#' @param deltaHP Jacquard cofficients under HP, row vector, length 9
#' @param deltaHT Jacquard cofficients under T, row vector, length 9
#'
#' @details The approach of HKB et al. (2019) is implemented
#' @return The expected value of LR
#' @export
#' @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)
#'

ELRH <- function(p, deltaHP, deltaHT = deltaHP){
  # Eq (1.3) in note:
  c("E(LR(HP; given T))" = as.double(deltaHP %*% LRH(p) %*% t(deltaHT)))
}