Skip Navigation
  Print Page

Biostatistics and Bioinformatics Branch (BBB)

Skip sharing on social media links
Share this:

Association analysis of complex diseases using triads, parent-child dyads and singleton monads

Unrestricted (General)

FT.UNR = function(triad, dyad, monad)
{
     Nsnp = nrow(triad)
     #5.General Ft
     ft.G <- 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[1]*x[3]+(1-2*x[1])*x[2]+x[1]-1)/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)
          f2 <- (c3+k3+t3)/x[2]-2*(N+M+S)*(x[1]-x[1]^2)/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)
          f3 <- (c4+k4+t4)/x[3]-(N+M+S)*x[1]^2/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)
          (f <-rbind(f1,f2,f3))
     }

     #5.General Jacobian
     Jac.G <- function(x,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
     {
          J <- matrix(0,ncol=3,nrow=3)
          J[1,1] <- -(c1+k1+t1)/x[1]^2-(c2+k2+t2)/(1-x[1])^2-2*(N+M+S)*((x[3]-2*x[2]+1)/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)-2*(x[1]*x[3]+(1-2*x[1])*x[2]+x[1]-1)^2/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)^2)
          J[1,2] <- -2*(N+M+S)*((1-2*x[1])/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)-2*x[1]*(1-x[1])*(x[1]*x[3]+(1-2*x[1])*x[2]+x[1]-1)/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)^2)
          J[1,3] <- -2*(N+M+S)*(x[1]/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)-x[1]^2*(x[1]*x[3]+(1-2*x[1])*x[2]+x[1]-1)/(x[1]^2*x[3]+(2*x[1]-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)/x[2]^2+4*(N+M+S)*x[1]^2*(1-x[1])^2/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)^2
          J[2,3] <- 2*(N+M+S)*x[1]^3*(1-x[1])/(x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)^2
          J[3,1] <- J[1,3]
          J[3,2] <- J[2,3]
          J[3,3] <- -(c4+k4+t4)/x[3]^2+(N+M+S)*x[1]^4/((x[1]^2*x[3]+(2*x[1]-2*x[1]^2)*x[2]+x[1]^2-2*x[1]+1)^2)
          J
     }

     #5.General Newton
     newton.G <- function(x,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
     {
          max <- 1000
          eps <- 1e-10
          xx <- x
          for (ii in 1:max)
          {
               JJ <- Jac.G(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.G(xx,c1,c2,c3,c4,k1,k2,k3,k4,t1,t2,t3,t4,N,M,S)
          }
          return(list(JJ,xx))
     }

     #5.General Solve function
     Solve.G <- function(triad, dyad, monad, model)
     {
          p0 <- matrix(0,Nsnp,1)
          p1 <- matrix(0,Nsnp,1)
          psi1 <- matrix(0,Nsnp,1)
          psi2 <- matrix(0,Nsnp,1)
          LR.G <- matrix(0,Nsnp,1)
          LRT.G <- matrix(0,Nsnp,1)
          pvalue.G <- 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)
                    {
                         n[j] <- triad[i,j]
                    }
                    N <- triad[i,11]
                    for (j in 1:7)
                    {
                         m[j] <- dyad[i,j]
                    }
                    M <- dyad[i,8]
                    for (j in 1:3)
                    {
                         s[j] <- monad[i,j]
                    }
                    S <- monad[i,4]
               }

               #2nd: Full Triad + Parent Child
               if (model==2)
               {
                    for (j in 1:10)
                    {
                         n[j] <- triad[i,j]
                    }
                    N <- triad[i,11]
                    for (j in 1:7)
                    {
                         m[j] <- dyad[i,j]
                    }
                    M <- dyad[i,8]
                    for (j in 1:3)
                    {
                         s[j] <- 0
                    }
                    S <- 0
               }

               #3rd: Full Triad Only
               if (model==3)
               {
                    for (j in 1:10)
                    {
                         n[j] <- triad[i,j]
                    }
                    N <- triad[i,11]
                    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]

               #General sol
               p0[i] <- (c1+k1+t1)/(c1+c2+k1+k2+t1+t2)
               x0 <- matrix(c(p0[i],0.5,0.5),nrow=3)
               rtn <- newton.G(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.G(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
                    psi2[i] <- NA
               }
               else
               {
                    p0[i] <- (c1+k1+t1)/(c1+c2+k1+k2+t1+t2)
                    p1[i] <- sol[1]
                    psi1[i] <- sol[2]
                    psi2[i] <- sol[3]
               }

               LR.G[i] <- (p0[i]/p1[i])^(c1+k1+t1)*((1-p0[i])/(1-p1[i]))^(c2+k2+t2)*(1/psi1[i])^(c3+k3+t3) *(1/psi2[i])^(c4+k4+t4)*(p1[i]^2*psi2[i]+2*p1[i]*(1-p1[i])*psi1[i]+(1-p1[i])^2)^(N+M+S)
               LRT.G[i] <- round(-2*log(LR.G[i]), digits=3)
               pvalue.G[i] <- round(1 - pchisq(LRT.G[i], df = 2), digits=3)
          }

          result.G <- cbind(round(p0,digits=3), round(p1,digits=3), round(psi1,digits=3), round(psi2,digits=3), LRT.G, pvalue.G)
          if (model==1) colnames(result.G) <- c("p0", "p","psi1","psi2", "FT+PC+CO LRT","p-value")
          if (model==2) colnames(result.G) <- c("p0", "p","psi1","psi2", "FT+PC LRT", "p-value")
          if (model==3) colnames(result.G) <- c("p0", "p","psi1","psi2", "FT LRT", "p-value")
          rownames(result.G) <- row.names(triad)
          return(result.G)
     }

     result.G1 = Solve.G(triad, dyad, monad, model=1)
     result.G2 = Solve.G(triad, dyad, monad, model=2)
     result.G3 = Solve.G(triad, dyad, monad, model=3)
     return(cbind(result.G1, result.G2, result.G3))
}
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
Vision National Institutes of Health Home BOND National Institues of Health Home Home Storz Lab: Section on Environmental Gene Regulation Home Machner Lab: Unit on Microbial Pathogenesis Home Division of Intramural Population Health Research Home Bonifacino Lab: Section on Intracellular Protein Trafficking Home Lilly Lab: Section on Gamete Development Home Lippincott-Schwartz Lab: Section on Organelle Biology