R/likJStates.R

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

#' Likelihood for for all 9 Jacquard states
#'
#' @param a Vector of positive integers (allele 1 of individual 1)
#' @param b Vector of positive integers (allele 2 of individual 1)
#' @param cc Vector of positive integers (allele 1 of individual 2)
#' @param d Vector of positive integers (allele 2 of individual 2)
#' @param pa Double vector of allele frequencies
#' @param pb Double vector of allele frequencies
#' @param pc Double vector of allele frequencies
#' @param pd Double vector of allele frequencies
#' @export
#' @return A matrix with 9 rows and columns equal to the number of markers.
#' Element (i,j) is the likelihood conditioned on Jacquard state i for marker j.
#' @examples
#' p = c(0.1, 0.2, 0.3, 0.4)
#' a  = c(1,1)
#' b  = c(2,2)
#' cc = c(3,1)
#' d  = c(4,2)
#' pa = p[a]
#' pb = p[b]
#' pc = p[cc]
#' pd = p[d]
#' Delta = c(0, 0, 0, 0, 0, 0, 1,0, 0) #unrelated
#' l2 = likJStates(a,b,cc,d, pa,pb,pc,pd)

likJStates = function(a,b,cc,d, pa, pb, pc, pd){
  I12 = a==b
  I13 = a==cc
  I14 = a==d
  I23 = b==cc
  I24 = b==d
  I34 = cc==d
  nMarkers = length(a)

  # create objects of correct size
  PG = matrix(ncol=9, nrow=nMarkers)
  case = character(nMarkers)

  index1 = I12 & I34 & I13 # case 1
  index2  = I12 & I34 & !I13 # case 2
  index3a = I12 & I13 & !I14 # case 3a
  index3b = I12 & !I13 & I14 # case 3b
  index4 = I12 & !I13 & !I14 & !I34 # case = 4
  index5a = I34 & I13 & !I12 # case 5a
  index5b = I34 & I23 & !I12  # case 5b
  index6 = I34 & !I12 & !I13 & !I23 # case 6
  index7a = I13 & I24 & !I12  # case 7a
  index7b = I14 & I23 & !I12  # case 7b
  index8a = I13 & !I24 & !I12 & !I34 # case 8a
  index8b = I14 & !I23 & !I12 & !I34 # case 8b
  index8c = I23 & !I14 & !I12 & !I34 # case 8c
  index8d = I24 & !I13 & !I12 & !I34 # case 8d
  index9 = !I12 & !I34 & !I13 & !I24 & !I14 & !I23 # case 9

  if(any(index1)) {
    case[index1] = '1'
    p = pa[index1]
    PG[index1,] = cbind(p,p^2,p^2,p^3,p^2,p^3,p^2, p^3,p^4)
  }

  if(any(index2)) {
    case[index2] = '2'
    p = pa[index2]
    q = pc[index2]
    PG[index2,] = cbind(0,p*q,0,p*q^2,0,p^2*q,0,0,p^2*q^2)
  }

  if(any(index3a)) {
    case[index3a] = '3a'
    p = pa[index3a]
    q = pd[index3a]
    PG[index3a,] = cbind(0,0,p*q,2*p^2*q,0,0,0,p^2*q,2*p^3*q)
  }

  if(any(index3b)) {
    case[index3b] = '3b'
    p = pa[index3b]
    q = pc[index3b]
    PG[index3b,] = cbind(0,0,p*q,2*p^2*q,0,0,0,p^2*q,2*p^3*q)
  }

  if(any(index4)) {
    case[index4] = '4'
    p = pa[index4]
    q = pc[index4]
    r = pd[index4]
    PG[index4,] = cbind(0,0,0,2*p*q*r,0,0,0,0,2*p^2*q*r)
  }

  if(any(index5a)) {
    case[index5a] = '5a'
    p  = pc[index5a]
    q =  pb[index5a]
    PG[index5a,] = cbind(0,0,0,0,p*q,2*p^2*q,0,p^2*q,2*p^3*q)
  }

  if(any(index5b)) {
    case[index5b] = '5b'
    p  = pc[index5b]
    q =  pa[index5b]
    PG[index5b,] = cbind(0,0,0,0,p*q,2*p^2*q,0,p^2*q,2*p^3*q)
  }

  if(any(index6)) {
    case[index6] = '6'
    p = pc[index6]
    q = pa[index6]
    r = pb[index6]
    PG[index6,] = cbind(0,0,0,0,0,2*p*q*r,0,0,2*p^2*q*r)
  }

  if(any(index7a)) {
    case[index7a] = '7a'
    p = pa[index7a]
    q = pb[index7a]
    PG[index7a,] = cbind(0,0,0,0,0,0,2*p*q,p*q*(p+q),4*p^2*q^2)
  }

  if(any(index7b)) {
    case[index7b] = '7b'
    p = pa[index7b]
    q = pb[index7b]
    PG[index7b,] = cbind(0,0,0,0,0,0,2*p*q,p*q*(p+q),4*p^2*q^2)
  }

  if(any(index8a)) {
    case[index8a] = '8a'
    p = pa[index8a]
    q = pb[index8a]
    r = pd[index8a]
    PG[index8a,] = cbind(0,0,0,0,0,0,0,p*q*r,4*p^2*q*r)
  }

  if(any(index8b)) {
    case[index8b] = '8b'
    p = pa[index8b]
    q = pb[index8b]
    r = pc[index8b]
    PG[index8b,] = cbind(0,0,0,0,0,0,0,p*q*r,4*p^2*q*r)
  }

  if(any(index8c)) {
    case[index8c] = '8c'
    p = pa[index8c]
    q = pb[index8c]
    r = pd[index8c]
    PG[index8c,] = cbind(0,0,0,0,0,0,0,p*q*r,4*q^2*p*r)
  }

  if(any(index8d)) {
    case[index8d] = '8d'
    p = pa[index8d]
    q = pb[index8d]
    r = pc[index8d]
    PG[index8d,] = cbind(0,0,0,0,0,0,0,p*q*r,4*q^2*p*r)
  }

  if(any(index9)) {
    case[index9] = '9'
    p = pa[index9]
    q = pb[index9]
    r = pc[index9]
    s = pd[index9]
    PG[index9,] = cbind(0,0,0,0,0,0,0,0,4*p*q*r*s)
  }
  liks = list()
  Delta = rep(0,9)
  for(s in 1:9){
    Delta[s] = 1
    liks[[s]] = as.numeric(PG %*% Delta)
    Delta[s] = 0
  }
  matrix(byrow = TRUE, nrow = 9, unlist(liks))
}