Using the FWER controlling z-test procedure of Stange, Bodnar, Dickhaus (2015)"

In MHTcop: Tests Controlling the FDR / FWER under Certain Copula Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(parallel)
library(pbapply)
library("MHTcop")

This vignette demonstrates the usage of the FWER controlling z-test described in J. Stange, T. Bodnar and T. Dickhaus (2015). The paper is available at https://doi.org/10.1007/s10182-014-0241-5. All following references to equations and theorems refer to this paper.

Performing FWER controlling z-tests using fwer.ztest

Let $X_1,\cdots,X_n$ denote an i.i.d. sample with values in $\mathbb{R}^m$. Furthermore let $\mu_j=\mathbb{E}\left[X_{1,j}\right]$ ($j=1,\cdots,m$) be the component-wise expectations. Then the multiple (two-sided) z-test simultaneously tests the hypotheses $H_{0,j}: \mu_j = \mu_j^$ versus the corresponding alternatives $H_{1,j}: \mu_j\not=\mu_j^$.

For values drawn from Model 1a) the test can be performed using fwer.ztest as follows:

set.seed(0)
Sigma_GaussianAR1<-function(rho,m) {
  exponents<-abs(outer(1:m,1:m,"-")) #exponents for AR(1) matrix
  Sigma<-rho^exponents               #computes the AR(1) matrix
  return(Sigma)
}
n <- 100
m <- 10
m0 <- 5
mu <- c(rep(0,m0),runif(m-m0))
sample <- mvtnorm::rmvnorm(n,mu,Sigma_GaussianAR1(0.1,m))
res <- fwer.ztest(sample,rep(0,m))
res$test$Empirical

Thus only the hypothesis $H_{0,7}$ is wrongly not rejected.

Figures 2 and 3: Empirical power and FWER for model 1a) (Gaussian with AR(1) covariance matrix)

The power and FWER of the procedure under model 1a) are estimated using a Monte carlo simulation. To this end samples are repeatedly drawn from a normal distribution

sig <- Sigma_GaussianAR1(rho,m)
samplefun <- function(mu) mvtnorm::rmvnorm(n,mu,sig)

and then the test is applied:

testfun <- function(X) fwer.ztest(X,rep(0,m))

To see the full code for generating the figures execute the following code:

v <- vignette('fwer-ztest',package = 'MHTcop')
file.edit(paste(v$Dir,'doc',v$R,sep='/'))

Figure 2

empiricalPerformance <- function(testfun,samplefun,m,m0) {
  mu <- c(rep(0,m0),runif(m-m0))
  X <- samplefun(mu)
  t <- testfun(X)
  m1 <- m-m0
  empPower <- c()
  if(m1 > 0) empPower <- sapply(t$test,function(t) sum(t[(m0+1):m])/m1)
  c(sapply(t$test,function(t) sum(t[1:m0])),empPower)
}

doPlot <- function(title,params,BB,SS,GG) {
  ue<-max(c(BB,SS,GG))*1.1
  le<-min(c(BB,SS,GG))*0.99

  plot(params,BB,main=title,ylim=c(le,ue),ylab="",xlab=expression(rho),col="green")
  points(params,SS,pch=2,col="blue")
  points(params,GG,pch=3,col="red")
  legend('topleft', legend=c("Bonferroni", "Sidak","Adaptive"),
         col=c("green", "blue", "red"), pch=1:3,lty=NULL)
}

rows <- c("Bonferroni_FWER","Sidak_FWER","Emp_FWER","Bonferroni_Power","Sidak_Power","Emp_Power")

# --------------------------- Figure 2
rhos <- seq(0,0.9,0.1)
K <- 2500
m <- 10
m0 <- 5
n <- 100

data <- matrix(0,6,length(rhos))
rownames(data) <- rows

if(!file.exists('./sim_data/fwer_figure2.rds')) {
  cat("Performing Simulations for Figure 2\n")
  for(i in 1:length(rhos)) {
    cat("\t for rho =",rhos[i],"\n")
    sig <- Sigma_GaussianAR1(rhos[i],m)
    samplefun <- function(mu) mvtnorm::rmvnorm(n,mu,sig)
    testfun <- function(X) fwer.ztest(X,rep(0,m))

    cl <- makeCluster(detectCores(logical=FALSE))
    parallel::clusterSetRNGStream(cl, iseed = 0L)
    parallel::clusterEvalQ(cl,library(MHTcop))
    parallel::clusterExport(cl,c('empiricalPerformance','testfun','samplefun','m','m0','sig','n'))

    res <- pbapply::pbreplicate(K,empiricalPerformance(testfun,samplefun,m,m0),cl=cl)

    stopCluster(cl)

    data[,i] <- rowMeans(res)
  }
  saveRDS(data,'./sim_data/fwer_figure2.rds')
} else {
  data <- readRDS('./sim_data/fwer_figure2.rds')
}

doPlot("Empirical FWER",rhos,data["Bonferroni_FWER",],
       data["Sidak_FWER",],data["Emp_FWER",])

doPlot("Empirical Power",rhos,data["Bonferroni_Power",],
       data["Sidak_Power",],data["Emp_Power",])

Figure 3

# --------------------------- Figure 3
rhos <- c(0.2, 0.5, 0.8)
K <- 2500
m <- 8
m0 <- 8
n <- 100
n.sims <- 300#200

data <- array(0,c(3,length(rhos),n.sims),list(rows[1:3],NULL,NULL))

