# Biostatistics and Bioinformatics Branch (BBB)

Skip sharing on social media links

### Dominant

{
#2. Dominant ft
ft.D <- function(x,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
{
f1 <- (c1+k1+t1)/x[1]-(c2+k2+t2)/(1-x[1])-2*(N+M+S)*(x[2]-x[1]*x[2]+x[1]-1)/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)
f2 <- (c3+k3+t3+c4+k4+t4)/x[2]-(N+M+S)*(2*x[1]-x[1]^2)/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)
(f <- rbind(f1,f2))
}

#2. Dominant Jacobian
Jac.D <- function(x,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
{
J <- matrix(0,ncol=2,nrow=2)
J[1,1] <- -(c1+k1+t1)/x[1]^2-(c2+k2+t2)/(1-x[1])^2-2*(N+M+S)*((1-x[2])/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)-2*(x[2]-x[1]*x[2]+x[1]-1)^2/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)^2)
J[1,2] <- -2*(N+M+S)*((1-x[1])/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)-(x[2]-x[1]*x[2]+x[1]-1)*(2*x[1]-x[1]^2)/(2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)^2)
J[2,1] <- J[1,2]
J[2,2] <- -(c3+k3+t3+c4+k4+t4)/x[2]^2+(N+M+S)*(2*x[1]-x[1]^2)^2/((2*x[1]*x[2]-x[1]^2*x[2]+x[1]^2-2*x[1]+1)^2)
J
}

#2.Dominant Newton
newton.D <- function(x,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
{
max <- 100
eps <- 1e-10
xx <- x
for (ii in 1:max)
{
JJ <- Jac.D(xx,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
if (kappa(JJ)>1e+10)
{
break
}
xx <- xx-solve(JJ)%*%ft.D(xx,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
}
return(list(JJ,xx))
}

#2.Dominant Solve function
{
p0 <- matrix(0,Nsnp,1)
p1 <- matrix(0,Nsnp,1)
psi1 <- matrix(0,Nsnp,1)
LR.D <- matrix(0,Nsnp,1)
LRT.D <- matrix(0,Nsnp,1)
pvalue.D <- matrix(0,Nsnp,1)

for (i in 1:Nsnp)
{
n <- array(0)
m <- array(0)
s <- array(0)

#1st: Full Triad + Parent Child + Case only
if (model==1)
{
for (j in 1:10)
{
}
for (j in 1:7)
{
}
for (j in 1:3)
{
}
}

#2nd: Full Triad + Parent Child
if (model==2)
{
for (j in 1:10)
{
}
for (j in 1:7)
{
}
for (j in 1:3)
{
s[j] <- 0
}
S <- 0
}

if (model==3)
{
for (j in 1:10)
{
}
for (j in 1:7)
{
m[j] <- 0
}
M <- 0
for (j in 1:3)
{
s[j] <- 0
}
S <- 0
}

c1 <- 4*n[1] + 3*n[2] + 3*n[3] + 2*n[4] + 2*n[5] + 2*n[6] + 2*n[7] + n[8] + n[9]
c2 <- n[2] + n[3] + 2*n[4] + 2*n[5] + 2*n[6] + 2*n[7] + 3*n[8] + 3*n[9] + 4*n[10]
c3 <- n[3] + n[4] + n[6] + n[8]
c4 <- n[1] + n[2] + n[5]
k1 <- 3*m[1] + 2*m[2] + 2*m[3] + m[4] + m[5] + m[6]
k2 <- m[2] + m[3] + m[4] + 2*m[5] + 2*m[6] + 3*m[7]
k3 <- m[2] + m[4] + m[6]
k4 <- m[1] + m[3]
t1 <- 2*s[3] + s[2]
t2 <- s[2] + 2*s[1]
t3 <- s[2]
t4 <- s[3]

#Dominant sol
p0[i] <- (c1+k1+t1)/(c1+c2+k1+k2+t1+t2)
x0 <- matrix(c(p0[i],0.5),nrow=2)
rtn <- newton.D(x0,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
JJ <- do.call(rbind, rtn[1])
sol <- do.call(rbind, rtn[2])
if (norm(ft.D(sol,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)) > 1e-10 | norm(sol) > 1e+10 | sol[1] > 1 | sol[2] < 0)
{
sol <- 0
p0[i] <- (c1+k1+t1)/(c1+c2+k1+k2+t1+t2)
p1[i] <- NA
psi1[i] <- NA
}
else
{
p0[i] <- (c1+k1+t1)/(c1+c2+k1+k2+t1+t2)
p1[i] <- sol[1]
psi1[i] <- sol[2]
}
LR.D[i] <- (p0[i]/p1[i])^(c1+k1+t1)*((1-p0[i])/(1-p1[i]))^(c2+k2+t2)*(1/psi1[i])^(c3+k3+t3+c4+k4+t4)*(2*p1[i]*psi1[i]-p1[i]^2*psi1[i]+p1[i]^2-2*p1[i]+1)^(N+M+S)
LRT.D[i] <- round(-2*log(LR.D[i]), digits=3)
pvalue.D[i] <- round(1 - pchisq(LRT.D[i], df = 1), digits=3)
}

result.D <- cbind(round(p0,digits=3), round(p1,digits=3), round(psi1,digits=3), LRT.D, pvalue.D)
if (model==1) colnames(result.D) <- c('p0', 'p', 'psi1', 'FT+PC+CO LRT', 'p-value')
if (model==2) colnames(result.D) <- c('p0', 'p', 'psi1', 'FT+PC LRT', 'p-value')
if (model==3) colnames(result.D) <- c('p0', 'p', 'psi1', 'FT LRT', 'p-value')
return(result.D)
}

return(cbind(result.D1, result.D2, result.D3))
}
Last Updated Date: 09/03/2013
Last Reviewed Date: 09/03/2013

## Contact Information

Name: Dr Paul Albert
Chief and Senior Investigator
Biostatistics and Bioinformatics Branch
Phone: 301-496-5582
E-mail: albertp@mail.nih.gov

Staff Directory