# Implementations of the likelihood ratio test of Hardy-weinberg equilibrium in sibships across strata with genotying error. # Qiong Li, Helene Massam and Xin Gao # July 11, 2013 # The likelihood ratio tests require the matrix of the numbers of sibships, with the strata in rows, and the genotype configurations of sibpairs in columns. # Example # c <- c(27,23,52,36,8,4,52,63,39,45,22,29) #with HWE # c <- c(25,31,21,37,3,4,75,53,32,37,44,38) #without HWE # num <- matrix(c, nrow=2, ncol=6) # num[1,2] is the number of sibships in strata 1 with the genotype (AA,Aa) # u: the probability of being incorrectly coded as a "a" allele when "A" is a true allele. # y: the probability of being incorrectly coded as a "A" allele when "a" is a true allele. # HWEtest.Strata(num,y,u) HWEtest.Strata(num,y,u) HWEtest.Strata <- function(num,y,u) { v <- 2 # the numbers of alleles n <- (v*v+v)/2 # the numbers of genotype config q <- matrix(0, nrow=6, ncol=6) # Penetrance matrix q[1,1] <- (1-y)^2*(1-y)^2 q[1,2] <- 1/4*(1-y)^2*(1-y)^2+2*1/4*(1-y)^2*u*(1-y)+1/4*u*(1-y)*u*(1-y) q[1,3] <- u*(1-y)*u*(1-y) q[1,4] <- 1/16*(1-y)^2*(1-y)^2+2*1/8*(1-y)^2*u*(1-y)+2*1/16*(1-y)^2*u^2+1/4*u*(1-y)*u*(1-y)+2*1/8*u*(1-y)*u^2+1/16*u^2*u^2 q[1,5] <- 1/4*u*(1-y)*u*(1-y)+2*1/4*u*(1-y)*u^2+1/4*u^2*u^2 q[1,6] <- u^2*u^2 q[2,1] <- (1-y)^2*2*y*(1-y) q[2,2] <- 1/4*(1-y)^2*2*y*(1-y)+1/4*(1-y)^2*(y*u+(1-y)*(1-u))+1/4*u*(1-y)*2*y*(1-y)+1/4*u*(1-y)*(y*u+(1-y)*(1-u)) q[2,3] <- u*(1-y)*(y*u+(1-y)*(1-u)) q[2,4] <- 1/16*(1-y)^2*2*y*(1-y)+1/8*(1-y)^2*(y*u+(1-y)*(1-u))+1/8*u*(1-y)*2*y*(1-y)+1/16*(1-y)^2*2*u*(1-u)+ 1/16*u^2*2*y*(1-y)+1/4*u*(1-y)*(y*u+(1-y)*(1-u))+1/8*u*(1-y)*2*u*(1-u)+1/8*u^2*(y*u+(1-y)*(1-u))+1/16*u^2*2*u*(1-u) q[2,5] <- 1/4*u*(1-y)*(y*u+(1-y)*(1-u))+1/4*u*(1-y)*2*u*(1-u)+1/4*u^2*(y*u+(1-y)*(1-u))+1/4*u^2*2*u*(1-u) q[2,6] <- u^2*2*u*(1-u) q[3,1] <- (1-y)^2*y^2 q[3,2] <- 1/4*(1-y)^2*y^2+1/4*(1-y)^2*y*(1-u)+1/4*u*(1-y)*y^2+1/4*u*(1-y)*y*(1-u) q[3,3] <- u*(1-y)*y*(1-u) q[3,4] <- 1/16*(1-y)^2*y^2+1/8*(1-y)^2*y*(1-u)+1/8*u*(1-y)*y^2+1/16*(1-y)^2*(1-u)^2+1/16*u^2*y^2+1/4*u*(1-y)*y*(1-u)+ 1/8*u*(1-y)*(1-u)^2+1/8*u^2*y*(1-u)+1/16*u^2*(1-u)^2 q[3,5] <- 1/4*u*(1-y)*y*(1-u)+1/4*u*(1-y)*(1-u)^2+1/4*u^2*y*(1-u)+1/4*u^2*(1-u)^2 q[3,6] <- u^2*(1-u)^2 q[4,1] <- 2*y*(1-y)*2*y*(1-y) q[4,2] <- 1/4*2*y*(1-y)*2*y*(1-y)+2*1/4*2*y*(1-y)*(y*u+(1-y)*(1-u))+1/4*(y*u+(1-y)*(1-u))^2 q[4,3] <- (y*u+(1-y)*(1-u))^2 q[4,4] <- 1/16*2*y*(1-y)*2*y*(1-y)+2*1/8*2*y*(1-y)*(y*u+(1-y)*(1-u))+2*1/16*2*y*(1-y)*2*u*(1-u)+1/4*(y*u+(1-y)*(1-u))^2+ 2*1/8*(y*u+(1-y)*(1-u))*2*u*(1-u)+1/16*2*u*(1-u)*2*u*(1-u) q[4,5] <- 1/4*(y*u+(1-y)*(1-u))^2+2*1/4*(y*u+(1-y)*(1-u))*2*u*(1-u)+1/4*2*u*(1-u)*2*u*(1-u) q[4,6] <- 2*u*(1-u)*2*u*(1-u) q[5,1] <- 2*y*(1-y)*y^2 q[5,2] <- 1/4*2*y*(1-y)*y^2+1/4*2*y*(1-y)*y*(1-u)+1/4*(y*u+(1-y)*(1-u))*y^2+1/4*(y*u+(1-y)*(1-u))*y*(1-u) q[5,3] <- (y*u+(1-y)*(1-u))*y*(1-u) q[5,4] <- 1/16*2*y*(1-y)*y^2+1/8*2*y*(1-y)*y*(1-u)+1/8*(y*u+(1-y)*(1-u))*y^2+1/16*2*y*(1-y)*(1-u)^2+ 1/16*2*u*(1-u)*y^2+1/4*(y*u+(1-y)*(1-u))*y*(1-u)+1/8*(y*u+(1-y)*(1-u))*(1-u)^2+1/8*y*(1-u)*2*u*(1-u)+ 1/16*2*u*(1-u)*(1-u)^2 q[5,5] <- 1/4*(y*u+(1-y)*(1-u))*y*(1-u)+1/4*(y*u+(1-y)*(1-u))*(1-u)^2+1/4*y*(1-u)*2*u*(1-u)+ 1/4*2*u*(1-u)*(1-u)^2 q[5,6] <- 2*u*(1-u)*(1-u)^2 q[6,1] <- y^2*y^2 q[6,2] <- 1/4*y^2*y^2+2*1/4*y^2*y*(1-u)+1/4*y*(1-u)*y*(1-u) q[6,3] <- y*(1-u)*y*(1-u) q[6,4] <- 1/16*y^2*y^2+2*1/8*y^2*y*(1-u)+2*1/16*y^2*(1-u)^2+1/4*y*(1-u)*y*(1-u)+2*1/8*y*(1-u)*(1-u)^2+1/16*(1-u)^2*(1-u)^2 q[6,5] <- 1/4*y*(1-u)*y*(1-u)+2*1/4*y*(1-u)*(1-u)^2+1/4*(1-u)^2*(1-u)^2 q[6,6] <- (1-u)^2*(1-u)^2 ln1 <- rep(0, 2) ln0 <- rep(0, 2) for (b in 1:2) { nu <- num[b,] m<- 100 # the number of iterations pt1<- matrix(0, m, 3) w<- matrix(0, nrow = 6, ncol = 6) pt1[1, ] <- c(0.1, 0.7, 0.2) # initial value t<- 1 while (t