Statistical Methods with R

Correlation & Inference

Unit B · Chapter 05 · Lecture 05b

Developed by Jeffrey M. Girard

Roadmap: Hypothesis Testing

1. Introductions

2. Calculating Correlations

3. Statistical Inference

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 relationship (a nonlinear one may still exist)
  • \rho \to \pm1: perfect (positive or negative) relationship

Positive linear

Positive linear

Negative linear

Negative linear

No linear / Nonlinear

No linear / Nonlinear

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

More on covariance

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

• \text{Cov} is positive when x and y vary together

• \text{Cov} is negative when x and y vary in opposition

• The \text{Cov} thus quantifies linear relationships…

  • But it has weird units and an unbounded range
  • So we standardize it to create the correlation

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("cigarettes.csv")
cigarettes

# A tibble: 46 × 4
   state log_packs log_price log_income
   <chr>     <dbl>     <dbl>      <dbl>
 1 AL         4.96     0.205       4.64
 2 AZ         4.66     0.166       4.68
 3 AR         5.11     0.234       4.59
 4 CA         4.50     0.364       4.88
 5 CT         4.67     0.321       5.09
 6 DE         5.05     0.219       4.87
 7 DC         4.66     0.289       5.06
 8 FL         4.80     0.287       4.81
 9 GA         4.98     0.128       4.73
10 ID         4.75     0.175       4.64
# ℹ 36 more rows

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

Anscombe’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 parametric 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 \rho\neq0 and specifically that \rho<0

Confidence intervals

• Because r is bounded [-1,1], its CI is often asymmetrical

• So we need a more complicated procedure to estimate its CI

1. Transform r to z,
   which is unbounded

2. Calculate a CI around z,
   which is symmetrical

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