Bivariate Data and Relationships

• We looked at single variables for descriptive statistics
• Most uses of statistics in economics and business investigate relationships between variables

• We will begin with bivariate data for relationships between X and Y
• Immediate aim is to explore associations between variables, quantified with correlation and linear regression
• Later we want to develop more sophisticated tools to argue for causation

Bivariate Data: Spreadsheets I

econfreedom <- read_csv("../data/econfreedom.csv")
head(econfreedom)

## # A tibble: 6 x 6
##      X1 Country   ISO      ef    gdp continent
##   <dbl> <chr>     <chr> <dbl>  <dbl> <chr>    
## 1     1 Albania   ALB    7.4   4543. Europe   
## 2     2 Algeria   DZA    5.15  4784. Africa   
## 3     3 Angola    AGO    5.08  4153. Africa   
## 4     4 Argentina ARG    4.81 10502. Americas 
## 5     5 Australia AUS    7.93 54688. Oceania  
## 6     6 Austria   AUT    7.56 47604. Europe

• Rows are individual observations
• Columns are variables on all individuals

Bivariate Data: Spreadsheets II

econfreedom %>%
  glimpse()

## Observations: 112
## Variables: 6
## $ X1        <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
## $ Country   <chr> "Albania", "Algeria", "Angola", "Argentina", "Australi…
## $ ISO       <chr> "ALB", "DZA", "AGO", "ARG", "AUS", "AUT", "BHR", "BGD"…
## $ ef        <dbl> 7.40, 5.15, 5.08, 4.81, 7.93, 7.56, 7.60, 6.35, 7.51, …
## $ gdp       <dbl> 4543.0880, 4784.1943, 4153.1463, 10501.6603, 54688.445…
## $ continent <chr> "Europe", "Africa", "Africa", "Americas", "Oceania", "…

Bivariate Data: Spreadsheets III

source("../files/summaries.R") # use my summary_table function
econfreedom %>%
    summary_table(ef, gdp)

Bivariate Data: Scatterplots

library("ggplot2")
ggplot(data = econfreedom)+
  aes(x = ef, y = gdp)+
  geom_point(aes(color = continent),
             size = 2)+
  labs(x = "Economic Freedom Index (2014)",
       y = "GDP per Capita (2014 USD)",
       color = "")+
  theme_bw(base_family = "Fira Sans Condensed",
           base_size=20)+
  theme(legend.position = "bottom")

Associations

Covariance

• For any two variables, we can measure their sample covariance, cov(X,Y) or sX,Y to quantify how they vary together¹

sX,Y=E[(X−ˉX)(Y−ˉY)]

• Intuition: if X is above its mean, would we expect Y:
  • to be above its mean also (X and Y covary positively)
  • to be below its mean (X and Y covary negatively)

• Covariance is a common measure, but the units are meaningless, thus we rarely need to use it so don't worry about learning the formula

¹ Henceforth we limit all measures to samples, for convenience. Population covariance is denoted σX,Y

Covariance, in R

# base R 
cov(econfreedom$ef,econfreedom$gdp)

## [1] 8922.933

# dplyr 
econfreedom %>%
  summarize(cov = cov(ef,gdp))

## # A tibble: 1 x 1
##     cov
##   <dbl>
## 1 8923.

Correlation

• More convenient to standardize covariance into a more intuitive concept: correlation, ρ or r ∈[−1,1]

rX,Y=sX,YsXsY=cov(X,Y)sd(X)sd(Y)

• Simply weight covariance by the product of the standard deviations of X and Y

• Alternatively, take the average¹ of the product of standardized (Z-scores for) each (xi,yi) pair:²

r=1n−1n∑i=1(xi−ˉXsX)(yi−ˉYsY)r=1n−1n∑i=1ZXZY

Correlation: Interpretation

Guess the Correlation!

Correlation and Covariance in R

# Base r: cov or cor(df$x, df$y)
cov(econfreedom$ef, econfreedom$gdp)

## [1] 8922.933

cor(econfreedom$ef, econfreedom$gdp)

## [1] 0.5867018

# dplyr 
econfreedom %>%
  summarize(covariance = cov(ef, gdp),
            correlation = cor(ef, gdp))

## # A tibble: 1 x 2
##   covariance correlation
##        <dbl>       <dbl>
## 1      8923.       0.587

Correlation and Covariance in R I

• corrplot is a great package (install and then load) to visualize correlations in data

library(corrplot) # see more at https://github.com/taiyun/corrplot
library(RColorBrewer) # for color scheme used here
library(gapminder) # for gapminder data
# need to make a corelation matrix with cor(); can only include numeric variables
gapminder_cor<- gapminder %>%
  dplyr::select(gdpPercap, pop, lifeExp)
# make a correlation table with cor (base R)
gapminder_cor_table<-cor(gapminder_cor)
# view it
gapminder_cor_table

Correlation and Covariance in R II

corrplot(gapminder_cor_table, type="upper", 
         method = "circle", # number for showing correlation coefficient
         order="alphabet", 
         col=brewer.pal(n=8, name="RdBu"))

