5 Hands-on Tutorial
Let’s start by load these packages:
library(plm)
library(lmtest)
library(texreg)Clear the environment
rm(list = ls())5.1 More Guns, Less Crimes
Download the guns dataset used by Stock and Watson.
Gun rights advocate John Lott argues in his book More Guns, Less Crimes that crime rates in the United States decrease when gun ownership restrictions are relaxed. The data used in Lott’s research compares violent crimes, robberies, and murders across 50 states to determine whether the so called “shall” laws that remove discretion from license granting authorities actually decrease crime rates. So far 41 states have passed these “shall” laws where a person applying for a licence to carry a concealed weapon doesn’t have to provide justification or “good cause” for requiring a concealed weapon permit.
Let’s load the dataset used by Lott and see if we can test the arguments made by gun rights advocates.
guns <- read.csv("guns.csv")The variables we’re interested in are described below.
| Variable | Definition |
|---|---|
| mur | Murder rate (incidents per 100,000) |
| shall | = 1 if the state has a shall-carry law in effect in that year = 0 otherwise |
| incarc_rate | Incarceration rate in the state in the previous year (sentenced prisoners per 100,000 residents; value for the previous year) |
| pm1029 | Percent of state population that is male, ages 10 to 29 |
| stateid | ID number of states (Alabama = 1, Alaska = 2, etc.) |
| year | Year (1977-1999) |
We will focus on murder rates in this example but you could try the same with variables measuring violent crimes or robberies as well.
Let’s create a factor variable representing whether a state has passed “shall” law or not. The variable already exists as 0 or 1 but we want to convert it to a factor for our analysis.
guns$shall <- factor(guns$shall, levels = c(0, 1), labels =c("NO", "YES"))5.2 Fixed Effects
Let’s estimate a fixed effect model on panel data using the plm() function with shall, incarc_rate, and pm1029 as the independent variables.
state_effects <- plm(mur ~ shall + incarc_rate + pm1029,
data = guns,
index = c("stateid", "year"),
model = "within",
effect = "individual")
summary(state_effects)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = mur ~ shall + incarc_rate + pm1029, data = guns,
## effect = "individual", model = "within", index = c("stateid",
## "year"))
##
## Balanced Panel: n=51, T=23, N=1173
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -21.102428 -0.958945 0.016047 1.082008 29.031961
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## shallYES -1.4513886 0.3154300 -4.6013 4.678e-06 ***
## incarc_rate 0.0174551 0.0011261 15.4998 < 2.2e-16 ***
## pm1029 0.9582993 0.0859610 11.1481 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 12016
## Residual Sum of Squares: 9800
## R-Squared: 0.18444
## Adj. R-Squared: 0.14581
## F-statistic: 84.3526 on 3 and 1119 DF, p-value: < 2.22e-16The state_effects model shows that all three of our independent variables are statistically significant, with shall decreasing murder rates by 1.45 incidents per 100000 members of the population. The effects of incarceration rate and percentage of male population between 10 and 29 years old are also statistically significant.
Before drawing any conclusions let’s make sure whether there are any state effects in our model using plmtest().
plmtest(state_effects, effect="individual")##
## Lagrange Multiplier Test - (Honda) for balanced panels
##
## data: mur ~ shall + incarc_rate + pm1029
## normal = 47.242, p-value < 2.2e-16
## alternative hypothesis: significant effectsThe p-value suggests the presence of state effects. In addition to state fixed effects, a number of factors could affect the murder rate that are not specific to an individual state. We can model these time fixed effects using the effect = "time" argument in plm().
time_effects <- plm(mur ~ shall + incarc_rate + pm1029,
data = guns,
index = c("stateid", "year"),
model = "within",
effect = "time")
summary(time_effects)## Oneway (time) effect Within Model
##
## Call:
## plm(formula = mur ~ shall + incarc_rate + pm1029, data = guns,
## effect = "time", model = "within", index = c("stateid", "year"))
##
## Balanced Panel: n=51, T=23, N=1173
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -21.68350 -2.04596 -0.31955 1.76758 35.20084
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## shallYES -0.2521605 0.3032163 -0.8316 0.4058
## incarc_rate 0.0412157 0.0007704 53.4993 <2e-16 ***
## pm1029 0.2148597 0.1407613 1.5264 0.1272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 65571
## Residual Sum of Squares: 18141
## R-Squared: 0.72334
## Adj. R-Squared: 0.71731
## F-statistic: 999.623 on 3 and 1147 DF, p-value: < 2.22e-16The incarc_rate variable is the only statistically significant variable in the time fixed effects model.
Now let’s run plmtest on the time_effects model to verify if time fixed effects are indeed present in the model.
plmtest(time_effects, effect="time")##
## Lagrange Multiplier Test - time effects (Honda) for balanced
## panels
##
## data: mur ~ shall + incarc_rate + pm1029
## normal = 16.104, p-value < 2.2e-16
## alternative hypothesis: significant effectsThe p-value tells us that we can reject the null hypothesis so we know that there are time fixed effects present in our model.
We already confirmed the presence of state fixed effects in the first model we estimated. Now, in order to control for both state AND time fixed effects, we need to estimate a model using the effect = "twoways" argument.
twoway_effects <- plm(mur ~ shall + incarc_rate + pm1029,
data = guns,
index = c("stateid", "year"),
model = "within",
effect = "twoways")
summary(twoway_effects)## Twoways effects Within Model
##
## Call:
## plm(formula = mur ~ shall + incarc_rate + pm1029, data = guns,
## effect = "twoways", model = "within", index = c("stateid",
## "year"))
##
## Balanced Panel: n=51, T=23, N=1173
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -19.2097691 -0.9748749 -0.0069663 1.0119176 27.1354552
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## shallYES -0.5640474 0.3325054 -1.6964 0.0901023 .
## incarc_rate 0.0209756 0.0011252 18.6411 < 2.2e-16 ***
## pm1029 0.7326357 0.2189770 3.3457 0.0008485 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 11263
## Residual Sum of Squares: 8519.4
## R-Squared: 0.24357
## Adj. R-Squared: 0.19186
## F-statistic: 117.746 on 3 and 1097 DF, p-value: < 2.22e-16In a twoway fixed effects model shall is no longer significant and the effect of male population between 10 and 29 years old has decreased from 0.95 to 0.73 incidents per 100,000 population.
The results of all three models are shown below.
screenreg(list(state_effects, time_effects, twoway_effects),
custom.model.names = c("State Fixed Effects",
"Time Fixed Effects",
"Twoway Fixed Effects"))##
## ==========================================================================
## State Fixed Effects Time Fixed Effects Twoway Fixed Effects
## --------------------------------------------------------------------------
## shallYES -1.45 *** -0.25 -0.56
## (0.32) (0.30) (0.33)
## incarc_rate 0.02 *** 0.04 *** 0.02 ***
## (0.00) (0.00) (0.00)
## pm1029 0.96 *** 0.21 0.73 ***
## (0.09) (0.14) (0.22)
## --------------------------------------------------------------------------
## R^2 0.18 0.72 0.24
## Adj. R^2 0.15 0.72 0.19
## Num. obs. 1173 1173 1173
## ==========================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.055.3 Serial Correlation
For time series data we need to address the potential for serial correlation in the error term. We will test for serial correlation with Breusch-Godfrey test using pbgtest() and provide solutions for correcting it if necessary.
pbgtest(twoway_effects)##
## Breusch-Godfrey/Wooldridge test for serial correlation in panel
## models
##
## data: mur ~ shall + incarc_rate + pm1029
## chisq = 765.16, df = 23, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errorsThe null hypothesis for the Breusch-Godfrey test is that there is no serial correlation. The p-value from the test tells us that we can reject the null hypothesis and confirms the presence of serial correlation in our error term.
We can correct for serial correlation using coeftest() similar to how we corrected for heteroskedastic errors. We’ll use the vcovHC() function for obtaining a heteroskedasticity-consistent covariance matrix, but since we’re interested in correcting for autocorrelation as well, we will specify method = "arellano" which corrects for both heteroskedasticity and autocorrelation.
twoway_effects_hac <- coeftest(twoway_effects,
vcov = vcovHC(twoway_effects,
method = "arellano",
type = "HC3"))
screenreg(list(twoway_effects, twoway_effects_hac),
custom.model.names = c("Twoway Fixed Effects",
"Twoway Fixed Effects (HAC)"))##
## =============================================================
## Twoway Fixed Effects Twoway Fixed Effects (HAC)
## -------------------------------------------------------------
## shallYES -0.56 -0.56
## (0.33) (0.48)
## incarc_rate 0.02 *** 0.02 *
## (0.00) (0.01)
## pm1029 0.73 *** 0.73
## (0.22) (0.54)
## -------------------------------------------------------------
## R^2 0.24
## Adj. R^2 0.19
## Num. obs. 1173
## =============================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05We can see that with heteroskedasticity and autocorrelation consistent (HAC) standard errors, the percent of male population (10 - 29 yr old) is no longer a significant predictor in our model.
5.4 Cross-sectional Dependence
If a federal law imposed restrictions on gun ownership or licensing requirements then the changes would likely affect all 50 states. This is an example of Cross Sectional Dependence and not accounted for in a fixed effect model. Other scenarios could also trigger cross sectional dependence that we should take into consideration. For example, security policies and law enforcement efforts might change after an extraordinary event (think of mass shootings or terrorist attacks) thus influencing law enforcement practices in all states. We can check for cross sectional dependence using the Pesaran cross sectional dependence test or pcdtest().
pcdtest(twoway_effects)##
## Pesaran CD test for cross-sectional dependence in panels
##
## data: mur ~ shall + incarc_rate + pm1029
## z = 3.9121, p-value = 9.148e-05
## alternative hypothesis: cross-sectional dependenceAs we’ve seen with other tests, the null hypothesis is that there is no cross sectional dependence. The p-value, however tells that there is indeed cross-sectional dependence and we need to correct it. There are two general approaches to correcting for cross sectional dependence.
Beck and Katz (1995) method or Panel Corrected Standard Errors (PCSE): We can obtain Panel Corrected Standard Errors (PCSE) by first obtaining a robust variance-covariance matrix for panel models with the Beck and Katz (1995) method using the vcovBK() and passing it to the familiar coeftest() function.
twoway_effects_pcse <- coeftest(twoway_effects,
vcov = vcovBK(twoway_effects,
type="HC3",
cluster = "group"))
twoway_effects_pcse##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## shallYES -0.5640474 0.7662556 -0.7361 0.4618
## incarc_rate 0.0209756 0.0028249 7.4253 2.254e-13 ***
## pm1029 0.7326357 0.5118496 1.4313 0.1526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The results from PCSE are sensitive to the ratio between the number of time periods in the dataset (T) and the total number of observations (N). When we’re dealing with large datasets (i.e. the T/N ratio is small), we use the Driscoll and Kraay method instead.
Driscoll and Kraay (1998) (SCC): The cross-sectional and serial correlation (SCC) method by Driscoll and Kraay addresses the limitations of Beck and Katz’s PCSE method is therefore preferred for obtaining heteroskedasticity and autocorrelation consistent errors that are also robust to cross-sectional dependence. We can get SCC corrected covariance matrix using the vcovSCC() function.
twoway_effects_scc <- coeftest(twoway_effects,
vcov = vcovSCC(twoway_effects,
type="HC3",
cluster = "group"))
twoway_effects_scc##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## shallYES -0.564047 0.542698 -1.0393 0.29888
## incarc_rate 0.020976 0.010321 2.0324 0.04236 *
## pm1029 0.732636 0.551066 1.3295 0.18396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1screenreg(list(state_effects,
time_effects,
twoway_effects,
twoway_effects_pcse,
twoway_effects_scc),
custom.model.names = c("State Effects",
"Time Effects",
"Twoway Fixed Effects",
"PCSE",
"SCC"))##
## ==================================================================================
## State Effects Time Effects Twoway Fixed Effects PCSE SCC
## ----------------------------------------------------------------------------------
## shallYES -1.45 *** -0.25 -0.56 -0.56 -0.56
## (0.32) (0.30) (0.33) (0.77) (0.54)
## incarc_rate 0.02 *** 0.04 *** 0.02 *** 0.02 *** 0.02 *
## (0.00) (0.00) (0.00) (0.00) (0.01)
## pm1029 0.96 *** 0.21 0.73 *** 0.73 0.73
## (0.09) (0.14) (0.22) (0.51) (0.55)
## ----------------------------------------------------------------------------------
## R^2 0.18 0.72 0.24
## Adj. R^2 0.15 0.72 0.19
## Num. obs. 1173 1173 1173
## ==================================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05