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")
})
}
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.