How t is calculated for pearson cor test?

Question

I am trying to understand how the t-statistic is calculated by the cor.test() function in R.

R> cor.test(seq(from=1, to=9, len=5), seq(from=2, to=10, len=5))

    Pearson's product-moment correlation

data:  seq(from = 1, to = 9, len = 5) and seq(from = 2, to = 10, len = 5)
t = 82191237, df = 3, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 1 1
sample estimates:
cor
  1

Wiki says the relation between t and r is the following:

$$ t = r \sqrt {{n-2}\over{1-r^2}} $$

But when r is 1, the above formula is not defined. What formula is used to compute the t statistic in the cor.test function?

Answer 1

Here is the relevant part from the code of cor.test():

r  <- cor(x, y)
df <- n - 2L
ESTIMATE  <- c(cor = r)
PARAMETER <- c(df = df)
STATISTIC <- c(t = sqrt(df) * r / sqrt(1 - r^2))

Which, when you run it with your generated sequences, produce STATISTIC equal to 82191237.

This happens because the computed r value is not exactly equal to 1:

r <- cor(seq(from=1, to=9, len=5), seq(from=2, to=10, len=5))

r == 1
[1] FALSE

print(r, digits=20)
[1] [1] 0.99999999999999977796

Plug in those numbers in the t-value calculation and you will get the observed result:

sqrt(3) * r / sqrt(1 - r^2)
[1] 82191237

This happens because of floating point representation errors. If the computed correlation would be exactly equal to one you would get t value equal to infinity due to division by zero:

sqrt(3) * 1 / sqrt(1 - 1^2)
[1] Inf

Now, R is seems a bit inconsistent with machine precision corrections. Some tests round the value if the difference is smaller than maximum machine precision, while others don't. My feeling is that more popular tests often get more rigorous treatment. So for example in the code of a t.test() you would find:

if(stderr < 10 *.Machine$double.eps * abs(mx))
  stop("data are essentially constant")

One of my side projects is re-writing R tests to work on rows/columns of matrices, and, as an example, in that version I added a check for correlation being indistinguishable from 1 (within the range of machine precision rounding) with the following substitutions:

rs[abs(rs - 1) < .Machine$double.eps^0.5] <- 1
rs[abs(rs + 1) < .Machine$double.eps^0.5] <- -1

Which then in this specific case returns the correct result:

matrixTests::row_cor_pearson(seq(from=1, to=9, len=5), seq(from=2, to=10, len=5))

  obs.paired cor df statistic pvalue conf.low conf.high alternative cor.null conf.level
1          5   1  3       Inf      0        1         1   two.sided        0       0.95

But you have to keep in mind that rounding machine precision errors is a double-edged thing. You cannot distinguish between two possibilities and have to favour one over the other. My reasoning was that correlation values extremely close to 1 will much more frequently occur due to low sample sizes, rather than any other cause (like using huge extremely correlated samples).