Statistical Methods with R

Correlation & Inference

Unit B · Chapter 05 · Lecture 05b

Developed by Jeffrey M. Girard

Introductions

Variable relationships

So far, we have focused on the mean of one variable
But we are often interested in variable relationships
- As one variable changes, does the other tend to change?
- People with more anxiety tend to have more depression
- Nations with more museums tend to have lower mortality
We can quantify such relationships using correlations

Correlations

A correlation \((\rho,r)\) quantifies the direction and strength of the linear relationships between two variables
Direction: do the variables change together or opposing?
- \(\rho>0\): the variables change together
- \(\rho<0\): the variables change in opposition
Strength: how well does one variable predict the other?
- \(\rho \to 0\): no linear (or nonlinear) relationship
- \(\rho \to \pm1\): perfect (positive or negative) relationship

Positive linear

Negative linear

No linear / Nonlinear

Calculating
Correlations

Variance and covariance

A single variable’s variance measures its spread/changes

\[\text{Var}(x) = \frac{1}{n-1}\sum_{i=1}^n(x_i - \bar{x})(x_i - \bar{x})\]

Two variables’ covariance measures their linked changes

\[\text{Cov}(x,y) = \frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})\]

Pearson’s correlation

In general, we standardize something by dividing it by its SD
Here we need to divide by the product of the SDs of \(x\) and \(y\)

\[r_{xy} = \frac{\text{Cov}(x,y)}{s_x s_y}\]

This is called the Pearson Correlation Coefficient (or PCC)
- Also called the product-moment or inter-class correlation
- The population parameter is \(\rho\), the sample statistic is \(r\)

Example dataset

library(tidyverse)
cigarettes <- read_csv("../../data/cigarettes.csv")

state	log_packs	log_price	log_income
AL	4.96213	0.20487	4.64039
AZ	4.66312	0.16640	4.68389
AR	5.10709	0.23406	4.59435
CA	4.50449	0.36399	4.88147
CT	4.66983	0.32149	5.09472
DE	5.04705	0.21929	4.87087
DC	4.65637	0.28946	5.05960
FL	4.80081	0.28733	4.81155
GA	4.97974	0.12826	4.73299
ID	4.74902	0.17541	4.64307
IL	4.81445	0.24806	4.90387
IN	5.11129	0.08992	4.72916
IA	4.80857	0.24081	4.74211
KS	4.79263	0.21642	4.79613
KY	5.37906	-0.03260	4.64937
LA	4.98602	0.23856	4.61461
ME	4.98722	0.29106	4.75501
MD	4.77751	0.12575	4.94692
MA	4.73877	0.22613	4.99998
MI	4.94744	0.23067	4.80620
MN	4.69589	0.34297	4.81207
MS	4.93990	0.13638	4.52938
MO	5.06430	0.08731	4.78189
MT	4.73313	0.15303	4.70417
NE	4.77558	0.18907	4.79671
NV	4.96642	0.32304	4.83816
NH	5.10990	0.15852	5.00319
NJ	4.70633	0.30901	5.10268
NM	4.58107	0.16458	4.58202
NY	4.66496	0.34701	4.96075
ND	4.58237	0.18197	4.69163
OH	4.97952	0.12889	4.75875
OK	4.72720	0.19554	4.62730
PA	4.80363	0.22784	4.83516
RI	4.84693	0.30324	4.84670
SC	5.07801	0.07944	4.62549
SD	4.81545	0.13139	4.67747
TN	5.04939	0.15547	4.72525
TX	4.65398	0.28196	4.73437
UT	4.40859	0.19260	4.55586
VT	5.08799	0.18018	4.77578
VA	4.93065	0.11818	4.85490
WA	4.66134	0.35053	4.85645
WV	4.82454	0.12008	4.56859
WI	4.83026	0.22954	4.75826
WY	5.00087	0.10029	4.71169

Scatterplots

Estimation in R

Covariance

x <- cigarettes$log_packs
y <- cigarettes$log_price
n <- nrow(cigarettes)

cov_xy <- 1 / (n - 1) * sum(
  (x - mean(x)) * (y - mean(y)))
cov_xy
## [1] -0.008910294

Correlation

s_x <- sd(x)
s_y <- sd(y)

r_xy <- cov_xy / (s_x * s_y)
r_xy
## [1] -0.5397069

Shortcut Functions:

cov_xy <- cov(x, y)
cov_xy
## [1] -0.008910294

r_xy <- cor(x, y)
r_xy
## [1] -0.5397069

Interpretation

\(r=0\) is no relationship, \(r=\pm1\) is a perfect relationship
How should we interpret correlations between \(0\) and \(\pm1\)?
It depends on the context (what is common in your field?)
But here are some heuristics to use as a starting place

\((0.0,0.1)\)	\([0.1,0.3)\)	\([0.3,0.5)\)	\([0.5,1.0)\)
Negligible	Small	Medium	Large

Ansombe’s quartet

All four data sets have the same correlation \((r=.816)\)
But they show very different relationships qualitatively!

Another effect size

We can square \(r\) to get the coefficient of determination \((r^2)\)
How much variability is explained by the linear relationship?

\(r\)	\(r^2\)	Explained
0.1	0.01	1%
0.3	0.09	9%
0.5	0.25	25%
0.7	0.49	49%
0.9	0.81	81%

Obligatory PSA

