Turner_Discussion_Heteroskedasticity

1 Setup

# Clear the workspace
  rm(list = ls()) # Clear environment
  gc()            # Clear unused memory

##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 520346 27.8    1152055 61.6         NA   669294 35.8
## Vcells 957707  7.4    8388608 64.0      16384  1839854 14.1

  cat("\f")       # Clear the console

  dev.off         # Clear the charts

## function (which = dev.cur()) 
## {
##     if (which == 1) 
##         stop("cannot shut down device 1 (the null device)")
##     .External(C_devoff, as.integer(which))
##     dev.cur()
## }
## <bytecode: 0x7fda4f474410>
## <environment: namespace:grDevices>

#Install packages
  packages <- c("knitr",
                "stargazer",
                "MASS",
                "tidyverse",
                "psych",
                "skedastic", # INSTALLS THE SKEDASTIC PACKAGE
                "lmtest",
                "dplyr",
                "ggplot2",
                "wooldridge",
                "sandwich",
                "AER")

  for (i in 1:length(packages)) {# Iterate through each package in the 'packages' vector.
    if (!packages[i] %in% rownames(installed.packages())) {#check if package is installed
      install.packages(packages[i], #install package if it isn't already
                       repos = "http://cran.rstudio.com/",
                       dependencies = TRUE)}
    library(packages[i], character.only = TRUE)} # Load the package into the R session.

## 
## Please cite as:

##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Attaching package: 'psych'
## 
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## 
## Loading required package: zoo
## 
## 
## Attaching package: 'zoo'
## 
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 
## 
## Attaching package: 'wooldridge'
## 
## 
## The following object is masked from 'package:MASS':
## 
##     cement
## 
## 
## Loading required package: car
## 
## Loading required package: carData
## 
## 
## Attaching package: 'car'
## 
## 
## The following object is masked from 'package:psych':
## 
##     logit
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## 
## Loading required package: survival

# Remove i and packages
  rm(i, packages)

2 Part I: Issue Summary

2.1 What is Heteroskedasticity?

Heteroskedasticity is the statistical issue where the variance of the errors is unequal throughout a regression model. A solid example of this would be age and income because most younger people make around the same lower level of income but, as they get older, the range between highest earners and lowest earners is likely to grow substantially. This is problematic as it violates the assumption of homoskedasticity in a classical linear regression and it can impact the accuracy of standard errors, leading to unreliable hypothesis tests and confidence intervals. In more plain terms, we know longer know if we are using the best, linear, unbiased estimator because we cannot trust the results of the significance tests.

2.2 Hypothesis Testing for Heteroskedasticity:

2.2.1 Breusch-Pagan Test:

Aux Regression:

\[ \hat{u}^2=\gamma_0+\gamma_1X_1+\gamma_2X_2 +...+\gamma_kX_k+v \]

• \(\hat{u}^2\) = residuals from original regression

• \(\gamma_i\) = estimators in aux regression

• \(X_i\) = value for \(\beta_i\) in original regression

• \(v\) = residuals in aux regression

Null Hypothesis: Consistent variance throughout model; no heteroskedasticity

• All of the betas in the auxiliary regression should = 0

\[ \sum_{i=1}^k\gamma_i=0 \]

Alternative: Heteroskedasticity in model

• There is at least one estimator that is not equal to 0

\[ \sum_{i=1}^k\beta_i\not=0 \]

2.2.2 White Test:

Aux Regression:

\[ \hat{u}^2=\gamma_0+\gamma_1X_1+\gamma_2X_2+...+\gamma_kX_k+ \delta_1 \epsilon^2_1 + \delta_2 \epsilon^2_2+...+\delta_k \epsilon^2_k+v \]

• \(\hat{u}^2\) = residuals from original regression

• \(\gamma_i\) = estimators for \(X_i\) aux regression

• \(X_i\) = value for \(\beta_i\) in original regression

• \(\delta_i\) = estimators for \(\epsilon_i\) in original regression

• \(\epsilon^2_i\) = squared residuals of \(X_i\) (the explanatory variable) in the original reg