Correlation and Endogeneity

Always Plot Your Data!

Fitting a Line to Data

Fitting a Line to Data

Fitting a Line to Data

Population Linear Regression Model

• Linear regression lets us estimate the slope of the population regression line between X and Y

• We can make inferences about the population slope coefficient

• eventually, a causal interpretation

• slope=ΔYΔX: for a 1-unit change in X, how many units will this cause Y to change?

Class Size Example

Class Size Example: Load the Data

# install.packages("haven") # install for first use
library("haven") # load for importing .dta files
CASchool<-read_dta("../data/caschool.dta")

Class Size Example: Look at the Data I

glimpse(CASchool)

## Observations: 420
## Variables: 21
## $ observat <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ dist_cod <dbl> 75119, 61499, 61549, 61457, 61523, 62042, 68536, 63834,…
## $ county   <chr> "Alameda", "Butte", "Butte", "Butte", "Butte", "Fresno"…
## $ district <chr> "Sunol Glen Unified", "Manzanita Elementary", "Thermali…
## $ gr_span  <chr> "KK-08", "KK-08", "KK-08", "KK-08", "KK-08", "KK-08", "…
## $ enrl_tot <dbl> 195, 240, 1550, 243, 1335, 137, 195, 888, 379, 2247, 44…
## $ teachers <dbl> 10.90, 11.15, 82.90, 14.00, 71.50, 6.40, 10.00, 42.50, …
## $ calw_pct <dbl> 0.5102, 15.4167, 55.0323, 36.4754, 33.1086, 12.3188, 12…
## $ meal_pct <dbl> 2.0408, 47.9167, 76.3226, 77.0492, 78.4270, 86.9565, 94…
## $ computer <dbl> 67, 101, 169, 85, 171, 25, 28, 66, 35, 0, 86, 56, 25, 0…
## $ testscr  <dbl> 690.80, 661.20, 643.60, 647.70, 640.85, 605.55, 606.75,…
## $ comp_stu <dbl> 0.34358975, 0.42083332, 0.10903226, 0.34979424, 0.12808…
## $ expn_stu <dbl> 6384.911, 5099.381, 5501.955, 7101.831, 5235.988, 5580.…
## $ str      <dbl> 17.88991, 21.52466, 18.69723, 17.35714, 18.67133, 21.40…
## $ avginc   <dbl> 22.690001, 9.824000, 8.978000, 8.978000, 9.080333, 10.4…
## $ el_pct   <dbl> 0.000000, 4.583333, 30.000002, 0.000000, 13.857677, 12.…
## $ read_scr <dbl> 691.6, 660.5, 636.3, 651.9, 641.8, 605.7, 604.5, 605.5,…
## $ math_scr <dbl> 690.0, 661.9, 650.9, 643.5, 639.9, 605.4, 609.0, 612.5,…
## $ aowijef  <dbl> 35.77982, 43.04933, 37.39445, 34.71429, 37.34266, 42.81…
## $ es_pct   <dbl> 1.000000, 3.583333, 29.000002, 1.000000, 12.857677, 11.…
## $ es_frac  <dbl> 0.01000000, 0.03583334, 0.29000002, 0.01000000, 0.12857…

Class Size Example: Look at the Data II

Class Size Example: Scatterplot

scatter<-ggplot(data = CASchool)+
  aes(x = str,
      y = testscr)+
  geom_point(color = "blue")+
  labs(x = "Student to Teacher Ratio",
       y = "Test Score")+
  theme_bw(base_family = "Fira Sans Condensed",
           base_size = 20)
scatter

Class Size Example: Slope I

Class Size Example: Slope II

Class Size Example: Slope II

Class Size Example: Slope and Average Effect

Class Size Example: Marginal Effects

Econometric Models Overview

Y=β0+β1X1+β2X2+u

• Y is the dependent variable of interest
  • AKA "response variable," "regressand," "Left-hand side (LHS) variable"

• X1 and X2 are independent variables
  • AKA "explanatory variables", "regressors," "Right-hand side (RHS) variables", "covariates"

• Our data consists of a spreadsheet of observed values of (X1i,X2i,Yi)

• To model, we "regress Y on X1 and X2"

• β0 and β1 are parameters that describe the population relationships between the variables

  • unknown! to be estimated!

• u is the random error term
  • 'U'nobservable, we can't measure it, and must model with assumptions about it

The Population Regression Model

Deriving OLS

Deriving OLS

Deriving OLS

Deriving OLS

Deriving OLS

Deriving OLS

O rdinary L east S quares Estimators

• The Ordinary Least Squares (OLS) estimators of the unknown population parameters β0 and β1, solve the calculus problem:

minβ0,β1n∑i=1[Yi−(β0+β1Xi⏟^Yi)⏟u]2

• Intuitively, OLS estimators minimize the average squared distance between the actual values (Yi) and the predicted values (ˆYi) along the estimated regression line

The OLS Regression Line

• The OLS regression line or sample regression line is the linear function constructed using the OLS estimators:

^Yi=^β0+^β1Xi