if(!file.exists('./sim_data/fwer_figure3.rds')){
  cat("Performing Simulations for Figure 3\n")
  for(i in 1:length(rhos)) {
    cl <- makeCluster(detectCores(logical=FALSE))
    parallel::clusterSetRNGStream(cl, iseed = 0L)
    parallel::clusterEvalQ(cl,library(MHTcop))
    for(j in 1:n.sims) {
      cat("\t for rho =",rhos[i],"rep",j,"\n")
      sig <- Sigma_GaussianAR1(rhos[i],m)
      samplefun <- function(mu) mvtnorm::rmvnorm(n,mu,sig)
      testfun <- function(X) fwer.ztest(X,rep(0,m))

      parallel::clusterExport(cl,c('empiricalPerformance','testfun','samplefun','m','m0','sig','n'))

      res <- pbapply::pbreplicate(K,empiricalPerformance(testfun,samplefun,m,m0),cl=cl)

      data[,i,j] <- rowMeans(res)
    }
    stopCluster(cl)
  }
  saveRDS(data,'./sim_data/fwer_figure3.rds')
} else {
  data <- readRDS('./sim_data/fwer_figure3.rds')
}

for(i in 1:3) {
  hist(data["Emp_FWER",i,],xlab=paste("rho =",rhos[i]),main="Empirical Distribution of FWER");
}

Figures 4 and 5: Empirical power and FWER for model 1b) (t-distribution with 9 degrees of freedom and AR(1) covariance matrix)

The power and FWER of the procedure under model 1b) are estimated using a Monte carlo simulation. To this end samples are repeatedly drawn from a t-distribution

Sigma_T_AR1<-function(m=2,rho=0,nu=3) {
  exponents<-abs(outer(1:m,1:m,"-"))  #exponents for AR(1) matrix
  Sigma<-(nu-2)/nu*rho^exponents      #computes the AR(1) matrix
  return(Sigma)
}
samplefun <- function(mu) {
  return(t(mu+t(mvtnorm::rmvt(n,sigma=sig,df=nu))))
}

and then the test is applied:

sig <- Sigma_T_AR1(m,rho,nu)
testfun <- function(X) fwer.ztest(X,rep(0,m))

Figure 4

# --------------------------- Figure 4
rhos <- seq(0,0.9,0.1)
K <- 2500
m <- 10
m0 <- 5
n <- 100

data <- matrix(0,6,length(rhos))
rownames(data) <- rows

if(!file.exists('./sim_data/fwer_figure4.rds')) {
  cat("Performing Simulations for Figure 4\n")
  for(i in 1:length(rhos)) {
    cat("\t for rho =",rhos[i],"\n")

    rho <- rhos[i]
    nu <- 9
    sig<-Sigma_T_AR1(m,rho,nu)

    samplefun <- function(mu) {
      return(t(mu+t(mvtnorm::rmvt(n,sigma=sig,df=nu))))
    }
    testfun <- function(X) fwer.ztest(X,rep(0,m))

    cl <- makeCluster(detectCores(logical=FALSE))
    parallel::clusterSetRNGStream(cl, iseed = 0L)
    parallel::clusterEvalQ(cl,library(MHTcop))
    parallel::clusterExport(cl,c('empiricalPerformance','testfun','samplefun','m','m0','sig','nu','n'))

    res <- pbapply::pbreplicate(K,empiricalPerformance(testfun,samplefun,m,m0),cl=cl)

    stopCluster(cl)

    data[,i] <- rowMeans(res)
  }
  saveRDS(data,'./sim_data/fwer_figure4.rds')
} else {
  data <- readRDS('./sim_data/fwer_figure4.rds')
}

doPlot("Empirical FWER",rhos,data["Bonferroni_FWER",],
       data["Sidak_FWER",],data["Emp_FWER",])

doPlot("Empirical Power",rhos,data["Bonferroni_Power",],
       data["Sidak_Power",],data["Emp_Power",])

Figure 5

# --------------------------- Figure 5
rhos <- c(0.2, 0.5, 0.8)
K <- 2500
m <- 8
m0 <- 8
n <- 100
n.sims <- 300 #200

data <- array(0,c(3,length(rhos),n.sims),list(rows[1:3],NULL,NULL))

if(!file.exists('./sim_data/fwer_figure5.rds')) {
  cat("Performing Simulations for Figure 5\n")
  for(i in 1:length(rhos)) {
    cl <- makeCluster(detectCores(logical=FALSE))
    parallel::clusterSetRNGStream(cl, iseed = 0L)
    parallel::clusterEvalQ(cl,library(MHTcop))

    nu <- 9
    rho <- rhos[i]
    sig <- Sigma_T_AR1(m,rho,nu)

    for(j in 1:n.sims) {
      cat("\t for rho =",rhos[i],"rep",j,"\n")
      samplefun <- function(mu) {
        return(t(mu+t(mvtnorm::rmvt(n,sigma=sig,df=nu))))
      }
      testfun <- function(X) fwer.ztest(X,rep(0,m))

      parallel::clusterExport(cl,c('empiricalPerformance','testfun','samplefun','m','m0','sig','nu','n'))

      res <- pbapply::pbreplicate(K,empiricalPerformance(testfun,samplefun,m,m0),cl=cl)

      data[,i,j] <- rowMeans(res)
      #saveRDS(data,'./sim_data/fwer_figure5.rds')
    }
    stopCluster(cl)
  }
} else {
  data <- readRDS('./sim_data/fwer_figure5.rds')
}

for(i in 1:3) {
  hist(data["Emp_FWER",i,],xlab=paste("rho =",rhos[i]),main="Empirical Distribution of FWER");
}

MHTcop documentation built on May 2, 2019, 7:59 a.m.