• \(v\) = residuals in aux regression

Null Hypothesis: Consistent variance and no serial correlation throughout model; no heteroskedasticity

\[ \sum_{i=1}^k\gamma_i= \sum_{i=1}^k\delta_i=0 \]

Alternative: Heteroskedasticity in model

\[ \sum_{i=1}^k\gamma_i\not=0 \space or \space \sum_{i=1}^k\delta_i\not=0 \]

2.2.3 Comparison of White and Breusch-Pagan (BP)

In both the White and BP tests, the conclusion from the hypothesis is either “no evidence of heteroskedasticity” or “there is evidence for heteroskedasticity.” The key difference is that the white test additionally checks for serial correlation. While the specific details of the auxiliary regression may differ slightly between the two tests, but their underlying logic is generally in line with assessing the homoscedasticity assumption in regression analysis.

3 Part 2: Coding

3.1 Selecting Package and cleaning

# Use credit card data from AER 
  data("CPS1988")
  wages <- CPS1988
  rm(CPS1988)
  
# Examine data set
  str(wages)

## 'data.frame':    28155 obs. of  7 variables:
##  $ wage      : num  355 123 370 755 594 ...
##  $ education : int  7 12 9 11 12 16 8 12 12 14 ...
##  $ experience: int  45 1 9 46 36 22 51 34 0 18 ...
##  $ ethnicity : Factor w/ 2 levels "cauc","afam": 1 1 1 1 1 1 1 1 1 1 ...
##  $ smsa      : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ region    : Factor w/ 4 levels "northeast","midwest",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ parttime  : Factor w/ 2 levels "no","yes": 1 2 1 1 1 1 1 1 1 1 ...

#  Convert ethnicity to numeric
  wages$ethnicity <- ifelse(wages$ethnicity == "afam", 0, 1)
                                
# Show summary stats of numeric data
  stargazer(wages, type = "text")

## 
## ====================================================
## Statistic    N     Mean   St. Dev.  Min      Max    
## ----------------------------------------------------
## wage       28,155 603.727 453.547  50.050 18,777.200
## education  28,155 13.068   2.900     0        18    
## experience 28,155 18.200   13.079    -4       63    
## ethnicity  28,155  0.921   0.270     0        1     
## ----------------------------------------------------

Wage data was chosen for its high potential to show heteroskedasticity. Specifically CPS1988 because it takes observations across industries.

Additionally, the summary stats do not reveal any significant problems in the data. It is odd that there is a high frequency of observations where the wage per week is below minimum wage in 1988 ($3.35/hr or $134/week). However, some jobs in the service industry can have tip based wages and be paid below minimum, so they will be included in model. Lastly, the -4 value for experience is illogical and will be omitted

3.2 Variable Definitions

Variable	Definition
wage	Wage (in dollars per hour)
education	Number of years of education.
experience	Number of years of potential work experience
ethnicity	Factor with levels "cauc" and "afam" (African-American).

3.3 Main Regression

\[ wage = \beta_0 + \beta_1 education_i + \beta_2 experience_i + \beta_3 ethnicity_i + \epsilon \]

# Omit negative values
  wages <- wages %>%
    filter(experience > 0)

