# Clear the workspace
rm(list = ls()) # Clear environment
gc() # Clear unused memory
cat("\f") # Clear the console
dev.off # Clear the charts#List all required packages
packages <- c("plm",
"tidyverse",
"dplyr",
"ggplot2",
"kableExtra",
"scales",
"psych",
"reshape2",
"corrplot",
"gridExtra",
"car",
"stargazer"
)
for (i in 1:length(packages)) {
if (!packages[i] %in% rownames(installed.packages())) {
install.packages(packages[i]
, repos = "http://cran.rstudio.com/"
, dependencies = TRUE
)
}
library(packages[i], character.only = TRUE)
}
rm(packages)# Set working directory (replace with your actual path)
setwd("~/Desktop/")
# Load Guns Dataset
guns <- read_csv("Guns-1.csv")
# Display sample of data
kable(head(guns), caption = "Guns Dataset")| year | violent | murder | robbery | prisoners | afam | cauc | male | population | income | density | state | law |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1977 | 414.4 | 14.2 | 96.8 | 83 | 8.384873 | 55.12291 | 18.17441 | 3.780403 | 9563.148 | 0.0745524 | Alabama | no |
| 1978 | 419.1 | 13.3 | 99.1 | 94 | 8.352101 | 55.14367 | 17.99408 | 3.831838 | 9932.000 | 0.0755667 | Alabama | no |
| 1979 | 413.3 | 13.2 | 109.5 | 144 | 8.329575 | 55.13586 | 17.83934 | 3.866248 | 9877.028 | 0.0762453 | Alabama | no |
| 1980 | 448.5 | 13.2 | 132.1 | 141 | 8.408386 | 54.91259 | 17.73420 | 3.900368 | 9541.428 | 0.0768288 | Alabama | no |
| 1981 | 470.5 | 11.9 | 126.5 | 149 | 8.483435 | 54.92513 | 17.67372 | 3.918531 | 9548.351 | 0.0771866 | Alabama | no |
| 1982 | 447.7 | 10.6 | 112.0 | 183 | 8.514000 | 54.89621 | 17.51052 | 3.925229 | 9478.919 | 0.0773185 | Alabama | no |
# Converting to panel data format
guns_panel <- pdata.frame(guns, index = c("state", "year"))
# Check if the dataset is balanced
is_balanced <- is.pbalanced(guns_panel)
print(is_balanced) # TRUE## [1] TRUE# Converting to panel data format
guns_panel <- pdata.frame(guns, index = c("state", "year"))
# Check if the dataset is balanced
is_balanced <- is.pbalanced(guns_panel)
print(is_balanced) # TRUE## [1] TRUE\[ \text{Violent Crime Rate}_i = \beta_0 + \beta_1 \cdot \text{Law}_i + \beta_2 \cdot \text{Prisoners}_i + \beta_3 \cdot \text{Income}_i + \beta_4 \cdot \text{Population}_i + \beta_5 \cdot \text{Male}_i + \beta_6 \cdot \text{Afam}_i + \beta_7 \cdot \text{Cauc}_i + \epsilon_i \]
# OLS model
ols_model <- lm(violent ~ law + prisoners + income + population + male + afam + cauc, data = guns)
# Summary of the OLS model
stargazer(ols_model, type = 'text')##
## ===============================================
## Dependent variable:
## ---------------------------
## violent
## -----------------------------------------------
## lawyes -87.029***
## (15.067)
##
## prisoners 1.102***
## (0.046)
##
## income 0.021***
## (0.003)
##
## population 13.782***
## (1.144)
##
## male 32.176***
## (4.764)
##
## afam -17.653**
## (7.477)
##
## cauc -15.237***
## (3.701)
##
## Constant 456.713*
## (248.122)
##
## -----------------------------------------------
## Observations 1,173
## R2 0.652
## Adjusted R2 0.649
## Residual Std. Error 197.926 (df = 1165)
## F Statistic 311.140*** (df = 7; 1165)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Direction:
Most variables align with expectations. For instance, as the number of prisoners or male population increases, violent crime rates also increase, which makes sense given the demographic association with violent crime.
Surprisingly, income has a small positive effect on violent crime. Generally, we expect rising income to lower crime rates, but it’s plausible that wealthier areas could attract more targeted violent crimes, like robbery.
LawYes (Shall-carry law) significantly reduces violent crime, as expected, suggesting that these laws may act as a deterrent.
Magnitude:
The most notable effects are from LawYes and Male population.
The LawYes variable indicates that states with shall-carry laws see a reduction of ~87 violent crimes, which is a strong effect.
Male population has a notable increase of ~32 units for every percentage point rise in the male share of the population, consistent with prior literature linking males to higher violent crime rates.
Statistical Significance:
Omitted Variable Bias:
\[ \text{Violent Crime Rate}_i = \beta_0 + \beta_1 \cdot \text{Law}_i + \beta_2 \cdot \text{Prisoners}_i + \beta_3 \cdot \text{Income}_i + \beta_4 \cdot \text{Population}_i \\ + \beta_5 \cdot \text{Male}_i + \beta_6 \cdot \text{Afam}_i + \beta_7 \cdot \text{Cauc}_i + \sum_{s=1}^{49} \delta_s \cdot \text{State}_s + \epsilon_i \]
ols_model_dummy <- lm(violent ~ law + prisoners + income + population + male + afam + cauc + state, data = guns)
# Summary
stargazer(ols_model_dummy, type = "text")##
## =====================================================
## Dependent variable:
## ---------------------------
## violent
## -----------------------------------------------------
## lawyes -21.991*
## (11.548)
##
## prisoners 0.212***
## (0.043)
##
## income -0.004
## (0.004)
##
## population 10.143*
## (5.319)
##
## male -21.564***
## (3.862)
##
## afam -3.949
## (9.364)
##
## cauc 9.644***
## (3.113)
##
## stateAlaska 91.021**
## (44.728)
##
## stateArizona -45.728
## (49.454)
##
## stateArkansas -158.466***
## (38.798)
##
## stateCalifornia 53.585
## (141.700)
##
## stateColorado -201.264***
## (64.347)
##
## stateConnecticut -233.689***
## (63.545)
##
## stateDelaware -12.543
## (39.159)
##
## stateDistrict of Columbia 1,750.092***
## (142.044)
##
## stateFlorida 292.074***
## (67.521)
##
## stateGeorgia 50.772
## (32.544)
##
## stateHawaii 142.710
## (136.247)
##
## stateIdaho -397.551***
## (67.286)
##
## stateIllinois 161.152**
## (64.941)
##
## stateIndiana -266.379***
## (62.828)
##
## stateIowa -437.006***
## (71.579)
##
## stateKansas -254.441***
## (57.293)
##
## stateKentucky -344.681***
## (58.422)
##
## stateLouisiana 275.271***
## (34.642)
##
## stateMaine -527.630***
## (75.152)
##
## stateMaryland 307.610***
## (34.177)
##
## stateMassachusetts -66.184
## (71.009)
##
## stateMichigan 6.534
## (58.301)
##
## stateMinnesota -395.492***
## (72.544)
##
## stateMississippi -73.221*
## (44.179)
##
## stateMissouri -83.513
## (52.870)
##
## stateMontana -467.233***
## (59.721)
##
## stateNebraska -347.897***
## (64.652)
##
## stateNevada 91.108*
## (48.346)
##
## stateNew Hampshire -541.407***
## (78.063)
##
## stateNew Jersey -95.270*
## (53.873)
##
## stateNew Mexico 140.667***
## (44.708)
##
## stateNew York 212.864**
## (89.014)
##
## stateNorth Carolina -56.672*
## (33.768)
##
## stateNorth Dakota -527.791***
## (63.772)
##
## stateOhio -290.137***
## (72.372)
##
## stateOklahoma -135.932***
## (38.762)
##
## stateOregon -190.559***
## (65.628)
##
## statePennsylvania -346.114***
## (80.803)
##
## stateRhode Island -276.096***
## (63.443)
##
## stateSouth Carolina 287.200***
## (35.105)
##
## stateSouth Dakota -444.237***
## (56.117)
##
## stateTennessee -43.990
## (42.494)
##
## stateTexas -146.706*
## (88.667)
##
## stateUtah -306.473***
## (63.207)
##
## stateVermont -537.282***
## (75.567)
##
## stateVirginia -269.442***
## (37.594)
##
## stateWashington -199.783***
## (60.850)
##
## stateWest Virginia -494.089***
## (68.843)
##
## stateWisconsin -454.844***
## (66.854)
##
## stateWyoming -374.911***
## (67.449)
##
## Constant 342.495
## (237.027)
##
## -----------------------------------------------------
## Observations 1,173
## R2 0.917
## Adjusted R2 0.913
## Residual Std. Error 98.856 (df = 1115)
## F Statistic 215.544*** (df = 57; 1115)
## =====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01\[ \text{Violent Crime Rate}_i - \overline{\text{Violent Crime Rate}}_s = \beta_1 \cdot (\text{Law}_i - \overline{\text{Law}}_s) + \beta_2 \cdot (\text{Prisoners}_i - \overline{\text{Prisoners}}_s) + \dots + \epsilon_i \]
# Convert 'law' column to numeric
guns <- guns %>% mutate(law_numeric = ifelse(law == "yes", 1, 0))
# Demean the data by subtracting the state-level means
guns_demeaned <- with(guns, data.frame(
violent = violent - ave(violent, state),
lawyes = law_numeric - ave(law_numeric, state),
prisoners = prisoners - ave(prisoners, state),
income = income - ave(income, state),
population = population - ave(population, state),
male = male - ave(male, state),
afam = afam - ave(afam, state),
cauc = cauc - ave(cauc, state)
))
# OLS regression on demeaned data (without intercept)
ols_model_demeaned <- lm(violent ~ lawyes + prisoners + income + population + male + afam + cauc - 1, data = guns_demeaned)
# Summary of the demeaned fixed effects model
stargazer(ols_model_demeaned, type = "text")##
## ===============================================
## Dependent variable:
## ---------------------------
## violent
## -----------------------------------------------
## lawyes -21.991*
## (11.292)
##
## prisoners 0.212***
## (0.042)
##
## income -0.004
## (0.004)
##
## population 10.143*
## (5.202)
##
## male -21.564***
## (3.777)
##
## afam -3.949
## (9.157)
##
## cauc 9.644***
## (3.044)
##
## -----------------------------------------------
## Observations 1,173
## R2 0.197
## Adjusted R2 0.193
## Residual Std. Error 96.670 (df = 1166)
## F Statistic 40.970*** (df = 7; 1166)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01This method automatically de-means the data to remove time-invariant state effects:
\[ \text{Violent Crime Rate}_{it} = \beta_1 \cdot \text{Law}_{it} + \beta_2 \cdot \text{Prisoners}_{it} + \beta_3 \cdot \text{Income}_{it} + \beta_4 \cdot \text{Population}_{it} + \dots + \mu_i + \lambda_t + \epsilon_{it} \]
ols_model_within <- plm(violent ~ law + prisoners + income + population + male + afam + cauc,
data = guns,
index = c("state", "year"),
model = "within")
# Summary of the within fixed effects model
stargazer(ols_model_within, type = "text")##
## ========================================
## Dependent variable:
## ---------------------------
## violent
## ----------------------------------------
## lawyes -21.991*
## (11.548)
##
## prisoners 0.212***
## (0.043)
##
## income -0.004
## (0.004)
##
## population 10.143*
## (5.319)
##
## male -21.564***
## (3.862)
##
## afam -3.949
## (9.364)
##
## cauc 9.644***
## (3.113)
##
## ----------------------------------------
## Observations 1,173
## R2 0.197
## Adjusted R2 0.156
## F Statistic 39.178*** (df = 7; 1115)
## ========================================
## Note: *p<0.1; **p<0.05; ***p<0.01In the two-way fixed effects model, we control for both state and year effects:
\[ \text{Violent Crime Rate}_{it} = \beta_0 + \beta_1 \cdot \text{Law}_{it} + \beta_2 \cdot \text{Prisoners}_{it} + \beta_3 \cdot \text{Income}_{it} + \dots + \sum_{s=1}^{49} \delta_s \cdot \text{State}_s + \sum_{t=1}^{T} \gamma_t \cdot \text{Year}_t + \epsilon_{it} \]
# OLS model with two-way fixed effects (state and year)
ols_model_SY <- lm(violent ~ law + prisoners + income + population + male + afam + cauc +
as.factor(state) + as.factor(year), data = guns)
# Summary
stargazer(ols_model_SY, type = "text")##
## ================================================================
## Dependent variable:
## ---------------------------
## violent
## ----------------------------------------------------------------
## lawyes -1.707
## (10.410)
##
## prisoners 0.267***
## (0.041)
##
## income 0.008**
## (0.004)
##
## population 1.205
## (4.759)
##
## male 45.848***
## (9.459)
##
## afam -42.638***
## (13.406)
##
## cauc -10.033**
## (4.753)
##
## as.factor(state)Alaska -49.491
## (41.274)
##
## as.factor(state)Arizona -78.796*
## (44.885)
##
## as.factor(state)Arkansas -157.102***
## (36.423)
##
## as.factor(state)California 210.636*
## (126.428)
##
## as.factor(state)Colorado -186.328***
## (56.221)
##
## as.factor(state)Connecticut -222.116***
## (56.713)
##
## as.factor(state)Delaware -42.178
## (34.460)
##
## as.factor(state)District of Columbia 1,618.748***
## (141.082)
##
## as.factor(state)Florida 392.468***
## (60.069)
##
## as.factor(state)Georgia 39.030
## (30.504)
##
## as.factor(state)Hawaii 74.334
## (130.172)
##
## as.factor(state)Idaho -452.400***
## (66.158)
##
## as.factor(state)Illinois 207.578***
## (57.163)
##
## as.factor(state)Indiana -258.949***
## (54.959)
##
## as.factor(state)Iowa -438.109***
## (66.451)
##
## as.factor(state)Kansas -286.839***
## (52.867)
##
## as.factor(state)Kentucky -323.605***
## (52.898)
##
## as.factor(state)Louisiana 206.075***
## (30.152)
##
## as.factor(state)Maine -491.267***
## (67.657)
##
## as.factor(state)Maryland 331.528***
## (39.500)
##
## as.factor(state)Massachusetts -36.225
## (62.049)
##
## as.factor(state)Michigan 44.623
## (50.742)
##
## as.factor(state)Minnesota -391.470***
## (64.809)
##
## as.factor(state)Mississippi -161.011***
## (39.003)
##
## as.factor(state)Missouri -57.173
## (46.516)
##
## as.factor(state)Montana -451.407***
## (54.720)
##
## as.factor(state)Nebraska -374.011***
## (60.919)
##
## as.factor(state)Nevada 101.841**
## (42.874)
##
## as.factor(state)New Hampshire -530.460***
## (68.700)
##
## as.factor(state)New Jersey -23.346
## (51.471)
##
## as.factor(state)New Mexico 116.858***
## (42.772)
##
## as.factor(state)New York 348.503***
## (82.516)
##
## as.factor(state)North Carolina -43.834
## (30.043)
##
## as.factor(state)North Dakota -623.353***
## (65.111)
##
## as.factor(state)Ohio -206.816***
## (62.890)
##
## as.factor(state)Oklahoma -148.461***
## (34.227)
##
## as.factor(state)Oregon -145.486**
## (57.762)
##
## as.factor(state)Pennsylvania -232.074***
## (70.683)
##
## as.factor(state)Rhode Island -281.459***
## (57.356)
##
## as.factor(state)South Carolina 246.525***
## (30.898)
##
## as.factor(state)South Dakota -513.456***
## (55.795)
##
## as.factor(state)Tennessee -1.456
## (36.843)
##
## as.factor(state)Texas -93.520
## (76.833)
##
## as.factor(state)Utah -498.669***
## (73.491)
##
## as.factor(state)Vermont -536.024***
## (68.062)
##
## as.factor(state)Virginia -254.251***
## (35.831)
##
## as.factor(state)Washington -185.820***
## (52.999)
##
## as.factor(state)West Virginia -429.218***
## (62.588)
##
## as.factor(state)Wisconsin -450.417***
## (59.985)
##
## as.factor(state)Wyoming -434.433***
## (62.624)
##
## as.factor(year)1978 22.764
## (17.084)
##
## as.factor(year)1979 72.057***
## (17.362)
##
## as.factor(year)1980 112.780***
## (17.652)
##
## as.factor(year)1981 123.398***
## (18.276)
##
## as.factor(year)1982 112.761***
## (19.252)
##
## as.factor(year)1983 87.623***
## (20.613)
##
## as.factor(year)1984 92.569***
## (22.722)
##
## as.factor(year)1985 112.612***
## (24.807)
##
## as.factor(year)1986 149.012***
## (27.176)
##
## as.factor(year)1987 146.649***
## (29.632)
##
## as.factor(year)1988 181.475***
## (32.378)
##
## as.factor(year)1989 208.459***
## (34.916)
##
## as.factor(year)1990 289.252***
## (42.716)
##
## as.factor(year)1991 322.070***
## (44.803)
##
## as.factor(year)1992 337.185***
## (47.294)
##
## as.factor(year)1993 349.598***
## (49.130)
##
## as.factor(year)1994 331.764***
## (51.368)
##
## as.factor(year)1995 318.127***
## (53.519)
##
## as.factor(year)1996 281.507***
## (55.575)
##
## as.factor(year)1997 258.594***
## (57.551)
##
## as.factor(year)1998 220.464***
## (59.810)
##
## as.factor(year)1999 186.109***
## (61.879)
##
## Constant 379.268
## (286.651)
##
## ----------------------------------------------------------------
## Observations 1,173
## R2 0.940
## Adjusted R2 0.935
## Residual Std. Error 85.026 (df = 1093)
## F Statistic 215.469*** (df = 79; 1093)
## ================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Generate a summary table with custom column labels
stargazer(ols_model_dummy,
ols_model_demeaned,
ols_model_within,
type = "text",
column.labels = c("Dummy Fixed Effects", "Demeaned Fixed Effects",
"Within Fixed Effects (PLM)"),
title = "Regression Results: Comparison of Fixed Effects Models",
dep.var.labels = "Violent Crime Rate",
covariate.labels = c("Shall Carry Law", "Prisoners", "Income",
"Population", "Male (10-29)", "African American",
"Caucasian"),
no.space = TRUE,
omit = c("state"), # Adjust this if you're facing issues with fixed effects
omit.stat = c("f", "ser", "adj.rsq"), # Optional: omit stats you don't need
notes = "Standard errors in parentheses") # Optional: custom notes##
## Regression Results: Comparison of Fixed Effects Models
## ======================================================================================
## Dependent variable:
## ---------------------------------------------------------------------
## Violent Crime Rate
## OLS panel
## linear
## Dummy Fixed Effects Demeaned Fixed Effects Within Fixed Effects (PLM)
## (1) (2) (3)
## --------------------------------------------------------------------------------------
## Shall Carry Law -21.991* -21.991* -21.991*
## (11.548) (11.292) (11.548)
## Prisoners 0.212*** 0.212*** 0.212***
## (0.043) (0.042) (0.043)
## Income -0.004 -0.004 -0.004
## (0.004) (0.004) (0.004)
## Population 10.143* 10.143* 10.143*
## (5.319) (5.202) (5.319)
## Male (10-29) -21.564*** -21.564*** -21.564***
## (3.862) (3.777) (3.862)
## African American -3.949 -3.949 -3.949
## (9.364) (9.157) (9.364)
## Caucasian 9.644*** 9.644*** 9.644***
## (3.113) (3.044) (3.113)
## Constant 342.495
## (237.027)
## --------------------------------------------------------------------------------------
## Observations 1,173 1,173 1,173
## R2 0.917 0.197 0.197
## ======================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## Standard errors in parentheses