moments {BookEKM}R Documentation

Univariate and bivariate expectation and variance of LR

Description

Expectation and variance of LR is calculated exactly. The true pedigree may differ from the one for which LR is calculated. The expected values depend only on allele frequencies and recombination fraction (for two linked markers) whereas the variance also depends on the allele frequencies.

Usage

moments(p1 = c(0.5, 0.5), km = c(0.5, 0.5, 0), p2 = NULL, km.marked = NULL)

Arguments

p1

Vector of allele frequencies for first marker.

km

The IBD probabilities if only one marker. A 3 by 3 matrix if two markers.

p2

Vector of allele frequencies for second marker marker.

km.marked

The IBD probabilities for true pedigree if only one marker. A 3 by 3 matrix if two markers.

Details

The current function replaces several earliers as the exact expected value and variance are calculated up to a pair of linked markers for arbitrary non-inbred pairwise relationships.

Value

E.LR

Expected LR

var.LR

Variance of LR

Author(s)

Thore Egeland <Thore.Egeland@gmail.com>

References

Egeland and Slooten (2016). "The likelihood ratio as a random variable for linked markers in kinship analysis". (Submitted).

Examples

### Comment: For the examples of the paper we check moments (the function)
#   against exact formula given in the paper.

### Example 1 (Example 3.1, i.e., first in paper).
#   Expectation and variance for half sibs; exact and numerically.
k <- c(0.5,0.5,0)
p <- c(0.4,0.6)
L <- length(p)
(mu.Exact <- (L+15)/16)
s <- 0
 for (a in 1:(L-1))
   for (b in (a+1):L)
    s <- s + p[a]/p[b]+p[b]/p[a]
(var.Exact <- (11*L+53)/64+s/128-((L+15)/16)^2) 
moments(p1=p,km=k)

### Example 2 (Example 3.2 in paper).
#   Expectation and variance for avuncular when PO is true
(mu.Exact <- (L+7)/8)
(var.Exact <- 0.25*((9*L+23)/8+s/16)-mu.Exact^2)
moments(p1=p,km=k, km.marked=c(0,1,0))

### Example 3 (Example 3.3 in paper).
#   Expectation and variance for base relationships 
#  (UNR, PO, MZ) for two linked markers
km <- matrix(c(1,0,0,0,0,0,0,0,0),3,3,byrow = TRUE)
unrelated <- moments(p1=p, p2=p, km=km, km.marked=km)

km <- matrix(c(0,0,0,0,1,0,0,0,0),3,3,byrow = TRUE)
PO <-  moments(p1=p, p2=p, km=km, km.marked=km)
L <- L1 <- L2 <- length(p)
mu.Exact <- (L1+3)/4*(L2+3)/4
var.Exact <- ((5*L1+3)/8+s/16)*((5*L2+3)/8+s/16)-mu.Exact^2

km <- matrix(c(0,0,0,0,0,0,0,0,1),3,3,byrow = TRUE)
MZ <-  moments(p1=p, p2=p, km=km, km.marked=km)
mu.Exact <- L1*(L1+1)/2*L2*(L2+1)/2

s1 <- sum(1/p^2)
s2 <- 0
 for (a in 2:L)
   for (b in 1:(a-1))
    s2 <- s2 + 1/(2*p[a]*p[b])
var.Exact <- (s1+s2)^2 - mu.Exact^2

var.Exact <- ((5*L1+3)/8+s/16)*((5*L2+3)/8+s/16)-mu.Exact^2

### Example 4 (Example 3.4 and Figure 2)
#   The below function calculates the expected values
Esymmetric <- function(rho,L1,L2,relationship="Grand"){
term1 <- (3/4)+(1/64)*(L1+3)*(L2+3)
if(relationship=="Grand")
term2 <- (1/64)*(L1-1)*(L2-1)*rho*(rho-2)
else if(relationship=="Half")
term2 <- (1/16)*(L1-1)*(L2-1)*rho*(rho-1)*(rho^2-rho+1)
else if(relationship=="Uncle")
term2 <- (1/256)*(L1-1)*(L2-1)*rho*(4*rho^2-8*rho+5)*(4*rho^3-8*rho^2+5*rho-4)
term1+term2
}
# Illustration
rho <- 0.25
rel <- c("Grand", "Half", "Uncle")
res <- NULL
for (i in 1:3){
 formula <- Esymmetric(rho,2,2, rel[i])
 km=kappaMat(rho,rel[i])
 mom <- moments(p1=p,p2=p,km=km,km.marked=km)$E.LR
 res <- rbind(res,c(rel[i],formula,mom))
}

rho <- 0.25;L1=10;L2=15
p1 <- rep(1,L1)/L1; p2 <- rep(1,L2)/L2;
y1 <- Esymmetric(rho, L1,L2, relationship="Grand")
y1.c <- moments(p1=p1,km=kappaMat(rho,relationship="Grand"), p2=p2)
y1==y1.c$E.LR

y1 <- Esymmetric(rho, L1,L2, relationship="Half")
y1.c <- moments(p1=p1,km=kappaMat(rho,relationship="Half"), p2=p2)
y1==y1.c$E.LR

y1 <- Esymmetric(rho, L1,L2, relationship="Uncle")
y1.c <- moments(p1=p1,km=kappaMat(rho,relationship="Uncle"), p2=p2)
y1==y1.c$E.LR

plot(c(0,0.5),c(0,8.5),type="n", xlab="rho", ylab="E[LR]")
rho <- seq(0,0.5,length=1000)
L1 <- 10; L2 <- 15
y1 <- Esymmetric(rho, L1,L2, relationship="Grand")
lines(rho,y1,lty=1)
y1 <- Esymmetric(rho, L1,L2, relationship="Half")
lines(rho,y1,lty=3)
y1 <- Esymmetric(rho, L1,L2, relationship="Uncle")
lines(rho,y1,lty=2)

