Nothing
#' Transmission/disequilibrium test of a multiallelic marker
#'
#' @param x the data table.
#' @param n.sim the number of simulations.
#'
#' @details
#' This function calculates transmission-disequilibrium statistics involving multiallelic marker.
#' Inside the function are tril and triu used to obtain lower and upper triangular matrices.
#'
#' @export
#' @return
#' It returned list contains the following components:
#' - SE Spielman-Ewens Chi-square from the observed data.
#' - ST Stuart or score Statistic from the observed data.
#' - pSE the simulated p value.
#' - sSE standard error of the simulated p value.
#' - pST the simulated p value.
#' - sST standard error of the simulated p value.
#'
#' @references
#' \insertRef{miller97}{gap}
#'
#' \insertRef{sham95}{gap}
#'
#' \insertRef{spielman96}{gap}
#'
#' \insertRef{zhao99}{gap}
#'
#' @seealso [`bt`]
#'
#' @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)
#'
#' # See note to bt for the score test obtained by SAS
#'
#' mtdt(x)
#' }
#'
#' @author Mike Miller, Jing Hua Zhao
#' @keywords models
mtdt <- function(x,n.sim=0)
{
#
# lower triangular matrix 5/6/2004
#
tril <- function(t)
{
a <- t
a[upper.tri(t)] <- 0
a
}
#
# upper triangular matrix 5/6/2004
#
triu <- function(t)
{
a <- t
a[lower.tri(t)] <- 0
a
}
#
# Spielman & Ewens' statistic (diagnoal elements kept)
# SAS 20/07/1999
#
T <- x
T <- T-diag(diag(T))
m <- dim(T)[1]
#
# Find non-zero off-diagonal elements of table.
#
b <- matrix(rep(0,m*m),nrow=m)
i <- j <- matrix(rep(0,m*(m-1)/2),ncol=1)
for (ii in 1:m)
{
for (jj in 1:(ii-1)) b[ii,jj]=T[ii,jj]+T[jj,ii]
}
l <- 0
for (ii in 1:(m-1))
{
for (jj in (ii + 1):m)
if (b[jj,ii]>0)
{
l <- l + 1
i[l] <- jj
j[l] <- ii
}
}
#
# Count number of different heterozygotes.
#
NH <- l
#
# Count number of each heterozygote.
#
t0 <- T
tc <- apply(t0,2,sum)
tr <- apply(t0,1,sum)
c0 <- (m-1)/m
se0 <- c0*sum((tc-tr)^2/(tc+tr))
if (n.sim>0)
{
C <- rep(0,NH)
for (k in 1:NH) C[k] <- T[i[k],j[k]] + T[j[k],i[k]]
X <- rep(0,n.sim)
for (k in 1:n.sim)
{
for (L in 1:NH)
{
T[i[L],j[L]] <- rbinom(1,C[L],0.5) # sum(runif(C(L))<0.5)
T[j[L],i[L]] <- C[L] - T[i[L],j[L]]
}
tc <- apply(T,2,sum)
tr <- apply(T,1,sum)
X[k] <- c0*sum((tc-tr)^2/(tc+tr))
}
MCp <- 0
for (k in 1:n.sim) if (X[k]>=se0) MCp <- MCp + 1
pSE <- MCp/n.sim
sSE <- sqrt(pSE*(1-pSE)/n.sim)
cat('Spielman-Ewens Chi-square and empirical p (se): ', se0, pSE, sSE, "\n")
} else cat('Spielman-Ewens Chi-square: ', se0, "\n")
#
# Simulate tables and compute TDT chi-square statistics.
#
# Should diag(T) is kept, the statistic will be similar to Spielman-Ewens'
# se.check
T <- x
# This is according to Mike Miller's Matlab program
# Produce inverse (IV) of variance-covariance matrix (V). This is
# constant across repeated samples in the Monte Carlo simulation.
V <- diag(apply(T,1,sum)+apply(T,2,sum))-tril(T)-t(triu(T))-t(tril(T))-triu(T)
IV <- solve(V[1:(m-1),1:(m-1)])
T0 <- T
d0=apply(T0[1:(m-1),1:(m-1)],1,sum)-apply(T0[1:(m-1),1:(m-1)],2,sum)
x0 <- t(d0)%*%IV%*%d0
se.check <- x0
T <- T - diag(diag(T))
V <- diag(apply(T,1,sum)+apply(T,2,sum))-tril(T)-t(triu(T))-t(tril(T))-triu(T)
IV <- solve(V[1:(m-1),1:(m-1)])
d0=apply(T0[1:(m-1),1:(m-1)],1,sum)-apply(T0[1:(m-1),1:(m-1)],2,sum)
st0 <- t(d0)%*%IV%*%d0
if (n.sim>0)
{
C <- rep(0,NH)
for (k in 1:NH) C[k] <- T[i[k],j[k]] + T[j[k],i[k]]
X <- rep(0,n.sim)
for (k in 1:n.sim)
{
for (L in 1:NH)
{
T[i[L],j[L]] <- rbinom(1,C[L],0.5) # sum(runif(C(L))<0.5)
T[j[L],i[L]] <- C[L]-T[i[L],j[L]]
}
d <- apply(T[1:(m-1),1:(m-1)],1,sum)-apply(T[1:(m-1),1:(m-1)],2,sum)
X[k] <- t(d)%*%IV%*%d
}
MCp <- 0
for (k in 1:n.sim) if (X[k]>st0) MCp <- MCp + 1
pST <- MCp/n.sim
sST <- sqrt(pST*(1-pST)/n.sim)
cat('Stuart Chi-square and p (se): ', st0, pST, sST,"\n")
} else {
cat('Stuart Chi-square ',st0, "\n")
cat('Value of Chi-square if diagonal elements are kept: ',se.check,"\n")
}
if (n.sim>0) list(SE=se0,pSE=pSE,sSE=sSE,ST=st0,pST=pST,sST=sST)
else list(SE=se0,ST=st0)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.