Skip Navigation
  Print Page

Biostatistics and Bioinformatics Branch (BBB)

Skip sharing on social media links
Share this:

Generalized Functional Linear Models for Gene-based Case-Control Association Studies

Beta Smooth Only

library(fda)
library(MASS)
library(SKAT)
library(Matrix)

### This Function is modified by Ruzong Fan, March 20, 2013 ###
gflm_beta_smooth_only = function(pheno, mode = "Additive", geno, pos, order, basis, covariate, base = "bspline", interaction = FALSE)
     {
     geno[is.na(geno)]=0
     covariate[is.na(covariate)] = 0
     ### define genotyping matrix for Dom and Rec modes ###
     ### For Dom mode, redefine geno[i,j] = 1 if geno[i,j] = 1 or 2
     ### For Rec mode, redefine geno[i,j] = 1 if geno[i,j] = 2
     geno_X = geno * 0
     for (i in 1:nrow(geno))
          for (j in 1:ncol(geno))
               {
               if (mode == "Dom")
                    {
                    if (geno[i, j] == 1 || geno[i, j] == 2)
                         geno_X[i,j] = 1
                    }
                    else if ( mode == "Rec")
                         if (geno[i, j] == 2)
                              geno_X[i, j] = 1
               }

     if ( mode == "Rec" || mode == "Dom")
          geno = geno_X

     idx = is.na(pheno)
     pheno = pheno[!idx]
     geno = geno[!idx,]
     if (is.vector(covariate)==FALSE )
          {
          covariate = covariate[!idx,]
          } else if (is.vector(covariate))
               {
               covariate[!idx]
               }
     dqr = qr(geno)
     index = dqr$pivot[1:dqr$rank]
     geno = geno[, index]
     pos = pos[index]
     nsample = nrow(geno)
     nsnp = ncol(geno)

     if(max(pos) > 1) {
                         pos = (pos - min(pos)) / (max(pos) - min(pos))
               }

     if (base == "bspline"){
          betabasis = create.bspline.basis(norder = order, nbasis = basis)
          } else if (base == "fspline"){
          betabasis = create.fourier.basis(c(0,1), nbasis = basis)
          }else { }

     B = eval.basis(pos, betabasis)

     UJ = geno %*% B

     ### Make sure UJ has full rank of bbasis or fbasis ###
     UJdqr = qr(UJ)
     UJindex = UJdqr$pivot[1:UJdqr$rank]
     UJ = UJ[, UJindex]
     ###

     #fitNull = glm (pheno ~ covariate, family = "binomial")
     pval = list()

     if (interaction == FALSE)
          {
          fit = glm (pheno ~ covariate + UJ, family = "binomial")
          pval$LRT = anova(fit, test = "LRT")[3,5]
          pval$Chisq = anova(fit, test = "Chisq")[3,5]
          pval$Rao = anova(fit, test = "Rao")[3,6]

          gtfit = gt(pheno ~ covariate, pheno ~ covariate + UJ, model = "logistic")
          pval$gt = p.value(gtfit)
          }else
               {
               Interaction = matrix (0, ncol = ncol(UJ), nrow =nrow(UJ) )
                    for (i in 1:nrow(Interaction) ) { Interaction[i,] = covariate[i,2] * UJ[i,]}

               fitInt = glm (pheno ~ covariate + UJ + Interaction, family = "binomial")
               pval$LRT = anova(fitInt, test = "LRT")[4,5]
               pval$Chisq = anova(fitInt, test = "Chisq")[4,5]
               pval$Rao = anova(fitInt, test = "Rao")[4,6]

               gtfitInt = gt(pheno ~ covariate + UJ, pheno ~ covariate + UJ + Interaction, model = "logistic")
               pval$gt = p.value(gtfitInt)
               }

     pval
     }

Last Updated Date: 04/25/2014
Last Reviewed Date: 04/25/2014

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