P-value for point biserial correlation in R Does anybody know of an R package that produces a p-value for point biserial correlations?
I've tried all of the major packages that I know (with some help from Google) and haven't found any. If a package doesn't come to mind, is there some way that I can intuitively calculate the p-value?
 A: The point-biserial correlation is equivalent to calculating the Pearson correlation between a continuous and a dichotomous variable (the latter needs to be encoded with 0 and 1). Therefore, you can just use the standard cor.test function in R, which will output the correlation, a 95% confidence interval, and an independent t-test with associated p-value:
set.seed(1)
x <- sample.int(100, 50, replace=TRUE)
y <- sample(c(0, 1), 50, replace=TRUE)
cor.test(x, y)

This yields a correlation of $r = 0.202$, which is not significant ($t = 1.429$, $\text{df} = 48$, $p = 0.1595$):
    Pearson's product-moment correlation

data:  x and y
t = 1.429, df = 48, p-value = 0.1595
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.08088534  0.45478598
sample estimates:
      cor 
0.2020105 

As @sal-mangiafico and @igor-p point out, the function biserial.cor from the ltm package produces slightly different results. This is because cor.test uses the population standard deviation, whereas biserial.cor uses the sample standard deviation. Furthermore, the result of biserial.cor has the opposite sign than the result of cor.test. This can be adjusted by specifying the argument level=2 in biserial.cor.
A: In response to @user9413061,
I think I discovered the source of the problem.
In the standard definition of biserial correlation, the population standard deviation is used.
ltm::biserial.cor uses the sample standard deviation.
In the following, a function is defined to calculate the population standard deviation.  The function biserial.cor.new is defined, which is the same as ltm::biserial.cor with sd.pop used instead of sd.
I think biserial.cor.new will return the same result as cor.test.
sd.pop = function(x){sd(x)*sqrt((length(x)-1)/length(x))}

biserial.cor.new = 
function (x, y, use = c("all.obs", "complete.obs"), level = 1) 
{
    if (!is.numeric(x)) 
        stop("'x' must be a numeric variable.\n")
    y <- as.factor(y)
    if (length(levs <- levels(y)) > 2) 
        stop("'y' must be a dichotomous variable.\n")
    if (length(x) != length(y)) 
        stop("'x' and 'y' do not have the same length")
    use <- match.arg(use)
    if (use == "complete.obs") {
        cc.ind <- complete.cases(x, y)
        x <- x[cc.ind]
        y <- y[cc.ind]
    }
    ind <- y == levs[level]
    diff.mu <- mean(x[ind]) - mean(x[!ind])
    prob <- mean(ind)
    diff.mu * sqrt(prob * (1 - prob))/sd.pop(x)
}

And an example:
x = c(3,4,5,6,7,5,6,7,8,9)
y = c(0,0,0,0,0,1,1,1,1,1)

library(ltm)

### DIFFERENT RESULTS WITH ltm::biserial.cor

biserial.cor(x,y, level=2)

   ### [1] 0.5477226

cor.test(x,y)

   ### Pearson's product-moment correlation
   ### sample estimates:
   ###       cor 
   ### 0.5773503

### SAME RESULTS WITH new function

biserial.cor.new(x,y, level=2)

    ### [1] 0.5773503

cor.test(x,y)

   ### Pearson's product-moment correlation
   ### sample estimates:
   ###       cor 
   ### 0.5773503

A: For my understanding you don't have to code the dichotome variable with 0 and 1. Therefore using other values results in exactly the same output. Try for example:
x <- 1:100
y <- rep(c(0,1), 50)
y2 <- rep(c(-786,345), 50)
cor.test(x, y)
cor.test(x, y2)

Both gives you an r of 0.01732137. The only thing that can happen by coding the dichotome variable differently is that you get -0.01732137, which will be the case if the first number is bigger than the second, e.g.
y3 <- rep(c(0,1), 50)
cor.test(x, y3)

results in -0.01732137.
Furthermore, I read on different pages that "the point-biserial correlation is equivalent to calculating the Pearson correlation between a continuous and a dichotomous variable", but in fact I get different results if I conduct a Pearson and a point-biserial correlation on same data.
An example:
x <- 1:100
y <- rep(c(0,1), 50)

cor.test(x, y)

gives me 0.01732137, but
biserial.cor(x, y) results in -0.01723455.
I understand that it is okay to get positive and negative values, but the absolute value should be the same, which is not the case. The results are also different if I use other data, e.g. x <- rnorm(100, 100, 15) instead of x <- 1:100.
For this reason I am unsure whether it is acceptable to use cor.test() and report that you have conducted a point-biserial correlation.
