Search code examples
rstatisticsbonferroni

How to apply p.adjust to each row in matrix in R?


I have a 100 x 10 000 matrix with p-values (pval) which corresponds to 10 000 repetitions of 100 hypotheses. Now I want to apply a bonferroni correction to each of the rows in the matrix, with the function p.adjust.

My code runs, but there are no changes in the p-values compared to the original p-value matrix, and the FWER is still at approximately 0.994 compared to the expected 0.05 level.

Earlier i tried to use apply instead, but realised as I wanted each row to be adjusted, sapply should be more appropriate.

Here is my code:

#we simulate n random variables and test the null hypothesis that all
#means are simultaneously 0. (case where all nulls are true.)
n <- 100
pval <- matrix(NA, nrow = 1e4, ncol = n) #making matrix for data
for(k in 1 : 1e4){
  X <- replicate(n=n, expr = rnorm(100)) #replicate makes a matrix with n columns and 100 rows of norm(1,0)
  pval[k, ] <-  apply(X=X, MARGIN = 2, FUN = function(x){ #apply applies a function on our X matrix. MARGIN = 2 means the function is applied on the columns.
    t.test(x = x, mu = 0)$p.value #the function being applied is t.test. (t test is taken on all columns in the X matrix (100 rows))
  }) #this returns a matrix with n rows, 10000 columns where each column represents a p-value.
} #the data is uncorrelated. all zero - hypotheses are true.

#now we apply the Bonferroni correction:
padjBonf <- matrix(NA, nrow = 1e4, ncol = n) #making matrix for adjusted p-vals
for(k in 1 : 1e4){
  padjBonf[k,] <- sapply(X = pval[k,], FUN = function(x){ #sapply applies a function on our X vector. MARGIN = 2 means the function is applied on the columns.
    p.adjust(x, method = "bonferroni") 
  }) 
} 

Solution

  • Using sapply like that you are running p.adjust on single values

    Observe the following:

    p.adjust( 0.05, method="bonferroni" ) # returns 0.05, unchanged!
    

    This is what you are experiencing.

    You likely meant to give p.adjust all the p-value of each experiment, hence the entire row of your p.val, like so:

    padjBonf[k,] = p.adjust( p.val[k,], method="bonferroni" )
    

    This should return all 1's as apropriate.

    Or you could continue to correct each p-value and tell it that n=100 as documented in the manual, but there's no need really, and p.adjust was written with the above usage in mind.