Introduction

Cochran-Armitage Trend Tests (CATT) is well used in genomics to test for association between a single marker and a disease (Zheng et al. 2012, sec. 3.3.1). When the true genetic model is respectively the recessive (REC), additive (ADD), or dominant (DOM) model, the trend test ZCATT(x), where x is equal to 0, 1/2, or 1 respectively, gives powerful tests. Yet, the true model is generally unknown and choosing one specific value of x can lead to poor powers for some alternative models.

Then, the MAX3 statistic defined by \[MAX3 = \max\{|ZCATT(0)|, |ZCATT(1/2)|, |ZCATT(1)|\}\] can be used to have a more robust test (the power of the test remains good whatever is the underlying unkonwn model).

Yet, we could make another robust test based on the idea of the first section.

Simulation of values for these three statistics

I followed the algorithm detailed in (Zheng et al. 2012, sec. 3.9.1) to simulate contingency tables under different parameters as, for example, the genotype relative risk (GRR) \(\lambda_2\), the genetic model, the minor allele frequency (MAF) \(p\), etc.

source("../code/simu_counts.R")
source("../code/ZCATT.R")

Let us plot simulated values of three statistics, ZCATT(0), ZCATT(1/2) and ZCATT(1), by pairs. I will add to these plot two decision boundaries corresponding to the rejection of the null hypothesis for the statistics \(MAX2 = \max\{|S_1|, |S_2|\}\) (the square) and \(DENS2 = \hat{f}_{S_1, S_2}\) (the oval).

source("../code/square.R")
pacman::p_load(ks) 

# Set of parameters
LWD <- 3; PCH <- 19; CEX <- 0.5
XLIM <- YLIM <-  c(-3, 3)
models <- c("NULL", "REC", "ADD", "DOM")
n <- length(models)
NSIM <- 200
lambda2 <- c(1.5, 1/1.5)
p <- 0.3

# Get lots of simulated values for the three statistics under H0
counts <- simu_counts(nsim = 1e5, p = p)             
simus <- sapply(c(0, 0.5, 1), function(x) ZCATT(counts, x = x))
simus.save <- replace(simus, is.na(simus), 0)
colnames(simus.save) <- paste0("ZCATT(", c(0, 0.5, 1), ")")
  
# Plot by pairs 
for (ind in -(1:3)) {
  simus2D <- simus.save[, ind]
  # DENS2
  k <- ks::kde(simus2D)
  plot(k, cont = 95, col = n+1, lwd = LWD,
       xlim = XLIM, ylim = YLIM)
  # MAX2
  q <- quantile(apply(abs(simus2D), 1, max), 0.95) 
  square(q, col = 5, lwd = LWD)
  # H0 + Alternatives
  for (lam2 in lambda2) {
    for (i in 1:n) {
      counts <- simu_counts(nsim = NSIM, model = models[i], 
                            lam2 = lam2, p = p)
      simus <- sapply(c(0, 0.5, 1), function(x) ZCATT(counts, x = x))
      simus <- replace(simus, is.na(simus), 0)[, ind]
      points(simus, col = i, cex = CEX, pch = PCH, 
             lwd = ifelse(i == 1, LWD, 1)) 
    }
  }
  # Legend
  legend(x = "bottomright", legend = models, pch = PCH, col = 1:n)
}

Let us plot these three statistics’ values in 3D:

pacman::p_load(rgl, rglwidget) 

rgl::plot3d(x = simus.save[1:NSIM, ], size = 5, xlab = "ZCATT(0)",
            ylab = "ZCATT(1/2)", zlab = "ZCATT(1)")
for (lam2 in lambda2) {
  for (i in 2:length(models)) {
    counts <- simu_counts(nsim = NSIM, p = p, model = models[i], lam2 = lam2)
    simus <- sapply(c(0, 0.5, 1), function(x) ZCATT(counts, x = x))
    rgl::plot3d(x = simus, col = i, add = TRUE)
  }
}
rglwidget()

We can see the three statistics’ values for the different models (H0: black, REC: red, ADD: green, DOM: blue) are almost on a same plane. So, it would be meaningless to estimate the 3D density.

Comparaison of MAX3 and this new idea

I rather projected these simulated values in 2D through a Principal Component Analysis and then made the new statistic with the 2D density of the first PCs:

  colnames(simus.save) <- NULL
  pca <- prcomp(simus.save, center = FALSE)
  simus2D <- pca$x[, 1:2]
  # DENS2 on PCs
  k <- ks::kde(simus2D)
  plot(k, cont = 95, col = n+1, lwd = LWD)
  # MAX3
  q <- quantile(apply(abs(simus.save), 1, max), 0.95) 
  for (j in 1:3) {
    cart <- pca$rotation[j, 1:2]
    abline(q/cart[2], -cart[1]/cart[2], col = j+1)
    abline(-q/cart[2], -cart[1]/cart[2], col = j+1)
  }
  legend(x = "topright", legend = paste0("ZCATT(", c(0, 0.5, 1), ")"), 
         lwd = 1, col = 2:4)
  # H0 + Alternatives
  for (lam2 in lambda2) {
    for (i in 1:n) {
      counts <- simu_counts(nsim = NSIM, model = models[i], 
                            lam2 = lam2, p = p)
      simus <- sapply(c(0, 0.5, 1), function(x) ZCATT(counts, x = x))
      simus <- replace(simus, is.na(simus), 0)
      simus <- predict(pca, simus)
      points(simus, col = i, cex = CEX, pch = PCH, 
             lwd = ifelse(i == 1, LWD, 1)) 
    }
  }
  # Legend
  legend(x = "bottomright", legend = models, pch = PCH, col = 1:n)

The lines correspond to the decision boundaries of MAX3: for each of the 3 studied statistic, \(H_0\) is rejected for a point which is outside the corresponding two lines. At this point, if I made myself clear enough, you should see that the decision boundary of MAX3 on the previous plot is defined by the convex hull of the 6 lines.

As before, the oval represents the DENS2 statistic (here, on the two first PCs). It does give a robust test, yet slightly less powerful than MAX3. Indeed, to compare these two tests, there is no need for computing their precise powers for each alternative. You just have to compare the number of alternatives (blue, green, and red points) for which \(H_0\) is rejected with one test but not the other, on the plot. At the top and bottom of the plot, we can see a large region where we would reject \(H_0\) with DENS2 and not MAX3, yet there are not many alternatives in there. On the contrary, there are 3 small areas on the left and on the right of the plot where we would reject \(H_0\) with MAX3 and not DENS2, yet there are many alternatives in there.

The problem of DENS2 being slightly less powerful, I think, is that the 2D or 3D distributions of the statistics’ values for the alternatives are highly correlated with the ones of the null hypothesis (they have the same shape and tends to be on the same diagonal). So, even if the “area of H0” is larger with MAX2 or MAX3 than DENS2 (see all 2D plots), the number of alternatives in there is lower, so they give a more powerful test than DENS2 in this case.