Testing the significance of the correlation coefficient r

The sample correlation coefficient is a measure of the strength of the linear relationship between two variables X and Y. Actually it is the coefficient of determination $r^2$ which is a better measure. It is computed by the formula

$$ r =\frac{\sum_{i = 1}^n (x_i - \overline{x} ) (y_i - \overline{y})}{S_x S_y} $$

where $S_x$ and $S_y$ are the standard deviations of the X and Y samples. We have For $S_x = \sqrt{\frac{\sum (x_i -\overline{x})^2}{n-1}}$ with a similar formula for $S_y$. The value of $r$ varies from -1 (perfect negative correlation) to +1 (perfect positive correlation). A negative correlation means that as one variable increases, the other other variable decreases, or in other words, both variables varies in opposite directions. On the other hand a positive correlation arises when both variables varies together in the same direction: if one variable increases, the other also increases and if one variable decreases, the other also decreases.

Correlation is generally used when it is difficult to determine the causative order. Does X causes Y or does Y causes X? It is heavily used in the Psychology field for example. In this post, we describe how to test the significance of the computed correlation coefficient.

The sample test statistic $t= \frac{r-\rho}{\sqrt{\frac{1-r^2}{n-2}}}$ follow assymptotically the t-distribution with $df=n-2$ degree of freedom and tables of the t-distribution may be used. Thus to write our Python program to use the ttest routine described previously in our blog, we specify the sample test statistic and the standard error.

def corrtest( rho, r, n, alpha = 0.5, side = 0):
     stattest = (r - rho)
     se = sqrt( (1-r*r)/(n-2.0))
     return ttest(r-rho, se, alpha, side)

Neat, is not it? We obtain the p-value, the distribution test statistic and the critical limits(values). For side = 0, a double sided test, the two limits do not strictly form the true confidence interval. Here is the complete testforcorr.py Python file where we included the ttest function. It calls the rstats.py module we have written before to use R-names for the Python statistical functions.

# -*- coding: utf-8 -*-
"""
file     testofcorr.py
author   Ernesto P. Adorio
         ernesto.adorio @ gmail.com
         UP at Clarkfield
         Angeles, Pampanga

version  0.0.1 april 4, 2010
"""

from math import *
from rstats import *



def ttest(samplestat, se, df, alpha= 0.05, side=0):
    """                                            
    T-test of a sample statistic.                  
    Arguments:                                     
     samplestat- sample statistic                  
     se -    standard error of sample statistic    
     df      degree of freedom                     
     alpha - significance level                    
     side  - -1 left-sided test                    
             0  double sided test                  
             +1 right-sided test                   
    """                                            
    Ttest = samplestat/se                          
    if side ==0:                                   
        pvalue = pt(Ttest, df)                     
        if Ttest > 0:                              
            pvalue = 1-pvalue                      
        pvalue *=2                                 
        tcrit1 =    qt(alpha/2,   df)              
        tcrit2 =    qt(1-alpha/2.0,df)             
        return pvalue, Ttest, (tcrit1,tcrit2)      
    elif side == -1:                               
       pvalue = pt(Ttest,  df)                     
       tcrit = qt(alpha,  df)                      
       return pvalue, Ttest, tcrit                 
    else:                                          
       pvalue = 1- pt(Ttest,  df)                  
       tcrit  = qt(1.0-alpha,  df)                 
       return pvalue, Ttest, tcrit                 


def cov(X,Y):
    n = len(X)
    if n != len(Y):
        raise "ArgumentError", "in cov: len(X) != len(Y)"
    xbar = float(sum(X))/ n
    ybar = float(sum(Y))/n
    return sum([(x -xbar)*(y-ybar) for x,y in zip(X,Y)])/(n-1)

def Sx(X):
    xbar = float(sum([x for x in X]))/n
    return sqrt( sum([(x - xbar)**2 for x in X]) /(n-1.0))
 
def cor(X,Y): 
    # correlation coefficient of X and Y. 
    return cov(X,Y)/(Sx(X)*Sx(Y))


def corrtest( rho, r, n, alpha = 0.5, side = 0):
     stattest = (r - rho)
     se = sqrt( (1-r*r)/(n-2.0))
     return ttest(r-rho, se, n-2,  alpha, side)


def xycorrtest(X,Y, rho, alpha = 0.5, side= 0):
     if len(X) != len(Y):
        raise ArgumentError, "in xycorrtest: X and Y must have same length"
     r = cor(X,Y)
     return corrtest(rho, r, len(X), alpha, side)


def Test():
     #Berenson , p. 546 
     rho = 0
     r = 0.951
     alpha = 0.05
     side = 1
     n = 14
     print "Berenson's example, p. 546"
     print corrtest( rho, r, n, alpha, side)

if __name__ == "__main__":
    Test()

When we invoke the Python interpreter it will run the Berenson, p. 546 example. The output is


toto@toto-laptop:~/Blogs/statistics$ python testofcorr.py
Berenson's example, p. 546
(8.9865653363219167e-08, 10.654779470653486, 1.7822875556491591)


The extremely low p-value means that we can reject the Null hypothesis that the correlation coefficient is zero.
To support this decision, the obtained distribution test statistic is greater then 1.78 critical value limit, i.e.,
is in the rejection region.


I need to specify a darker background for <code> blocks and a central download repository for
the files. Stay tuned. We will add more useful technical contents soon!