Pearson’s Product Moment Correlation
The Pearson product moment correlation coefficient r estimates population correlation \(\rho\). The test evaluates \(H_0: \rho = 0\) and constructs a \((1 - \alpha)\%\) confidence interval around \(r\):
\(r = \frac{\sum_{i=1}^n{(X_i - \bar{X})(Y_i - \bar{Y})}}{\sqrt{\sum_{i=1}^n{(X_i - \bar{X})^2 \sum_{i=1}^n{(Y_i - \bar{Y})^2}}}}\)
\(r\) ranges from \(-1\) (perfect negative linear relationship) to \(+1\) (perfect positive linear relationship, and \(r = 0\) when there is no linear relationship. A correlation in the range \((.1, .3)\) is small, \((.3, .5)\) is medium, and \((.5, 1.0)\) is large.
Only use the Pearson correlation when
- the variables are interval or ratio,
- the variables are normally distributed,
- the variables are linearly related,
- there are minimal outliers (or they are removed entirely), and
- the residuals are homoscedastic.
Test \(H_0: \rho = 0\) with test statistic \(T = r \sqrt{\frac{n-2}{1-r^2}}\). The test also defines a \((1 - \alpha)\%\) confidence interval.3 I found an explanation of the confidence interval here.
Pearson’s Product Moment Correlation Example
The nascard data set consists of 898 races from 1975 - 2003. In race 1 of 1975, what was the correlation between the driver’s finishing position and prize?
library(dplyr) # for glimpse()
nascard <- read.fwf(file = url("http://jse.amstat.org/datasets/nascard.dat.txt"),
widths = c(5, 6, 4, 4, 4, 5, 9, 4, 11, 30),
col.names = c('series_race', 'year', 'race_year',
'finish_pos', 'start_pos',
'laps_comp', 'winnings', 'num_cars',
'car_make', 'driver'))
nascard_sr1 <- nascard[nascard$series_race == 1,]
glimpse(nascard_sr1)## Observations: 35
## Variables: 10
## $ series_race <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ year <dbl> 1975, 1975, 1975, 1975, 1975, 1975, 1975, 1975, 19...
## $ race_year <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ finish_pos <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,...
## $ start_pos <dbl> 1, 2, 25, 27, 22, 35, 3, 21, 33, 10, 17, 24, 6, 28...
## $ laps_comp <dbl> 191, 191, 184, 183, 178, 175, 172, 168, 167, 166, ...
## $ winnings <dbl> 12035, 8135, 6535, 5035, 3835, 2885, 3185, 2485, 2...
## $ num_cars <dbl> 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35...
## $ car_make <fct> Matador , Mercury , Chevrolet , Dodge ...
## $ driver <fct> BobbyAllison , DavidPearson ...Explore the relationship with a scatterplot. As one would expect, there is a negative relationship between finish position and winnings, but it is not linear.
library(ggplot2)
ggplot(data = nascard_sr1, aes(x = finish_pos, y = winnings)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
However, 1/winnings may be linearly related to finish position.4 I tried several variable tranformations from Stat Trek, and chose the tranformation with the highest \(R^2\) that also had a normal distribution in each variable.
ggplot(data = nascard_sr1, aes(x = finish_pos, y = 1/winnings)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)
#summary(lm(winnings ~ finish_pos, data = nascard_sr1)) # r2 = .5474
#summary(lm(log(winnings) ~ finish_pos, data = nascard_sr1)) # r2 = .9015
summary(lm(1/winnings ~ finish_pos, data = nascard_sr1)) # r2 = .9509##
## Call:
## lm(formula = 1/winnings ~ finish_pos, data = nascard_sr1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.343e-04 -7.449e-05 -1.710e-06 4.575e-05 3.623e-04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.829e-05 3.682e-05 -0.497 0.623
## finish_pos 4.512e-05 1.784e-06 25.292 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0001066 on 33 degrees of freedom
## Multiple R-squared: 0.9509, Adjusted R-squared: 0.9495
## F-statistic: 639.7 on 1 and 33 DF, p-value: < 2.2e-16#summary(lm(winnings ~ log(finish_pos), data = nascard_sr1)) # r2 = .883
#summary(lm(log(winnings) ~ log(finish_pos), data = nascard_sr1)) # r2 = .9766
#summary(lm(winnings ~ log(finish_pos), data = nascard_sr1)) # r2 = .883The finish position and prize are quantitative ratio variables, so the Pearson’s product moment correlation should apply. First check whether each variable is normally distributed. The normal distribution plots look good and the Anderson-Darling tests both fail to reject the normality null hypothesis.
qqnorm(nascard_sr1$finish_pos)
qqline(nascard_sr1$finish_pos)
qqnorm(1/(nascard_sr1$winnings))
qqline(1/(nascard_sr1$winnings))
library(nortest) # for Anderson-Darling
ad.test(nascard_sr1$finish_pos)##
## Anderson-Darling normality test
##
## data: nascard_sr1$finish_pos
## A = 0.37344, p-value = 0.3989ad.test(1/nascard_sr1$winnings)##
## Anderson-Darling normality test
##
## data: 1/nascard_sr1$winnings
## A = 0.3292, p-value = 0.505The Pearson product moment correlation is \(r = \frac{\sum_{i=1}^n{(X_i - \bar{X})(Y_i - \bar{Y})}}{\sqrt{\sum_{i=1}^n{(X_i - \bar{X})^2 \sum_{i=1}^n{(Y_i - \bar{Y})^2}}}} = .9752\)
x <- nascard_sr1$finish_pos
y <- 1/nascard_sr1$winnings
(r = sum((x - mean(x)) * (y - mean(y))) /
sqrt(sum((x - mean(x))^2) * sum((y - mean(y))^2)))## [1] 0.9751633The test statistic is \(T = r \sqrt{\frac{n-2}{1-r^2}} = .9752 \sqrt{\frac{35-2}{1-.9752^2}} = 25.292\). \(P(T>.9752) < .0001\), so reject \(H_0\) that the correlation is zero.
n <- nrow(nascard_sr1)
(t = r * sqrt((n-2) / (1 - r^2)))## [1] 25.29218pt(q = t, df = n-2, lower.tail = FALSE) * 2## [1] 3.516944e-23R function cor.test performs these calculations. Specify method = "pearson".
cor.test(x = nascard_sr1$finish_pos,
y = 1/nascard_sr1$winnings,
alternative = "two.sided",
method = "pearson")##
## Pearson's product-moment correlation
##
## data: nascard_sr1$finish_pos and 1/nascard_sr1$winnings
## t = 25.292, df = 33, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9509453 0.9875017
## sample estimates:
## cor
## 0.9751633
