moments {BookEKM} | R Documentation |
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.
moments(p1 = c(0.5, 0.5), km = c(0.5, 0.5, 0), p2 = NULL, km.marked = NULL)
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. |
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.
E.LR |
Expected LR |
var.LR |
Variance of LR |
Thore Egeland <Thore.Egeland@gmail.com>
Egeland and Slooten (2016). "The likelihood ratio as a random variable for linked markers in kinship analysis". (Submitted).
### 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))