R/mtdt2.R

In gap: Genetic Analysis Package

#' Transmission/disequilibrium test of a multiallelic marker by Bradley-Terry model
#'
#' @param x the data table.
#' @param verbose To print out test statistics if TRUE.
#' @param n.sim Number of simulations.
#' @param ... other options compatible with the BTm function.
#'
#' @details
#' This function calculates transmission-disequilibrium statistics involving
#' multiallelic marker according to Bradley-Terry model.
#'
#' @export
#' @return
#' It returned list contains the following components:
#' - c2b A data frame in four-column format showing transmitted vs nontransmitted counts.
#' - BTm A fitted Bradley-Terry model object.
#' - X2 Allele-wise, genotype-wise and goodness-of-fit Chi-squared statistics.
#' - df Degrees of freedom.
#' - p P value.
#' - pn Monte Carlo p values when n.sim is specified.
#'
#' @references
#' \insertRef{firth05}{gap}
#'
#' Turner H, Firth D (2010) Bradley-Terry models in R: The BradleyTerry2 package.
#' https://cran.r-project.org/web/packages/BradleyTerry2/vignettes/BradleyTerry.pdf.
#'
#' @seealso [`mtdt`]
#'
#' @examples
#' \dontrun{
#' x <- matrix(c(0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0,
#'               0,0, 1, 3, 0,0, 0, 2, 3, 0, 0, 0,
#'               2,3,26,35, 7,0, 2,10,11, 3, 4, 1,
#'               2,3,22,26, 6,2, 4, 4,10, 2, 2, 0,
#'               0,1, 7,10, 2,0, 0, 2, 2, 1, 1, 0,
#'               0,0, 1, 4, 0,1, 0, 1, 0, 0, 0, 0,
#'               0,2, 5, 4, 1,1, 0, 0, 0, 2, 0, 0,
#'               0,0, 2, 6, 1,0, 2, 0, 2, 0, 0, 0,
#'               0,3, 6,19, 6,0, 0, 2, 5, 3, 0, 0,
#'               0,0, 3, 1, 1,0, 0, 0, 1, 0, 0, 0,
#'               0,0, 0, 2, 0,0, 0, 0, 0, 0, 0, 0,
#'               0,0, 1, 0, 0,0, 0, 0, 0, 0, 0, 0),nrow=12)
#'
#' xx <- mtdt2(x,refcat="12")
#' }
#'
#' @author Jing Hua Zhao
#' keywords models
#' keywords htest

mtdt2 <- function(x,verbose=TRUE,n.sim=NULL,...)
{
  for(p in c("BradleyTerry2")) {
     if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
        if (!requireNamespace(p, quietly = TRUE))
        warning(paste("mtdt2 needs package `", p, "' to be fully functional; please install", sep=""))
     }
  }
  dims <- dim(x)[1]
  colnames(x) <- paste(1:dims,sep="")
  rownames(x) <- paste(1:dims,sep="")
  c2b <- as.data.frame(BradleyTerry2::countsToBinomial(x))
  names(c2b) <- c("allele1","allele2","transmitted","nontransmitted")
  allele1 <- with(c2b, allele1)
  allele2 <- with(c2b, allele2)
  transmitted <- with(c2b, transmitted)
  nontransmitted <- with(c2b, nontransmitted)
  btx <- BradleyTerry2::BTm(cbind(transmitted,nontransmitted), allele1, allele2, ~ allele, id="allele", data=c2b, ...)
  t1 <- btx$null.deviance
  t2 <- btx$deviance
  t3 <- t1-t2
  f1 <- btx$df.residual
  f2 <- btx$df.null
  f3 <- f2-f1
  ctest <- c(t3,t1,t2)
  df <- c(f3,f2,f1)
  pexp <- pchisq(ctest,df,lower.tail=FALSE,log.p=TRUE)/log(10)
  p <- 10^pexp
  transmissions <- with(c2b,transmitted+nontransmitted)
  if(!missing(n.sim))
  {
    transmissionlen <- length(transmissions)
    transmittedn <- with(c2b,transmitted)
    nontransmittedn <- with(c2b,nontransmitted)
    pn <- rep(0,3)
    for(i in 1:n.sim)
    {
       for(j in 1:transmissionlen) transmittedn[j] <- rbinom(1,transmissions[j],0.5)
       nontransmittedn <- transmissions-transmittedn
       c2bn <- data.frame(c2b,transmittedn,nontransmittedn)
       btn <- BradleyTerry2::BTm(cbind(transmittedn,nontransmittedn), allele1, allele2, ~ allele, id="allele", data=c2bn, ...)
       t1 <- btn$null.deviance
       t2 <- btn$deviance
       t3 <- t1-t2
       ctestn <- c(t3,t1,t2)
       if(ctestn[1]>=ctest[1]) pn[1] <- pn[1]+1
       if(ctestn[2]>=ctest[2]) pn[2] <- pn[2]+1
       if(ctestn[3]>=ctest[3]) pn[3] <- pn[3]+1
    }
    pn <- (pn+1)/(n.sim+1)
  }
  if(verbose)
  {
    print(cbind(c2b,transmissions,fitted=btx$fitted.values*transmissions))
    print(summary(btx,corr=TRUE))
    cat("Chi-square for allele-wise TDT =",ctest[1],"df =",df[1],"p =",p[1],"\n")
    cat("Chi-square for genotype-wise TDT =",ctest[2],"df =",df[2],"p =",p[2],"\n")
    cat("Chi-square for goodness-of-fit of allele-wise model =",ctest[3],"df =",df[3],"p =",p[3],"\n")
    if(!missing(n.sim)) cat("The corresponding Monte Carlo p values are",pn,"\n")
  }
  if(missing(n.sim)) invisible(list(c2b=c2b,BTm=btx,X2=ctest,df=df,p=p))
  else invisible(list(c2b=c2b,BTM=btx,X2=ctest,df=df,p=p,pn=pn))
}

# created on 20-4-2010
# last updated on 22-4-2010
# As BTm(data=) is not working, we need get around with attachment or global assignment.