6.4 Diagnostic checking in RStudio

Diagnostic checking in R Studio requires installation of additional packages

TABLE 6.3: RStudio packages and commands that support diagnostic checking
Package	Command	Description
`lmtest`	`bptest()`	Breusch-Pagan test for heteroscedasticity
`lmtest`	`bgtest()`	Breusch-Godfrey test for autocorrelation
`lmtest`	`coeftest()`	significance testing by robust standard errors
`car`	`hccm()`	heteroscedasticity corrected covariance matrix
`car`	`vif()`	variance inflation factors
`tseries`	`jarque.bera.test()`	Jarque-Bera test for normality
`sandwich`	`vcovHAC()`	heteroscedasticity and autocorrelation corrected covariance matrix

Exercise 31. Using sample data from a newdata object estimate a multivariate model by OLS method: \[gdp_i=\beta_0+\beta_1population_i+\beta_2unemployment_i+u_i\] Perform diagnostic checking of the multivariate model (significance level \(\alpha=0.05\)). Compute White’s robust standard errors. Also, apply WLS method (take the inverse values of the population as the weights \(w_{ii}\)). Summarize all results together in a single table using modelsummary() command.

Solution

Copy the code lines below to the clipboard, paste them into an R Script file opened in RStudio, and run them. In this example a multivariate model passes all diagnostic tests, except for BP test, which indicates the presence of heteroscedasticity problem (th null hypothesis is rejected in both versions of BP test at significance level \(5\%\)). Two solutions of heteroscedasticity problem are considered: a) using robust standard errors which are heteroscedasticity corrected (so called White’s robust standard errors), and b) re-estimating the model using WLS method (assuming that the source of heteroscedasticity is in population).

# OLS estimation of multivariate model by command lm()
multivariate=lm(gdp~population+unemployment,data=newdata) 

# Reporting only estimates, their standard errors, t-statistics and p-values
library(lmtest) # loading already installed package from the library that supports coeftest()
coeftest(multivariate)

# Checking for multicolinearity by VIF (Variance Inflation Factor)
library(car) # loading already installed package from the library that supports vif()
vif(multivariate)

# BP test using ordinary residuals
bptest(multivariate,studentize=FALSE)
# BP test (White's version) using squared RHS variables and interaction terms

bptest(multivariate,~population+unemployment+I(population^2)+I(unemployment^2)+population:unemployment,studentize=FALSE,data=newdata)

# First order BG test
bgtest(multivariate) # argument order=1 is omitted as it is default
# Second order BG test
bgtest(multivariate,order=2)

# Normality test of residuals
install.packages("tseries") # this package should be installed (for the first time) 
library(tseries) # loading the same package from the library
jarque.bera.test(resid(multivariate)) # normality test based on skewness and kurtosis of residuals

# Reporting the matrix "Gamma-hat" and classical (non-robust) standard errors 
vcov(multivariate) # covariance matrix "Gamma-hat"
sqrt(diag(vcov(multivariate))) # square root of it's diagonal elements are standard errors

# Reporting matrix "Gamma-tilde" and robust standard errors (corrected for heteroscedasticity)
hccm(multivariate,type="hc0") # argument type="hc0" sets White covariance matrix 
sqrt(diag(hccm(multivariate,type="hc0")))
# Significance testing using robust standard errors
coeftest(multivariate,vcov=hccm(multivariate,type="hc0"))

# WLS estimation method
multivariate.wls=lm(gdp~population+unemployment,weights=1/population,data=newdata)
coeftest(multivariate.wls)

# Summarizing all results in a single table
install.packages("sandwich") # installing "sandwich" package which arguments of the vcov() command can be passed directly into modelsummary()
library(sandwich) # loading installed package from the library
library(modelsummary) # loading installed package from the library
modelsummary(list("OLS with classical standard errors"=multivariate,"OLS with robust standard errors"=multivariate,"WLS with classical standard errors"=multivariate.wls),vcov=c("classical", "HC0", "classical"), stars=TRUE,fmt=4)