qkappa {BookEKM}R Documentation

Calculates joint distribution for a pair of individals given IBD probabilities

Description

Based on conditional distribution given IBD from q012, the joint probability distribution for two individuals are given for specified IBD probabilities. Mutations have been removed from some previous examples.

Usage

qkappa(kappa = c(0, 1, 0), q = NULL)

Arguments

kappa

Three reals summing to 1 giving IBD (0,1,2) probabilities

q

The joint probability distribution for two individuals

Value

A matrix giving the joint distribution

Author(s)

Thore Egeland

References

To appear

See Also

FamiliasLocus to specify mutation models.

Examples

require(Familias)
require(paramlink)
#SNP example
p <- c(0.2, 0.8)
L1<- FamiliasLocus(p, maleMutationRate = 0.001, allelenames = 1:length(p))
M <- L1$maleMutationMatrix
qq <- q012(p)
res1 <- qkappa(kappa=c(0,1,0), q=qq)

# Check, Familias
persons <- c("CH", "AF")
sex     <- c("male", "male")
ped1 <- FamiliasPedigree(id = persons, dadid = c("AF", NA),
      momid = c(NA, NA), sex=sex)
pedigrees <- list(isFather= ped1)
CH <- c(2,2)
AF <- c(1,1)
datamatrix <- rbind(AF,CH)
res2 <-FamiliasPosterior(ped1, L1, datamatrix)$likelihoodsPerSystem

#Sibs. One SNP marker with 
qkappa(kappa=c(0.25,0.5,0.25),q012(p=c(0.2,0.8)))

# Example. Table 8.1 in book (Egeland, Kling Mostad, 2016) reproduced and checked.
# The code can easily be modified to deal with general pairwise relationships
# (by changing kappa), markers having more than two alleles 
# (by changing afreq)  and also general mutation models
# (by changing the parameters for FamiliasLocus) 
p <- 0.7; q <- 1-p;m <- 0.05
afreq <- c(p,q)
M <- FamiliasLocus(afreq, allelenames = 1:length(afreq), 
                   MutationRate = m)$maleMutationMatrix
qq <- q012(afreq) #Genotype probs conditioned on IBD
res1 <- qkappa(kappa=c(0,1,0), q=qq) #Unconditional genotype probs
res2 <- qkappa(kappa=c(1,0.0,0), q=qq)
order1 <- c(1,7,4,3,9,6,2,8,5) #Sort as in Table 8.1
LR.mut <- as.vector(res1/res2)[order1]
M <- FamiliasLocus(afreq, allelenames = 1:length(afreq), 
                   MutationRate = 0.00)$maleMutationMatrix
qq <- q012(afreq) #Genotype probs conditioned on IBD
res1.no <- qkappa(kappa=c(0,1,0), q=qq) #Unconditional genotype probs
res2.no <- qkappa(kappa=c(1,0.0,0), q=qq)
LR <- as.vector(res1.no/res2.no)[order1]
AF <- c(rep("AA",3),rep("AB",3),rep("BB",3))
CH <- rep(c("AA","AB","BB"),3)
tab <- data.frame(AF,CH,"H1.no"=as.vector(res1.no)[order1],
                  "H1"=as.vector(res1)[order1],
                  "H2"=as.vector(res2.no)[order1],LR,LR.mut)
# Check against formulae in Table 8.1
LRE <- c(1/p,1/(2*p),0,1/(2*p),1/(4*p*q),1/(2*q),0,1/(2*q),1/q)
all(round(LR,12)==round(LRE,12))
LR.mut.E <- c((1-m)/p,((1-m)*q+m*p)/(2*p*q),m/q,
              1/(2*p),1/(4*p*q),1/(2*q),
              m/p,((1-m)*p+m*q)/(2*p*q),(1-m)/q)
all(round(LR.mut,12)==round(LR.mut.E,12))
#
sum(tab[,3]*tab[,6]) #E(LR|Hp)=1.25 as it should
sum(tab[,5]*tab[,6]) #E(LR|Hd)=1 as it should

[Package BookEKM version 1.0 Index]