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:

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?

  • $\begingroup$ When $r=1$ you can take $t$ to be arbitrarily large. $\endgroup$
    – whuber
    Aug 1 at 17:05
  • $\begingroup$ But how r calculates it so that it is a finite but very large number? $\endgroup$ Aug 1 at 17:17
  • 1
    $\begingroup$ R appears to exploit limitations of double-precision floating point representation to conclude $t=82191237:$ you can take this to be an arbitrarily large value. The correct p-value, BTW, is $0$ (exactly). $\endgroup$
    – whuber
    Aug 1 at 17:26
  • 1
    $\begingroup$ It's not wrong when you understand that it's working with double-precision floats. That's one reason the p-value is printed as "<2.2e-16": that's one of the smallest values distinguishable from $1.$ You also know that any huge t-value like this one is going to be associated with an extremely small (tail) probability, no matter how many degrees of freedom it has. If you need exact computation (which, for practical applications, you don't), then use an appropriate mathematical platform like Mathematica. $\endgroup$
    – whuber
    Aug 1 at 17:37
  • 1
    $\begingroup$ Look at the source code. R computes t from the correlation of the inputs. That, in turn, is computed by the cor function. Because of floating point rounding error, the correlation of the inputs is not exactly $1:$ it is a tiny bit less. Thus, t is not computed as infinity. $\endgroup$
    – whuber
    Aug 1 at 18:42

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

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.