• ^β0 and ^β1 ("beta 0 hat" & "beta 1 hat") are the OLS estimators of population parameters β0 and β1 using sample data

• The predicted value of Y given X, based on the regression, is E(Yi|Xi)=^Yi

• The residual or prediction error for the ith observation is the difference between observed Yi and its predicted value, ^ui=Yi−^Yi

The OLS Regression Estimators

• The solution to the SSE minimization problem yields:¹

ˆβ0=ˉY−^β1ˉX

ˆβ1=n∑i=1(Xi−ˉX)(Yi−ˉY)n∑i=1(Xi−ˉX)2=sXYs2X=cov(X,Y)var(X)

¹ See tomorrow's class notes page for proofs.

Class Size Scatterplot (Again)

scatter

Class SIze Scatterplot with Regression Line

scatter+
  geom_smooth(method = "lm", color = "red")

OLS in R

# run regression of testscr on str
school_reg <- lm(testscr ~ str, 
                 data = CASchool)

OLS in R II

# look at reg object
school_reg

## 
## Call:
## lm(formula = testscr ~ str, data = CASchool)
## 
## Coefficients:
## (Intercept)          str  
##      698.93        -2.28

OLS in R III

summary(school_reg) # get full summary

## 
## Call:
## lm(formula = testscr ~ str, data = CASchool)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.727 -14.251   0.483  12.822  48.540 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 698.9330     9.4675  73.825  < 2e-16 ***
## str          -2.2798     0.4798  -4.751 2.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared:  0.05124,    Adjusted R-squared:  0.04897 
## F-statistic: 22.58 on 1 and 418 DF,  p-value: 2.783e-06

Tidy OLS in R: broom I

# load packages
library(broom)
# tidy regression output
tidy(school_reg)

## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)   699.       9.47      73.8  6.57e-242
## 2 str            -2.28     0.480     -4.75 2.78e-  6

¹ See more at broom.tidyverse.org.

Tidy OLS in R: broom II

# load packages
library(broom)
# tidy regression output (with confidence intervals!)
tidy(school_reg,
     conf.int = TRUE)

## # A tibble: 2 x 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)   699.       9.47      73.8  6.57e-242   680.      718.  
## 2 str            -2.28     0.480     -4.75 2.78e-  6    -3.22     -1.34

¹ See more at broom.tidyverse.org.

More broom Tools: glance

• glance() shows us a lot of overall regression statistics and diagnostics
  • We'll interpret these in the next lecture and beyond

# look at regression statistics and diagnostics
glance(school_reg)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
## 1    0.0512        0.0490  18.6      22.6 2.78e-6     2 -1822. 3650. 3663.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

More broom Tools: augment

# add regression-based values to data
augment(school_reg)

## # A tibble: 420 x 9
##    testscr   str .fitted .se.fit .resid    .hat .sigma  .cooksd .std.resid
##      <dbl> <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>    <dbl>      <dbl>
##  1    691.  17.9    658.   1.24    32.7 0.00442   18.5 0.00689       1.76 
##  2    661.  21.5    650.   1.28    11.3 0.00475   18.6 0.000893      0.612
##  3    644.  18.7    656.   1.01   -12.7 0.00297   18.6 0.000700     -0.685
##  4    648.  17.4    659.   1.42   -11.7 0.00586   18.6 0.00117      -0.629
##  5    641.  18.7    656.   1.02   -15.5 0.00301   18.6 0.00105      -0.836
##  6    606.  21.4    650.   1.24   -44.6 0.00446   18.5 0.0130       -2.40 
##  7    607.  19.5    654.   0.909  -47.7 0.00239   18.5 0.00794      -2.57 
##  8    609   20.9    651.   1.09   -42.3 0.00343   18.5 0.00895      -2.28 
##  9    612.  19.9    653.   0.919  -41.0 0.00244   18.5 0.00597      -2.21 
## 10    613.  20.8    652.   1.07   -38.9 0.00329   18.5 0.00723      -2.09 
## # … with 410 more rows

Class Size Regression Result I

• Using OLS, we find: ^test score=689.9−2.28×str

Class Size Regression Result II

• There's a great package called equatiomatic that prints this equation in markdown or LATEX.

testscr=698.93−2.28(str)+ϵ

Here was my code:

# install.packages("equatiomatic") # install for first use
library(equatiomatic) # load it
extract_eq(school_reg, # regression lm object
           use_coefs = TRUE, # use names of variables
           coef_digits = 2, # round to 2 digits
           fix_signs = TRUE) # fix negatives (instead of + -)

## $$
## \text{testscr} = 698.93 - 2.28(\text{str}) + \epsilon
## $$

• In R chunk in R markdown, set {r, results="asis"} to print this raw output to be rendered

Class Size Regression: A Data Point

CASchool %>%
  filter(district=="Richmond Elementary") %>%
  dplyr::select(district, testscr, str)

## # A tibble: 1 x 3
##   district            testscr   str
##   <chr>                 <dbl> <dbl>
## 1 Richmond Elementary    672.    22