p-value output now matched aldex.glm

[kyle-mcgovern]

This code should show that the ALDEx3 and ALDEx2 glm methods produce roughly the same output

TODO: Add tests, probably somehow make below code into the unit test, add more tests for p-values at each step.

Example:

library(ALDEx2)
library(ALDEx3)
set.seed(43254)

## Sim params
N <- 400
D <- 300
DE <- 120

disease <- c(rep(0, N/2), rep(1, N/2))
metadata <- cbind(disease)

## Simulation
sim_A <- matrix(0, nrow=D, ncol=N)
DE_taxa <- sample(1:D, DE)
lfcs <- rep(0, D)
lfcs[DE_taxa] <- c(rnorm(DE-30, 1.5, 1), rnorm(30, -1.5, 1))
for(d in 1:D) {
  for(n in 1:N) {
    sim_A[d,n] <- rnorm(1,5,2)+rnorm(1, metadata[n,1]*lfcs[d], 1)
  }
}
sim_A <- 2^sim_A
sim_Y <- apply(sim_A, 2, function(col) rmultinom(1, 10000, col/sum(col)))

## True LFC in Scale
log_scales <- log2(colSums(sim_A))
mean(log_scales[((N/2)+1):N])-mean(log_scales[1:((N/2))])

## Benchmark
mc.samples <- 2000

gamma_func <- function(X, Y, logWpara) {
  z <- replicate(mc.samples, {
    rnorm(1, 0.5, 0.5)*X[2,]
  })
  return(z)
}

aldex3.res <- aldex.lm(sim_Y, t(cbind(1, metadata)),
                       nsample=mc.samples,
                       GAMMA=gamma_func)

gamma <- replicate(mc.samples, {
  rnorm(1, 0.5, 0.5)*metadata[,1]
})
gamma <- 2^gamma
glm_meta <- cbind(disease)
aldex.obj <- aldex.clr(sim_Y, glm_meta, mc.samples=mc.samples, gamma=gamma)
aldex.res <- aldex.glm(aldex.obj, fdr.method="BH")

all.res <- rbind(apply(aldex3.res$estimate, c(1,2), FUN=`mean`),
                 aldex.res[,":Est"], aldex3.res$p.val.adj,
                 aldex.res[,":pval.padj"], lfcs)
row.names(all.res) <- c("ALDEx3 Intercept", "ALDEx3 Est", "ALDEx2 Est",
                        "ALDEx3 Intercept Adj p-val", "ALDEx3 Adj p-val",
                        "ALDEx2 Adj p-val", "True LFC")
all.res

[jsilve24]

Wonderful, thank you!

Could you please write a very minimial unit test for this? It would be great if it didn't depend on ALDEx2 but was just a calculation of the p-values based on a hand-coded calculation. This would help ensure future changes don't break things without relying on ALDEx2.

[kyle-mcgovern]

I could hard code the p-values or only put aldex2 under Suggests in the Namespace file, which I think(?) means it won't be downloaded by default.

[jsilve24]

hard-code please. Suggests stuff is ok when necessary but can be a pain. Justin

[kyle-mcgovern]

Ok, test added

[jsilve24]

Sorry I was just thinking, we probably want to add BH correction into this as well as it would be impossible to do after the summary p-values are calculated right? i.e., its not something the user could easily do, it needs to be done in the function. If so, could you add?

[kyle-mcgovern]

There is already an option to pass a p.adjust.method="BH" to the function and the function now returns corrected p-values and uncorrected p-values.

[kyle-mcgovern]

Adjusted p-values returned here @jsilve24

[jsilve24]

This looks great! Thank you.