L1 <- 15; L2 <- 15
y1 <- Esymmetric(rho, L1,L2, relationship="Grand")
lines(rho,y1,lty=1)
y1 <- Esymmetric(rho, L1,L2, relationship="Half")
lines(rho,y1,lty=3)
y1 <- Esymmetric(rho, L1,L2, relationship="Uncle")
lines(rho,y1,lty=2)

L1 <- 10; L2 <- 30
y1 <- Esymmetric(rho, L1,L2, relationship="Grand")
lines(rho,y1,lty=1)
y1 <- Esymmetric(rho, L1,L2, relationship="Half")
lines(rho,y1,lty=3)
y1 <- Esymmetric(rho, L1,L2, relationship="Uncle")
lines(rho,y1,lty=2)

legend("bottomleft",lty=1:3,c("Grand","Uncle","Half"))




### Example 5 (Example 3.5 in paper)  
# Two markers. Expectation and variance. Grandparent, true PO
rho <- 0.25
km <- kappaMat(rho, relationship="Grand")
km.marked <- matrix(c(0,0,0,0,1,0,0,0,0),3,3)
p1 <- c(0.98,0.01,0.01); p2 <- rep(1,10)/10
#p1 <- c(0.5,0.5); p2 <- c(0.5,0.5)
foo1 <- moments(p1=p1, p2=p2,km=km, km.marked=km.marked)

#Formula implemented below
gr <- function(p1,p2, rho, E){
 L1 <- length(p1)
 L2 <- length(p2)
 r <- rho
 sum1 <- sum2 <- 0
 if(L1>1){
 for (a in 1:(L1-1))
   for (b in (a+1):L1)
    sum1 <- sum1 + p1[a]/p1[b]+p1[b]/p1[a]
	}
 if(L2>1){
 for (a in 1:(L2-1))
   for (b in (a+1):L2)
    sum2 <- sum2 + p2[a]/p2[b]+p2[b]/p2[a]
	}
z2 <- (5*L2+3)/8+sum2/16
z1 <- (5*L1+3)/8+sum1/16
A <- (L1+3)/4
B <- (L2+3)/4
term1 <- (1-r)^2/4
term2 <-(1/4)*(2*r*(1-r)*(A+B)+(2*(1-r)^2+2*r^2)*A*B)
term3 <- (1/4)*(r^2*(z1+z2)+2*r*(1-r)*(z1*B+z2*A) + (1-r)^2*z1*z2)-E^2
list(var.alt=term1+term2+term3)
}
foo2 <- gr(p1,p2, rho, E=foo1$E.LR)$var.alt
abs(foo1$var.LR - foo2) < 1e-12 

### Example 6 (Example 3.7 in paper)
LR2sibs <- function(rho=0.5,L1=2,L2=3){
R <- rho^2+(1-rho)^2
t1<- (1/64)*(2*L1^2+L1*L2+13*L1+2*L2^2+13*L2+33)
t2 <- (1/64)*(L1-1)*(L2-1)*(4*R+2*(L1+L2+2)*R^2+L1*L2*R^4)
t1+t2
}
LR2sibs(L1=50,L2=17,rho=0.25)
km <- kappaMat(0.25,relationship="Sibs")
p1 <- rep(1,50); p1 <- p1/sum(p1)
p2 <- rep(1,17); p2 <- p2/sum(p2)
moments(p1,p2,km=km)

### Example 7 (Table in Egeland and Slooten (2016))
table1 <- function(rho = c(0, 0.25, 0.5),p1 = c(0.5,0.5),p2 = c(1,1,1)/3){
tab <- NULL
rels <- c("UN", "PO", "MZ", "Grand", "Half", "Uncle", "Sibs")
for (re in rels){
line0 <- NULL
for (r in rho){
  k <- kappaMat(rho=r,relationship = re)
  foo <- moments(p1=p1,p2=p2,km=k)
  foo[[2]] <- sqrt(foo[[2]])
  line0 <- c(line0, foo)
  }
  tab <- rbind(tab, line0)
}
tab
}
tab <- table1()
rels <- c("UN", "PO", "MZ", "Grand", "Half", "Uncle", "Sibs")
rownames(tab) <- rels
colnames(tab) <- c("E(0)","sd(0)","E(0.25)","sd(0.25)","E(0.5)","sd(0.5)")

### Example 8 (Table 3 in Egeland and Slooten (2016) for SEE3 and D6S1043)
library(Familias)
data(NorwegianFrequencies)
p1 <- as.double(NorwegianFrequencies[["SE33"]]) #Next from 
p2 <- c(0.0009,0.0054,0.2165,0.2301,0.1060,0.0643,0.0389, 0.0082,0.0516,0.0978,0.0009,0.0879,0.0018,0.0525,0.0018,0.0109,0.0100,0.0018,0.0045,0.0036,0.0009,0.0037)
tab <- table1(p1=p1,p2=p2,rho=c(0,0.044,0.5))
rownames(tab) <- rels
colnames(tab) <- c("E(0)","sd(0)","E(0.044)","sd(0.044)","E(0.5)","sd(0.5)")

### Example 9 (Appendix B Expectation and variance for first cousins; exact and numerically.)
k <- c(0.75,0.25,0)
p <- c(0.4,0.6)
L <- length(p)
(E <- (L+63)/64)
s <- 0
 for (a in 1:(L-1))
   for (b in (a+1):L)
    s <- s + p[a]/p[b]+p[b]/p[a]
(var.Exact <- (23*L+489)/512+s/1024-E^2)
(foo <- moments(p1=p,km=k))



[Package BookEKM version 1.0 Index]