library(psych)library(plm)library(lmtest)library(sandwich)library(ggplot2)library(dplyr)library(GGally)library(car)data <- read.csv('handguns.csv')df <- pdata.frame(data, index = c("state", "year"))reg1 <- plm(log(vio) ~ shall, data = df, model = "pooling")vcov1 <- vcovHC(reg1, method = "arellano", type = "HC1", cluster = "group")coeftest(reg1, vcov = vcov1)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.134919 0.078282 78.370 < 2.2e-16 ***
shall -0.442965 0.155538 -2.848 0.004477 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg2 <- plm(log(vio) ~ shall + incarc_rate + density + avginc + pop + pb1064 + pw1064 + pm1029, data=df, model = "pooling")vcov2 <- vcovHC(reg2, method = "arellano", type = "HC1", cluster = "group")coeftest(reg2, vcov = vcov2)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.98173825 2.14608239 1.3894 0.1649809
shall -0.36838695 0.11286253 -3.2640 0.0011303 **
incarc_rate 0.00161263 0.00059429 2.7135 0.0067551 **
density 0.02668847 0.04109968 0.6494 0.5162340
avginc 0.00120512 0.02385371 0.0505 0.9597156
pop 0.04270983 0.01161835 3.6761 0.0002477 ***
pb1064 0.08085260 0.07071429 1.1434 0.2531201
pw1064 0.03120051 0.03376852 0.9240 0.3557025
pm1029 0.00887088 0.03377483 0.2626 0.7928686
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1The coefficient is \(-0.368\) means that if such parameters as the incarceration rate in the state in previous year, population density, average income, population, percentage of population that is male and aged from 10 to 29, percentage of population that is white and aged from 10 to 64, percentage of population that is black and aged from 10 to 64 are held constant, then having “shall-issue” laws decreases the violence by \(36.8\)%, which is significant in real-life sense.
In the regression (i) the coefficient equals to \(-0.4429\), which means that the estimated decrease is \(44.3\)%. In the regression (ii) the coefficient equals to \(-0.368\), which means that the estimated decrease is \(36.8\)%. Both coefficients are significant at 5% level and thus the difference in the results might be explained by omitted variable biases. In both regressions an estimated decrease of violent crimes due to introduction of “shall-issue” laws is significant. Moreover, the difference between coefficients in these models is \(7.5\) percent points, which is relatively high in real-life sense.
We can think of additional state specific laws that bring more restrictions on having a weapon or on the contrary make it easier to get one. These laws almost don’t vary over years but they are different in each state. An introduction of these additional laws affects crime rate since more or less strict laws affect people’s behavior and at the same time it correlates with shall variable, because it directly affects how many people are able to get a weapon.
Table 1
reg3 <- plm(log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029, data = df, model = "within", effect = "individual")vcov3 <- vcovHC(reg3, method = "arellano", type = "HC1", cluster = "group")coeftest(reg3, vcov = vcov3)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall -4.6142e-02 4.1350e-02 -1.1159 0.264716
incarc_rate -7.1008e-05 2.4788e-04 -0.2865 0.774581
density -1.7229e-01 1.3626e-01 -1.2645 0.206332
avginc -9.2037e-03 1.2837e-02 -0.7170 0.473547
pop 1.1525e-02 1.4084e-02 0.8183 0.413367
pw1064 4.0861e-02 1.3326e-02 3.0663 0.002219 **
pb1064 1.0428e-01 3.2363e-02 3.2222 0.001309 **
pm1029 -5.0273e-02 2.0491e-02 -2.4534 0.014303 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1pFtest(reg3, reg2)
F test for individual effects
data: log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 142.57, df1 = 50, df2 = 1114, p-value < 2.2e-16
alternative hypothesis: significant effectsreg4 <- plm(log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029 + factor(year), data = df, model = "within", effect = "individual")vcov4 <- vcovHC(reg4, method = "arellano", type = "HC1", cluster = "group")coeftest(reg4, vcov = vcov4)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall -2.7994e-02 4.0315e-02 -0.6944 0.4876001
incarc_rate 7.5995e-05 2.0584e-04 0.3692 0.7120521
density -9.1555e-02 1.2264e-01 -0.7465 0.4555075
avginc 9.5865e-04 1.6330e-02 0.0587 0.9531993
pop -4.7544e-03 1.5079e-02 -0.3153 0.7525956
pw1064 9.2501e-03 2.3522e-02 0.3933 0.6942092
pb1064 2.9186e-02 4.9052e-02 0.5950 0.5519639
pm1029 7.3325e-02 5.1956e-02 1.4113 0.1584397
factor(year)78 5.8526e-02 1.5996e-02 3.6588 0.0002656 ***
factor(year)79 1.6395e-01 2.4217e-02 6.7701 2.098e-11 ***
factor(year)80 2.1708e-01 3.3089e-02 6.5604 8.268e-11 ***
factor(year)81 2.1726e-01 3.8809e-02 5.5981 2.741e-08 ***
factor(year)82 1.9463e-01 4.6115e-02 4.2206 2.639e-05 ***
factor(year)83 1.5864e-01 5.8799e-02 2.6981 0.0070812 **
factor(year)84 1.9299e-01 7.6243e-02 2.5312 0.0115054 *
factor(year)85 2.4448e-01 9.1312e-02 2.6774 0.0075314 **
factor(year)86 3.2409e-01 1.0784e-01 3.0052 0.0027148 **
factor(year)87 3.2437e-01 1.2376e-01 2.6210 0.0088885 **
factor(year)88 3.8674e-01 1.3833e-01 2.7958 0.0052679 **
factor(year)89 4.4221e-01 1.5202e-01 2.9089 0.0037006 **
factor(year)90 5.4305e-01 1.9415e-01 2.7970 0.0052481 **
factor(year)91 5.9595e-01 2.0206e-01 2.9494 0.0032516 **
factor(year)92 6.2752e-01 2.1489e-01 2.9202 0.0035702 **
factor(year)93 6.4974e-01 2.2240e-01 2.9215 0.0035555 **
factor(year)94 6.3542e-01 2.3094e-01 2.7514 0.0060321 **
factor(year)95 6.2768e-01 2.3997e-01 2.6157 0.0090280 **
factor(year)96 5.7134e-01 2.5091e-01 2.2771 0.0229726 *
factor(year)97 5.5012e-01 2.5877e-01 2.1259 0.0337393 *
factor(year)98 4.9329e-01 2.7195e-01 1.8139 0.0699627 .
factor(year)99 4.3288e-01 2.8340e-01 1.5275 0.1269357
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg5 <- plm(log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029 + factor(year), data = df, model = "pooling", effect = "individual")pFtest(reg4, reg5)
F test for individual effects
data: log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 181.42, df1 = 50, df2 = 1092, p-value < 2.2e-16
alternative hypothesis: significant effectslinearHypothesis(reg4, matchCoefs(reg4, "year"), vcov = vcov4)Linear hypothesis test
Hypothesis:
factor(year)78 = 0
factor(year)79 = 0
factor(year)80 = 0
factor(year)81 = 0
factor(year)82 = 0
factor(year)83 = 0
factor(year)84 = 0
factor(year)85 = 0
factor(year)86 = 0
factor(year)87 = 0
factor(year)88 = 0
factor(year)89 = 0
factor(year)90 = 0
factor(year)91 = 0
factor(year)92 = 0
factor(year)93 = 0
factor(year)94 = 0
factor(year)95 = 0
factor(year)96 = 0
factor(year)97 = 0
factor(year)98 = 0
factor(year)99 = 0
Model 1: restricted model
Model 2: log(vio) ~ shall + incarc_rate + density + avginc + pop + pw1064 +
pb1064 + pm1029 + factor(year)
Note: Coefficient covariance matrix supplied.
Res.Df Df Chisq Pr(>Chisq)
1 1114
2 1092 22 485.2 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Table 2
reg6 <- plm(log(rob) ~ shall, data = df, model = "pooling")vcov6 <- vcovHC(reg6, method = "arellano", type = "HC1", cluster = "group")coeftest(reg6, vcov = vcov6)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.87305 0.11455 42.5393 < 2.2e-16 ***
shall -0.77332 0.22311 -3.4662 0.0005471 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg7 <- plm(log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029, data = df, model = "pooling")vcov7 <- vcovHC(reg7, method = "arellano", type = "HC1", cluster = "group")coeftest(reg7, vcov = vcov7)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.90413831 3.03262931 0.2981 0.7656520
shall -0.52882023 0.15935939 -3.3184 0.0009333 ***
incarc_rate 0.00100573 0.00063404 1.5862 0.1129556
density 0.09050478 0.04554604 1.9871 0.0471448 *
avginc 0.04073249 0.02789125 1.4604 0.1444489
pop 0.07781765 0.02230700 3.4885 0.0005039 ***
pw1064 0.02752086 0.04458436 0.6173 0.5371733
pb1064 0.10218815 0.08856445 1.1538 0.2488075
pm1029 0.02725647 0.04133195 0.6595 0.5097355
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg8 <- plm(log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029, data = df, model = "within", effect = "individual")vcov8 <- vcovHC(reg8, method = "arellano", type = "HC1", cluster = "group")coeftest(reg8, vcov = vcov8)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall -7.8190e-03 5.4622e-02 -0.1431 0.88620
incarc_rate -7.6342e-05 3.1786e-04 -0.2402 0.81024
density -1.8609e-01 1.6470e-01 -1.1299 0.25877
avginc -1.7519e-02 2.1818e-02 -0.8030 0.42216
pop 1.6333e-02 2.7316e-02 0.5979 0.55000
pw1064 2.7181e-02 1.6272e-02 1.6704 0.09513 .
pb1064 1.1154e-01 5.0650e-02 2.2022 0.02786 *
pm1029 1.1182e-02 2.8811e-02 0.3881 0.69801
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1pFtest(reg8, reg7)
F test for individual effects
data: log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 164.06, df1 = 50, df2 = 1114, p-value < 2.2e-16
alternative hypothesis: significant effectsreg9 <- plm(log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029 + factor(year), data = df, model = "within", effect = "individual")vcov9 <- vcovHC(reg9, method = "arellano", type = "HC1", cluster = "group")coeftest(reg9, vcov = vcov9)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall 2.6830e-02 5.1661e-02 0.5193 0.6036242
incarc_rate 3.1370e-05 3.4425e-04 0.0911 0.9274087
density -4.4745e-02 1.9626e-01 -0.2280 0.8196970
avginc 1.4357e-02 2.4523e-02 0.5854 0.5583737
pop 1.6388e-05 2.5682e-02 0.0006 0.9994910
pw1064 -1.2832e-02 3.2439e-02 -0.3956 0.6924960
pb1064 1.4108e-02 8.3232e-02 0.1695 0.8654348
pm1029 1.0460e-01 7.2277e-02 1.4473 0.1481059
factor(year)78 3.2850e-02 2.1476e-02 1.5296 0.1264015
factor(year)79 1.3759e-01 3.1800e-02 4.3268 1.652e-05 ***
factor(year)80 2.4341e-01 4.5016e-02 5.4072 7.861e-08 ***
factor(year)81 2.7371e-01 5.0377e-02 5.4332 6.824e-08 ***
factor(year)82 2.1599e-01 6.3776e-02 3.3867 0.0007325 ***
factor(year)83 1.2082e-01 8.5851e-02 1.4073 0.1596326
factor(year)84 7.8831e-02 1.0538e-01 0.7481 0.4545869
factor(year)85 1.1315e-01 1.2601e-01 0.8980 0.3694062
factor(year)86 1.8957e-01 1.5064e-01 1.2584 0.2085230
factor(year)87 1.5722e-01 1.6722e-01 0.9402 0.3473427
factor(year)88 1.9276e-01 1.8603e-01 1.0362 0.3003544
factor(year)89 2.4873e-01 2.1195e-01 1.1736 0.2408270
factor(year)90 3.5098e-01 2.6423e-01 1.3283 0.1843509
factor(year)91 4.6685e-01 2.7642e-01 1.6889 0.0915220 .
factor(year)92 4.6332e-01 2.9222e-01 1.5856 0.1131306
factor(year)93 4.7970e-01 3.0519e-01 1.5718 0.1162908
factor(year)94 4.9438e-01 3.2022e-01 1.5439 0.1229143
factor(year)95 4.9402e-01 3.3055e-01 1.4945 0.1353297
factor(year)96 4.3416e-01 3.4698e-01 1.2513 0.2111054
factor(year)97 3.6524e-01 3.5464e-01 1.0299 0.3032928
factor(year)98 2.6771e-01 3.6540e-01 0.7327 0.4639202
factor(year)99 1.8947e-01 3.8075e-01 0.4976 0.6188516
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg10 <- plm(log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029 + factor(year), data = df, model = "pooling", effect = "individual")pFtest(reg9, reg10)
F test for individual effects
data: log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 174.38, df1 = 50, df2 = 1092, p-value < 2.2e-16
alternative hypothesis: significant effectslinearHypothesis(reg9, matchCoefs(reg9, "year"), vcov = vcov9)Linear hypothesis test
Hypothesis:
factor(year)78 = 0
factor(year)79 = 0
factor(year)80 = 0
factor(year)81 = 0
factor(year)82 = 0
factor(year)83 = 0
factor(year)84 = 0
factor(year)85 = 0
factor(year)86 = 0
factor(year)87 = 0
factor(year)88 = 0
factor(year)89 = 0
factor(year)90 = 0
factor(year)91 = 0
factor(year)92 = 0
factor(year)93 = 0
factor(year)94 = 0
factor(year)95 = 0
factor(year)96 = 0
factor(year)97 = 0
factor(year)98 = 0
factor(year)99 = 0
Model 1: restricted model
Model 2: log(rob) ~ shall + incarc_rate + density + avginc + pop + pw1064 +
pb1064 + pm1029 + factor(year)
Note: Coefficient covariance matrix supplied.
Res.Df Df Chisq Pr(>Chisq)
1 1114
2 1092 22 580.35 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Table 3
reg11 <- plm(log(mur) ~ shall, data = df, model = "pooling")vcov11 <- vcovHC(reg11, method = "arellano", type = "HC1", cluster = "group")coeftest(reg11, vcov = vcov11)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.897556 0.092263 20.5669 < 2.2e-16 ***
shall -0.473372 0.147464 -3.2101 0.001363 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg12 <- plm(log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029, data = df, model = "pooling")vcov12 <- vcovHC(reg12, method = "arellano", type = "HC1", cluster = "group")coeftest(reg12, vcov = vcov12)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.485593 1.973297 -1.2596 0.2080610
shall -0.313174 0.098108 -3.1921 0.0014500 **
incarc_rate 0.002097 0.000456 4.5986 4.718e-06 ***
density 0.039667 0.039517 1.0038 0.3156848
avginc -0.077258 0.026789 -2.8839 0.0039997 **
pop 0.041617 0.011814 3.5229 0.0004435 ***
pw1064 0.047080 0.028322 1.6623 0.0967198 .
pb1064 0.130764 0.060614 2.1573 0.0311859 *
pm1029 0.065531 0.035823 1.8293 0.0676121 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg13 <- plm(log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029, data = df, model = "within", effect = "individual")vcov13 <- vcovHC(reg13, method = "arellano", type = "HC1", cluster = "group")coeftest(reg13, vcov = vcov13)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall -0.06080998 0.03659889 -1.6615 0.09689 .
incarc_rate -0.00036000 0.00041898 -0.8592 0.39040
density -0.67071318 0.39187400 -1.7116 0.08726 .
avginc 0.02431136 0.01552339 1.5661 0.11761
pop -0.02570538 0.02014519 -1.2760 0.20222
pw1064 0.01033127 0.01275067 0.8103 0.41797
pb1064 0.03070086 0.07735459 0.3969 0.69153
pm1029 0.03923836 0.02138361 1.8350 0.06678 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1pFtest(reg13, reg12)
F test for individual effects
data: log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 62.003, df1 = 72, df2 = 1092, p-value < 2.2e-16
alternative hypothesis: significant effectsvcov14 <- vcovHC(reg14, method = "arellano", type = "HC1", cluster = "group")coeftest(reg14, vcov = vcov14)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
shall -0.01495237 0.03786313 -0.3949 0.6929895
incarc_rate -0.00011640 0.00035950 -0.3238 0.7461542
density -0.54426353 0.31607167 -1.7220 0.0853595 .
avginc 0.05664917 0.01639211 3.4559 0.0005696 ***
pop -0.03207693 0.02077490 -1.5440 0.1228724
pw1064 -0.00048933 0.01990607 -0.0246 0.9803929
pb1064 0.02198327 0.07506729 0.2928 0.7696943
pm1029 0.06919411 0.04138223 1.6721 0.0947964 .
factor(year)78 -0.00071950 0.03195386 -0.0225 0.9820399
factor(year)79 0.05924812 0.03080724 1.9232 0.0547162 .
factor(year)80 0.09018140 0.04065306 2.2183 0.0267383 *
factor(year)81 0.10215430 0.05055996 2.0205 0.0435796 *
factor(year)82 0.02240978 0.05761214 0.3890 0.6973693
factor(year)83 -0.03143848 0.06343026 -0.4956 0.6202492
factor(year)84 -0.13591915 0.07095514 -1.9156 0.0556815 .
factor(year)85 -0.08661433 0.08485126 -1.0208 0.3075856
factor(year)86 -0.01227518 0.09181401 -0.1337 0.8936675
factor(year)87 -0.02903376 0.09895500 -0.2934 0.7692694
factor(year)88 -0.01745941 0.11850879 -0.1473 0.8829020
factor(year)89 -0.01456168 0.13080044 -0.1113 0.9113771
factor(year)90 0.05999806 0.16334464 0.3673 0.7134593
factor(year)91 0.10530717 0.17375992 0.6060 0.5446076
factor(year)92 0.06810019 0.18103181 0.3762 0.7068576
factor(year)93 0.15442971 0.18793912 0.8217 0.4114266
factor(year)94 0.04426479 0.19524584 0.2267 0.8206893
factor(year)95 0.05566014 0.19694628 0.2826 0.7775249
factor(year)96 -0.01570898 0.21044014 -0.0746 0.9405083
factor(year)97 -0.12218236 0.21651375 -0.5643 0.5726544
factor(year)98 -0.18633807 0.23099549 -0.8067 0.4200299
factor(year)99 -0.25542854 0.23965595 -1.0658 0.2867436
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1reg15 <- plm(log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 + pb1064 + pm1029 + factor(year), data = df, model = "pooling")pFtest(reg14, reg15)
F test for individual effects
data: log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 + ...
F = 76.862, df1 = 50, df2 = 1092, p-value < 2.2e-16
alternative hypothesis: significant effectslinearHypothesis(reg14, matchCoefs(reg14, "year"), vcov = vcov14)Linear hypothesis test
Hypothesis:
factor(year)78 = 0
factor(year)79 = 0
factor(year)80 = 0
factor(year)81 = 0
factor(year)82 = 0
factor(year)83 = 0
factor(year)84 = 0
factor(year)85 = 0
factor(year)86 = 0
factor(year)87 = 0
factor(year)88 = 0
factor(year)89 = 0
factor(year)90 = 0
factor(year)91 = 0
factor(year)92 = 0
factor(year)93 = 0
factor(year)94 = 0
factor(year)95 = 0
factor(year)96 = 0
factor(year)97 = 0
factor(year)98 = 0
factor(year)99 = 0
Model 1: restricted model
Model 2: log(mur) ~ shall + incarc_rate + density + avginc + pop + pw1064 +
pb1064 + pm1029 + factor(year)
Note: Coefficient covariance matrix supplied.
Res.Df Df Chisq Pr(>Chisq)
1 1114
2 1092 22 440.06 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1