# Plot wage on education
ggplot(data = wages, 
       mapping = aes(x = education, y = wage)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_x_continuous(breaks = seq(min(wages$education), max(wages$education), by = 2))+

  labs(x = "Education, years", y = "Wage, $ per week")

## `geom_smooth()` using formula = 'y ~ x'

  #plot wage on experience
  ggplot(data = wages, 
         mapping = aes(x = experience, y = wage)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    scale_x_continuous(breaks = seq(min(wages$experience), max(wages$experience), by = 5))+
  
    labs(x = "Experience, years", y = "Wage, $ per week")

## `geom_smooth()` using formula = 'y ~ x'

The scatter plot above for education does show that there may be heteroskedasticity in the model and some potential outliers. Individuals with less than 9 years of education all seem to have much less variation in wages than those who went to high school (9th grade and above). The plot for experience, however, shows no clear evidence of heteroskedasticity or even a relationship between experience and wage.

# Run main regression
  main_reg <- lm(wage ~ education + experience + ethnicity, data = wages)
  
  #display results
  stargazer(main_reg, 
            title = "Main Regression", align = TRUE, type = "text", header = FALSE)

## 
## Main Regression
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                                 wage            
## ------------------------------------------------
## education                    59.901***          
##                               (0.895)           
##                                                 
## experience                    9.184***          
##                               (0.206)           
##                                                 
## ethnicity                    135.065***         
##                               (9.288)           
##                                                 
## Constant                    -457.762***         
##                               (15.511)          
##                                                 
## ------------------------------------------------
## Observations                   26,895           
## R2                             0.172            
## Adjusted R2                    0.172            
## Residual Std. Error     412.765 (df = 26891)    
## F Statistic         1,856.822*** (df = 3; 26891)
## ================================================
## Note:                *p<0.1; **p<0.05; ***p<0.01

  #Show res plots
  par(mfrow = c(2, 2))
  plot(main_reg)

The sheer number of observations makes it difficult to assess heteroskedasticity by eyeballing the scale-location plot. Thus, a test is necessary.

3.4 White Test

  # Conduct White test for heteroscedasticity
  skedastic_package_white  <- white(mainlm =  main_reg, 
                                    interactions = TRUE
                                    )         
  
  # View the test results
  skedastic_package_white

## # A tibble: 1 × 5
##   statistic   p.value parameter method       alternative
##       <dbl>     <dbl>     <dbl> <chr>        <chr>      
## 1      34.8 0.0000645         9 White's Test greater

With a test-stat of 34.81 and a p-value of 6.45e-05, we can reject the null hypothesis and conclude that there is evidence for heteroskedasticity.

3.5 Manually Testing

3.5.1 Create Residuals squared

# Extract model residuals from objects returned by modeling functions.
  wages$residuals <- resid(object = main_reg)
    summary(wages$residuals) # mean should be zero

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1116.6  -216.1   -49.9     0.0   153.0 18113.9

# Square the residuals
  wages$squared_residuals <- (wages$residuals)^2 
    summary(wages$squared_residuals) # should not have any negative values

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##         0      8087     36501    170350    111543 328114429

    wages$residuals <- NULL  # remove residuals from table 

3.5.2 Create the X variables

#Run white aux reg
  white_auxillary_reg <- lm(formula = squared_residuals ~ 
                              #base variables
                              education + experience + ethnicity  + 
                              
                              #squared values of base variables
                              I(education^2)  + I(experience^2)   + I(ethnicity^2) +
                              
                              # Interaction terms
                              education:experience + experience:ethnicity +education:ethnicity,
                               data = wages
                     ) 
  
  white_auxillary_reg_summary <- summary(white_auxillary_reg)
  white_auxillary_reg_summary

## 
## Call:
## lm(formula = squared_residuals ~ education + experience + ethnicity + 
##     I(education^2) + I(experience^2) + I(ethnicity^2) + education:experience + 
##     experience:ethnicity + education:ethnicity, data = wages)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
##   -595816   -155501    -82383    -25136 327916662 
## 
## Coefficients: (1 not defined because of singularities)
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          369875.0   446647.7   0.828   0.4076  
## education            -58631.4    43354.1  -1.352   0.1763  
## experience           -10563.3    10385.0  -1.017   0.3091  
## ethnicity            -21420.7   352712.2  -0.061   0.9516  
## I(education^2)         2906.4     1222.1   2.378   0.0174 *
## I(experience^2)         176.5      105.4   1.674   0.0942 .
## I(ethnicity^2)             NA         NA      NA       NA  
## education:experience    353.1      507.4   0.696   0.4865  
## experience:ethnicity   3128.5     5059.9   0.618   0.5364  
## education:ethnicity    2946.0    23892.9   0.123   0.9019  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2684000 on 26886 degrees of freedom
## Multiple R-squared:  0.001294,   Adjusted R-squared:  0.000997 
## F-statistic: 4.355 on 8 and 26886 DF,  p-value: 2.883e-05