Finding a correlation does not prove any causal story
There are many ways for \(x\) and \(y\) to become correlated
- Maybe \(x\) causes \(y\) (theorized model)
- Maybe \(y\) causes \(x\) (reverse causation)
- Maybe \(x\) and \(y\) cause each other (bidirectional)
- Maybe \(z\) causes both \(x\) and \(y\) (third variable)
- Maybe it is just sampling error (spurious)

Statistical Inference

Sampling error simulation

Samples from \(\rho=0.0\)

Samples from \(\rho=0.7\)

Statistical inference

I tend to treat \(r\) as a descriptive summary statistic
- I estimate it and maybe its CI but don’t use NHST
- I use more sophisticated methods to test relationships
But we can test the significance of a correlation
- Usually the null hypothesis is that it equals zero

\[ H_0: \rho=0\\ H_1: \rho\neq0 \]

Test statistic

Test statistics incorporate the effect size and sample size

\[t_{obs}=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}\]

What is its sampling distribution under the null?
- The frequentist parameteric approach assumes:

\[t_{obs} \sim t(\textit{df}=n-2, \mu=0, \sigma=1)\]

Significance testing

r_xy
## [1] -0.5397069

t_obs <- (r_xy * sqrt(n - 2)) / sqrt(1 - r_xy^2)
t_obs
## [1] -4.252537

t_crit <- qt(c(0.025, 0.975), df = n - 2)
t_crit
## [1] -2.015368  2.015368

Because \(t_{obs}\) is more extreme than \(t_{crit}\), we reject \(H_0\)
Thus, we conclude that \(r\neq0\) and specifically that \(r<0\)

Confidence intervals

Because \(r\) is bounded \([0,1]\), its CI is often asymmetrical
So we need a more complicated procedure to estimate its CI

Transform \(r\) to \(z\),
which is unbounded
Calculate a CI around \(z\),
which is symmetrical
Transform the bounds
of the \(z\) CI back to \(r\)

Confidence intervals

Step 1

Transform \(r\) to \(z\)

\(z = \frac{1}{2}\ln\left(\frac{1+r}{1-r}\right)\)

Step 2

Estimate \(z\) CI

\(z \pm \sqrt{\frac{1}{n-3}} \times 1.96\)

Step 3

Transform CI to \(r\)

\(r=\frac{e^{2z}-1}{e^{2z}+1}\)

r_xy
## [1] -0.5397069

z <-  (1 / 2) * log((1 + r_xy) / (1 - r_xy), base = exp(1))
z
## [1] -0.603742

zmult <- qnorm(0.975)
ci_z <- c(lo = z - sqrt(1 / (n - 3)) * zmult, hi = z + sqrt(1 / (n - 3)) * zmult)
ci_z
##         lo         hi 
## -0.9026337 -0.3048503

ci_r <- (exp(2 * ci_z) - 1) / (exp(2 * ci_z) + 1)
ci_r
##         lo         hi 
## -0.7175778 -0.2957450

A simple function

The cor.test() function will do all this work for you
- Note that our results are identical to before

cor.test(x, y)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = -4.2525, df = 44, p-value = 0.0001085
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7175778 -0.2957450
## sample estimates:
##        cor 
## -0.5397069

Prettier output

The correlation() function is more convenient

library(easystats)
correlation(cigarettes, select = "log_packs", select2 = "log_price")
## # Correlation Matrix (pearson-method)
## 
## Parameter1 | Parameter2 |     r |         95% CI | t(44) |         p
## --------------------------------------------------------------------
## log_packs  |  log_price | -0.54 | [-0.72, -0.30] | -4.25 | < .001***
## 
## p-value adjustment method: Holm (1979)
## Observations: 46

Many correlations

r_all <- correlation(cigarettes)
r_all
## # Correlation Matrix (pearson-method)
## 
## Parameter1 | Parameter2 |     r |         95% CI | t(44) |         p
## --------------------------------------------------------------------
## log_packs  |  log_price | -0.54 | [-0.72, -0.30] | -4.25 | < .001***
## log_packs  | log_income | -0.17 | [-0.44,  0.13] | -1.14 | 0.262    
## log_price  | log_income |  0.49 | [ 0.24,  0.68] |  3.75 | 0.001**  
## 
## p-value adjustment method: Holm (1979)
## Observations: 46

This also applies the Holm method to adjust the \(p\)-values
- We will talk more about this in a future lecture

Correlation matrix

r_mat <- summary(r_all)
r_mat
## # Correlation Matrix (pearson-method)
## 
## Parameter | log_income | log_price
## ----------------------------------
## log_packs |      -0.17 |  -0.54***
## log_price |     0.49** |          
## 
## p-value adjustment method: Holm (1979)

Plotting a matrix

plot(r_mat)

Correlation & Inference

Introductions

Variable relationships

Correlations

Positive linear

Positive linear

Negative linear

Negative linear

No linear / Nonlinear

No linear / Nonlinear

No linear / Nonlinear

Calculating Correlations

Variance and covariance

More on covariance

Pearson’s correlation

Example dataset

Scatterplots

Estimation in R

Interpretation

Ansombe’s quartet

Another effect size

Obligatory PSA

Statistical Inference

Sampling error simulation

Statistical inference

Test statistic

Significance testing

Confidence intervals

Confidence intervals

Step 1

Step 2

Step 3

A simple function

Prettier output

Many correlations

Correlation matrix

Plotting a matrix

Calculating
Correlations