compMean: Permutation Test for the Mean

In bozenne/butils: Various useful functions

Permutation Test for the Mean

Description

Permutation test for comparing mean between groups with efficient adjustment for multiple comparisons

Usage

compMean(
  Y,
  group,
  time = NULL,
  permutation.type = "group",
  n.perm = 1000,
  cl = NULL,
  trace = TRUE
)

Arguments

Y	[matrix] Matrix where each column correspond to a difference outcome.
group	[character vector] Group variable.
time	[character vector] Optional time variable. Can only take two different values. The mean difference and the mean change over time between the two groups will then be tested.
permutation.type	[character] Should the group labels ("group") or the outcome labels ("Y") be permuted. In the latter case, exactly two outcomes should be specified and the function will test whether the group differences are the same for both outcomes.
n.perm	[integer, >0] number of permutations.
cl	[integer, >0] Optional cluster or number of cores used to run the permutations in parallel.
trace	[logical] Should a progress bar be displayed to follow the permutations.

Details

Several adjustment for multiple comparisons are available via the print function

• Holm: only accounts for the ordering of the tests to gain efficiency compared to a Bonferroni adjustment.

• max: only accounts for the correlation between the tests to gain efficiency compared to a Bonferroni adjustment.

• Step-down max: accounts both for the correlation between the tests and their ordering to gain efficiency compared to a Bonferroni adjustment.

References

Westfall, P. H. and Troendle, J. F. Multiple Testing with Minimal Assumptions, Biometrical Journal 50 (2008) 5, 745–755 DOI: 10.1002/bimj.200710456b. Dudoit S., Shaffer J. P. and Boldrick J. C. Multiple Hypothesis Testing in Microarray Experiments. Statistical Science 2003, Vol. 18, No. 1, 71–103.

Examples

library(data.table)
set.seed(11)
n <- 100

#### 1- H0 under the global null ####
dt <- rbind(data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "W1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T", time = "W1"))

dt.mean <- dt[,.(meanY1 = mean(Y1)),by=c("group","time")]
dt.mean[group=="T", diff(meanY1)] - dt.mean[group=="C", diff(meanY1)]

res1 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                permutation.type = "group")
print(res1)
print(res1, method = "max")
print(res1, method = "holm")
summary(lm(Y1~time*group, data = dt))$coef
summary(lm(Y2~time*group, data = dt))$coef

res11 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                permutation.type = "Y")
print(res11)
summary(lm(I(Y2-Y1)~group*time, data = dt))$coef

#### 2- H1 under the alternative for the group effect ####
dt <- rbind(data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "W1"),
            data.table(Y1 = rnorm(n)+0.5, Y2 = rnorm(n)+0.5, group = "T", time = "W1"))

dt.mean <- dt[,.(meanY1 = mean(Y1)),by=c("group","time")]
dt.mean[group=="T", diff(meanY1)] - dt.mean[group=="C", diff(meanY1)]

set.seed(10)
res2 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                 permutation.type = "group")
print(res2, method = "max")
summary(lm(Y1~time*group, data = dt))$coef
summary(lm(Y2~time*group, data = dt))$coef

res22 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                permutation.type = "Y")
print(res22, method = "max")
summary(lm(I(Y2-Y1)~group*time, data = dt))$coef

#### 3- H1 under the alternative for the marker effect ####
dt <- rbind(data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "W1"),
            data.table(Y1 = rnorm(n) + 0.5, Y2 = rnorm(n) + 1, group = "T", time = "W1"))

res3 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                permutation.type = "group")
print(res3, method = "max")
summary(lm(Y1~time*group, data = dt))$coef
summary(lm(Y2~time*group, data = dt))$coef

res33 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group, time = dt$time,
                 permutation.type = "Y", n.perm = 1e4, cl = 3)
print(res33)
summary(lm(I(Y2-Y1)~group*time, data = dt))

#### 4- Without time ####
n <- 500
set.seed(10)
dt <- rbind(data.table(Y1 = rnorm(n), Y2 = rnorm(n), Y3 = rnorm(n), group = "C"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), Y3 = rnorm(n), group = "T"))
res4 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2, Y3 = dt$Y3), group = dt$group,
                permutation.type = "group", n.perm = 1e4, cl = 3, trace = TRUE)
print(res4, method = "max")
print(res4, method = "max-step-down")

res44 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), group = dt$group,
                permutation.type = "Y", n.perm = 1e3, trace = TRUE)
res44

## comparison to multtest (BiocManager::install("multtest"))
library(multtest)
mt.maxT(X = t(cbind(Y1 = dt$Y1, Y2 = dt$Y2, Y3 = dt$Y3)), classlabel = as.numeric(as.factor(dt$group))-1)

## comparison to multcomp
library(multcomp)
dt$group <- as.factor(dt$group)
lll <- mmm(Y1 = lm(Y1 ~ group, data = dt),
          Y2 = lm(Y2 ~ group, data = dt),
          Y3 = lm(Y3 ~ group, data = dt))
summary(glht(lll, linfct = mlf(mcp(group = "Dunnett"))))

#### 5- With more than two groups ####
n <- 500
set.seed(10)
dt <- rbind(data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "D1"),
            data.table(Y1 = rnorm(n)+0.1, Y2 = rnorm(n), group = "T1", time = "D1"),
            data.table(Y1 = rnorm(n)+0.2, Y2 = rnorm(n), group = "T2", time = "D1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "C", time = "W1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T1", time = "W1"),
            data.table(Y1 = rnorm(n), Y2 = rnorm(n), group = "T2", time = "W1"))
dtD1 <- dt[time=="D1"]

res5 <- compMean(Y = cbind(Y1 = dtD1$Y1, Y2 = dtD1$Y2), group = dtD1$group,
                permutation.type = "group", n.perm = 1e3, trace = TRUE)
print(res5)

summary(lm(Y1~group, data = dtD1))
summary(lm(Y2~group, data = dtD1))

res55 <- compMean(Y = cbind(Y1 = dt$Y1, Y2 = dt$Y2), time = dt$time, group = dt$group,
                  permutation.type = "group", n.perm = 1e3, trace = TRUE)
print(res55)

summary(lm(Y1~group*time, data = dt))
summary(lm(Y2~group*time, data = dt))

bozenne/butils documentation built on Oct. 14, 2023, 6:19 a.m.