On R
, I used the boostrap method to get a correlation coefficient estimation and the confidence intervals.
To get the p-value, I thought, I can calculate the proportion of the confidence intervals which do not contain zero. But this is not the solution.
How can I get the p-value in this case ?
I am using cor.test
to get the coefficient estimation. cor.test
may also gives me the p-value from every test. But how can I get the bootstrapped p-value ?
Thank you very much !
Below an example :
n=30
data = matrix (data = c (rnorm (n), rnorm (n),rnorm (n), rpois(n,1),
rbinom(n,1,0.6)), nrow = n, byrow = F)
data= as.data.frame(data)
z1 = replicate( Brep, sample(1:dim(data)[1], dim(data)[1], replace = T))
res = do.call ( rbind, apply(z1, 2, function(x){ res=cor.test(data$V1[x], data$V2[x]) ; return ((list(res$p.value,res$estimate))) }))
coeffcorr = mean(unlist(res[,2]), na.rm = T) #bootstrapped coefficient
confInter1 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[1] #confidence interval 1
confInter2 = quantile(unlist(res[,2]), c(0.025, 0.975), na.rm = T)[2] #confidence interval 2
p.value = mean (unlist(res[,1]), na.rm = T ) # pvalue
The standard way of bootstrapping in R is to use base package boot
. You start by defining the bootstrap function, a function that takes two arguments, the dataset and an index into the dataset. This is function bootCorTest
below. In the functionyou subset the dataset selecting just the rows defined by the index.
The rest is straightforward.
library(boot)
bootCorTest <- function(data, i){
d <- data[i, ]
cor.test(d$x, d$y)$p.value
}
# First dataset in help("cor.test")
x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8)
dat <- data.frame(x, y)
b <- boot(dat, bootCorTest, R = 1000)
b$t0
#[1] 0.10817
mean(b$t)
#[1] 0.134634
boot.ci(b)
For more information on the results of functions boot
and boot.ci
see their respective help pages.
EDIT.
If you want to return several values from the boot statistic function bootCorTest
, you should return a vector. In the following case it returns a named vector with the values required.
Note that I set the RNG seed, to make the results reproducible. I should already have done it above.
set.seed(7612) # Make the results reproducible
bootCorTest2 <- function(data, i){
d <- data[i, ]
res <- cor.test(d$x, d$y)
c(stat = res$statistic, p.value = res$p.value)
}
b2 <- boot(dat, bootCorTest, R = 1000)
b2$t0
# stat.t p.value
#1.841083 0.108173
colMeans(b2$t)
#[1] 2.869479 0.133857