Report
Christina De Cesaris
3/13/2021
Abstract
Using panel data regarding violent crime, demographics, and gun control laws in the United States between 1977-1999, the effect of shall-carry gun laws on rates of violent crime among 50 states and the District of Colombia were investigated using a Propensity Weighted Linear Mixed Effects model. The model attempted to answer if a causal relationship could be determined between shall-carry laws and violent crime using robbery rates, murder rates, population density, and the percentage of young males between 10-29 years old as confounders. State and year were defined as random effects by the model. To reduce the selection bias introduced by confounders, propensity scores were calculated using a logistic regression model and weights were derived using the Inverse Probability of Treatment Weights method (IPTW). The final weighted model indicated a negative correlation between law=‘yes’ (p < 0.1) and violent crime rates. The significance of law was reduced after adding the weights suggesting some selection bias was removed from weighting. Murder (p < 0.05), robbery (p < 0.001), and young male population percent (p < 0.001) were found positively correlated with violent crime. Because selection bias was reduced but not eliminated, a causal relationship could not be determined.
Introduciton and Background
The controversy surrounding guns has become a hotbed for political discourse in recent years. Some have linked gun ownership to increases in homicides, suicides, and organized criminal activity (Miller et. al, 2002). A 2004 publication found that owning a gun makes a person three to five times more likely to die by violence (Zimring). Others have suggested that gun ownership does not lead to a net increase in crime but rather increased crime leads to more gun ownership. In the case of this argument, authors argued that gun control laws would be ineffective in reducing violent crime (Kleck et. al, 2002).
In this project, we investigate the relationship between population demographics, various crime rates, and shall carry gun laws across different regions of the United States between 1977-1999. Shall carry gun laws are laws which allow only licensed, trained, and background checked citizens to carry concealed guns in specific states.
We seek to answer the following questions:
What is the effect of shall carry laws on violent crime in the United States?
Is this relationship casual?
It is important to note that a key point regarding the theory of causal interference is the factor variable in question should be changeable. In the case of this project, the shall-carry law in effect is a factor which can be implemented or rescinded by law makers making it ideal for this analysis.
Dataset
The dataset was downloaded using the AER package in R. It consists of 13 variables and of 1173 observations. The data spanned twenty-two years and contained for fifty American states and the District of Colombia. The variables are as described from the source as follows: (Stock & Watson, 2007)
| variable | desciption |
|---|---|
| state | factor indicating state. |
| year | factor indicating year. |
| violent | violent crime rate (incidents per 100,000 members of the population). |
| murder | murder rate (incidents per 100,000). |
| robbery | robbery rate (incidents per 100,000). |
| prisoners | incarceration rate in the state in the previous year (sentenced prisoners per 100,000 residents; value for the previous year). |
| afam | percent of state population that is African-American, ages 10 to 64. |
| cauc | percent of state population that is Caucasian, ages 10 to 64. |
| male | percent of state population that is male, ages 10 to 29. |
| population | state population, in millions of people. |
| income | real per capita personal income in the state (US dollars). |
| density | population per square mile of land area, divided by 1,000. |
| law | factor. Does the state have a shall carry law in effect in that year? |
Exploratory Analysis
The full, raw dataset is provided as an interactive table for viewer purposes:
Summary averages of all numerical variables were grouped by state and year and constructed into respective interactive data tables. Washington DC and the year 1993 have the highest average rates for violent crime.
The summary statistics: mean, median, standard deviation, minimum, and maximum of our response: violent crime rates were calculated for year and state respectively.
A scatter plot matrix was created for numerical variables along with their calculated pearson correlation coefficients to assess the relationships between potential predictors and the response.
From the matrix, violent crime has high positive correlation with murder rates, robbery rates, and prison rates. Violent crime has a weak negative correlation with young male populations.
A plot comparing violent crime rates to whether or not a shall-carry law is in effect shows lower average violent crime rates compared to when the law is not in effect. But whether this is significant or not is unclear for the graph alone.
Interactive time-series plots of year violent crime rates for each state reveal Washington DC as an outlier. We considered dropping Washington DC, but this data was not a result of misrecording or other oddity. Washington DC did have disproportional amounts of crime in the 1980s and 1990s so the level was kept.
In the case of all interactive plots in this project, singularly clicking on a legend item will remove it from the graph. Double clicking on a legend item will plot it individually.
A boxplot of violent crime rates for each year again shows significant outliers which all happen to be the values from Washington DC. Violent crime rates appear to fluctuate over time with the sharpest increase occuring in the early 1990s.
A plot of average crime rates (violent, robbery, murder) and incarceration rates shows the positive correlation between crime rates as they fluctuate overtime. Interestingly, the average incarceration rates continue to increase even when crime rates decrease leading us to believe there are other reasons besides crime causing this vast increase in prisioners.
Ethnic demographic percentages other than Caucasian and African American were not included in this dataset. We manually calculated the percentages of other ethnicities using those provided. As a result, the ‘other’ ethnicities category mirrors that of the Caucasian category. Because of these inconsistencies, ethnicities percentages were not considered for our model.
However, the young male demographic percentage was chosen for our model. Isolating the young male demographic into its own graph (done by double clicking ‘Males…’ in the legend) shows a decrease in the young male demographic between 1977 and 1999.
Inferential Analysis
Linear Mixed Model Statistical Form
We propose a propensity weighted linear mixed model to determine if a casual relationship exists between a state’s shall-carry law status and violent crime rates.
It is under our assumption that multiple factors contribute to the implementation of a shall-carry law, and we believe density, robbery rate, murder rate, young male demographics, and incarceration rates have confounding effects on both law implementation and violent crime rates.
The lmerTest package was used to fit our Linear Mixed-Effects Models. Since the data contained a multitude of observations for different states over the course of twenty-two years, we assumed existing random, immeasurable effects were present in our data. In this model, we define state and year effects as random effects we cannot control.
Our Linear Mixed-Effects Model takes the following form:
\[ \vec y = { \boldsymbol X \vec \beta} + {\boldsymbol Z \vec u} + {\vec \epsilon}\] Where
\(\vec y\) represents the \(1173 \times 1\) vector of violent crime rates.
\(X\) is a \(1173 \times 7\) matrix of our fixed effects variables:law, male, density, murder rates, robbery rates and incarceration rates. The first entry of \(X\) is \(1\) to include the intercept.
$ $ is a \(7 \times 1\) vector of our fixed effects coefficients and intercept.
\(Z\) is a \(1173 \times 2\) matrix for the observed values for the \(2\) random variables:state and year.
\(\vec u\) is a \(2 \times 1\) represents the unobserved random effects for our random variables (analogous to \(\vec \beta\) for fixed variables)
Finally, \(\vec \epsilon\) is a \(1173 \times 1\) vector which captures the random portion of \(y\) not included by the rest of the model.
Note that the weights do not appear in the mathematical form of this model. This is because the type of weight applied is conditional on the observed value of \(X_{law}\). More is explained in the following weights section.
Linear Mixed Model Assumptions:
Homogeneity of variance
Errors are normally distributed \(\epsilon \sim N(0, \sigma_e^2)\)
Random variable effects are normally distributed \(u \sim N(0, \sigma_u^2)\)
No extreme multicoliniarity
(McCullagh & Nelder, 1989)
Weights
Propensity score methods attempt to account for the probability of a certain treatment occurring given other observed covariates in the model. In an ideal situation, the use of propensity score methods will reduce selection bias introduced by confounding effects in our model.
We assume that young male demographics, density, murder rates, robbery rates and incarceration rates influence both the value of law and violence rates. For example, it is possible that higher murder rates will influence both violent crime and whether law makers decided to enact a shall-carry law.
And so, the weights derived from our propensity scores seek to reduce the selection bias young male demographics, density, murder rates, robbery rates and incarceration rates have on law, so we have a clearer idea of the average effect of law on violent crime rates.
(Olmos & Govindasmy, 2015)
To calculate weights, first the probability of law = ‘yes’or law=’no’ was estimated given our covaraites. In this project, probabilities were estimated using a logistic regression model defined as follows:
\[P(y_{law}^{(i)}=yes)=\frac{1}{1+exp(-(\beta_{0}+\beta_{1}x^{(i)}_{1}+\ldots+\beta_{p}x^{(i)}_{p}))}\] where \(1 \ldots p\) represent \(male, robbery, murder, \& \ density\) respectively.
The estimated probabilities from logistic regression are then converted from probabilities into weights using the Inverse Probability of Treatment Weights method which follows as:
Scores for where law = ‘yes’: \(w_i = \frac{1}{\hat{p_i}}\)
Scores for where law = ‘no’: \(w_i = \frac{1}{(1-\hat{p_i})}\)
These weights were then applied to the linear mixed model.
Applying the Model
The distribution of the response was plotted and found to be right skewed. As a result, a log transformation was applied to the response to make the distribution more normal.
An unweighted linear mixed model was initially fit using law, male, robbery, murder, and density as predictors. As mentioned previously, state and year were considered random.
From the unweighted model:
the intercept, robbery, and male were all significantly (p<0.001)and positively related with log(violent)
law = ‘yes’ was significantly (level: 0.05) and negatively correlated with log(violent)
density and murder were not significant
#before weights
mod_pre = lmer(log(violent)~law+male+robbery+murder+density+(1|state)+(1|year), data=Guns)
summary(mod_pre)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(violent) ~ law + male + robbery + murder + density + (1 |
## state) + (1 | year)
## Data: Guns
##
## REML criterion at convergence: -1078.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9155 -0.5158 0.0310 0.5247 5.1815
##
## Random effects:
## Groups Name Variance Std.Dev.
## state (Intercept) 0.2386 0.4885
## year (Intercept) 0.0290 0.1703
## Residual 0.0160 0.1265
## Number of obs: 1173, groups: state, 51; year, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.065e+00 1.632e-01 3.419e+02 31.031 < 2e-16 ***
## lawyes -3.752e-02 1.497e-02 1.124e+03 -2.505 0.0124 *
## male 4.660e-02 9.057e-03 4.652e+02 5.145 3.96e-07 ***
## robbery 1.350e-03 8.816e-05 1.121e+03 15.311 < 2e-16 ***
## murder 3.213e-03 1.347e-03 9.881e+02 2.384 0.0173 *
## density -6.052e-02 3.446e-02 2.039e+02 -1.756 0.0806 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lawyes male robbry murder
## lawyes 0.150
## male -0.874 -0.191
## robbery 0.014 -0.105 -0.062
## murder -0.039 0.089 -0.031 -0.337
## density 0.010 0.024 -0.095 -0.305 0.421Before constructing weights, balancing analysis was performed to evaluate the effect of our covariates with regard to our treatment group law.
From the pre-weighted balance analysis, we see there is significant imbalance between law and all other covariates. In an ideal situation, the addition of propensity weights to these models will decrease the imbalance, and therefore make the relationship between our treatment and covariates insignficant.
Pre-Weighted Balance Analysis [[1]]
|
[[2]]
|
[[3]]
|
[[4]]
|
As previously stated, the probabilities use in weighting were calculated by a logistic regression model using our treatment (law) as the response and our other fixed effects as predictors.
##
## Call:
## glm(formula = law ~ male + robbery + murder + density, family = binomial,
## data = Guns)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6513 -0.7063 -0.3844 0.0000 2.3971
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.942724 0.928093 10.713 < 2e-16 ***
## male -0.606987 0.058052 -10.456 < 2e-16 ***
## robbery -0.004295 0.001580 -2.718 0.00657 **
## murder -0.069072 0.035449 -1.949 0.05135 .
## density -4.164473 0.721933 -5.769 8e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1300.8 on 1172 degrees of freedom
## Residual deviance: 1013.8 on 1168 degrees of freedom
## AIC: 1023.8
##
## Number of Fisher Scoring iterations: 8After the probability scores were converted into pscores, the pscores were trimmed to reduce instability. The distribution of all pscores was plotted, but this in itself is not particularly useful besides showing the distribution of scores.
On the other hand, the plot of the overlapping distributions of law=‘yes’ and law=‘no’ reveal there is similarities between the distributions. The observations corresponding to the pscores in this purple area have similar covariate observations for different values of law. In some ways, the weights allow us to look at the model with priority on these covaraite ‘balanced’ observations. However, given the amount of non-overlap, we can infer that not all covaraites will be balanced.
prob= mod_log$fitted #probabilities
pscore = ifelse(Guns$law=='yes', prob, 1-prob)
Guns$pscore=pscore
weight=1/pscore
#distribution left skewed, not terrible extremes, consider trimming [0.01,0.99]
Guns$pscore=pscore
Guns$pscore=ifelse(pscore>0.99, 0.99, Guns$pscore)
Guns$pscore=ifelse(pscore<0.01, 0.01, Guns$pscore)
summary(Guns$pscore)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05653 0.54452 0.78654 0.71547 0.93182 0.99000
Unsurprisingly, the addition of weights did not provide perfect balance. Balance was achieved with regard to the covarite male, and imbalance was reduced with regard to the covariate murder. The imbalance resulting from robbery and density was decreased (as seen by the increased p-values) but this change was not significant.
Weighted Balance Assessment [[1]]
|
[[2]]
|
[[3]]
|
[[4]]
|
The final, weighted linear mixed model was fit using law, male, robbery, murder, and density as predictors. Again, state and year were considered random varaibles.
From the weighted model:
the intercept, robbery, and male were all significantly (p<0.001)and positively related with log(violent)
law = ‘yes’ was significantly (level: 0.1) and negatively correlated with log(violent). The effect of law was found to be less signficant with the addition of weights. This indicates some bias was effectively removed from the model
density and murder were not significant
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(violent) ~ law + male + robbery + murder + density + (1 |
## state) + (1 | year)
## Data: Guns
## Weights: weight
##
## REML criterion at convergence: -990.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8448 -0.4846 0.0077 0.4710 4.7924
##
## Random effects:
## Groups Name Variance Std.Dev.
## state (Intercept) 0.21428 0.4629
## year (Intercept) 0.03180 0.1783
## Residual 0.02663 0.1632
## Number of obs: 1173, groups: state, 51; year, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.919e+00 1.605e-01 3.999e+02 30.654 < 2e-16 ***
## lawyes -2.380e-02 1.405e-02 1.125e+03 -1.694 0.0905 .
## male 5.320e-02 8.887e-03 5.952e+02 5.987 3.71e-09 ***
## robbery 1.579e-03 1.024e-04 1.131e+03 15.418 < 2e-16 ***
## murder 3.834e-03 1.592e-03 1.019e+03 2.409 0.0162 *
## density -5.614e-02 3.737e-02 1.341e+02 -1.502 0.1354
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lawyes male robbry murder
## lawyes 0.111
## male -0.875 -0.146
## robbery -0.004 -0.152 -0.050
## murder -0.043 0.088 -0.028 -0.360
## density -0.025 0.038 -0.057 -0.308 0.336Sensitivity Analysis
Sensitivity analysis allows us to evaluate the sufficiency of the model according to how well it follows initial assumptions.
For this model, the assumption of homogeneity of variance is followed given the residuals plot. As well, the linear assumption appears to hold as there is no distinct or clear pattern indicated by the residuals plot.
## Loading required namespace: testthat
Another important linear assumption is the normality assumption. We assumed the randomness of our model follows a normal distribution. From the QQ plot of residuals, we see this assumption is roughly met. The plot indicates there are heavy tails possibly as a result of outliers. The distribution of residuals as well appears roughly normal. 
While the random effects from our model are in some way incorporated into the previous QQ plot, it is a good idea to check if our random variables year and state hold up to the normality assumption as well. From the quantile plots below, it appears both year and state do not deviate greatly from the normal distribution.
We also assumed the independence of fixed variables in our model. After calculating the variance inflation factors (VIF) for our model, we see there aren’t any instances of high multicoliniarity (usually indicated by VIF > 2.5 in conservative cases)
vif(mod_post)## law male robbery murder density
## 1.050742 1.034305 1.235785 1.230632 1.186698Conclusions
Controversy regarding gun control laws and their effectiveness has been the source of several political conflicts. The purpose of this project is to analyze the issue from a statistical perspective in order to see if there exists a casual relationship between shall-carry gun laws and violent crime rates in the USA.
This project used the panel Guns dataset from the AER package and analyzed the relationship between young male population percentage, population density, murder rates, shall-carry laws, and robbery rates and violent crime. The data was collected at the state level over the period between 1977-1999.
It was assumed that young male population percentage, population density, murder rates, and robbery rates had confounding effects on both violent crime and the implementation of shall-carry laws.
The main question of interest was whether the relationship between shall-carry laws and violent crime was casual. A propensity score weighted linear mixed effects model was fit to answer this question. State and year were regarded as random effects. Violent crime rate was the response and young male population percentage, population density, murder rates, shall-carry laws, and robbery rates were predictors.
To reduce selection bias introduced by our confounders on shall-carry laws, a logistic regression model using law as response and confounders as predictors was fit. These fitted values (propensity scores) were then converted into weights and applied to the final linear mixed model.
Balancing analysis suggests the weights only sucessfully achieved balanced between law and young male percentage. While imbalance may have been reduced among the treatment (law) and other confounders, it was far from eliminated.
The final weighted model indicated a negative correlation between law=‘yes’ and violent crime rates. The significance of law was reduced after adding the weights suggesting some selection bias was removed from weighting. Murder, robbery, and young male population percent were positively correlated with violent crime.
From the sensitivity analysis, model assumptions were roughly followed and the model was deemed adequate. It is possible that the assumption of no unobserved confounders was violated as not all data variables were included in the model.
Overall, we can NOT infer a causal relationship between shall-carry laws and violent crime rates because selection bias remains in the model.
The shall-carry law was also not highly, significantly correlated with violent crime rates which may grant merit to the cited studies which claimed gun control laws were ineffective in reducing violence (Kleck et. al, 2002).
References
Brookhart, M. A., Wyss, R., Layton, J. B., & Stürmer, T. (2013). Propensity score methods for confounding control in nonexperimental research. Circulation. Cardiovascular quality and outcomes, 6(5), 604–611. https://doi.org/10.1161/CIRCOUTCOMES.113.000359
Kleck, G., Kovandzic, T., & Bellows, J. (2016). Does Gun Control Reduce Violent Crime? Criminal Justice Review, 41(4), 488–513. https://doi.org/10.1177/0734016816670457
Kleiber, C. & Zeileis, A. (2008). Applied Econometrics with R. Springer-Verlag.
McCullagh, P, & Nelder, J. A. (1989). Generalized Linear Models, 2nd ed. Chapman & Hall/CRC Press.
Miller, M., Azrael, D., & Hemenway, D. (2002). Rates of household firearm ownership and homicide across US regions and states, 1988-1997. American journal of public health, 92(12), 1988–1993. https://doi.org/10.2105/ajph.92.12.1988
Olmos, A., & Govindasamy, P. (2015). A Practical Guide for Using Propensity Score Weighting in R. Practical Assessment, Research & Evaluation, 20(13).
Stock, J.H. and Watson, M.W. (2007). Introduction to Econometrics, 2nd ed. Boston: Addison Wesley.
Zimring, F. (2004). Firearms, Violence, and the Potential Impact of Firearms Control. Journal of Law, Medicine & Ethics, 32(1), 34-37. doi:10.1111/j.1748-720X.2004.tb00446.x
Code Appendix
# plotting
library(redres)
library(ggplot2)
library(plotly)
# analysis
library(MASS)
library(tidyverse)
library(lmerTest)
library(scales)
# formatting
library(grid)
library(gridExtra)
library(stargazer)
library(knitr)
# data
library(AER)
library(DT)
data(Guns)
variable = c("state",
"year","violent","murder","robbery","prisoners","afam","cauc","male","population","income","density","law")
desciption=c("factor indicating state.","factor indicating year.",'violent crime rate (incidents per 100,000 members of the population).',"murder rate (incidents per 100,000).","robbery rate (incidents per 100,000).",'incarceration rate in the state in the previous year (sentenced prisoners per 100,000 residents; value for the previous year).',"percent of state population that is African-American, ages 10 to 64.","percent of state population that is Caucasian, ages 10 to 64.","percent of state population that is male, ages 10 to 29.","state population, in millions of people.","real per capita personal income in the state (US dollars).","population per square mile of land area, divided by 1,000.","factor. Does the state have a shall carry law in effect in that year?")
vars=cbind(variable, desciption)
kable(vars)
av_st = Guns %>%
group_by(state) %>% # average all state
summarise_all(mean) %>%
mutate_if(is.numeric,funs(round(.,2)))
av_st
#av_st
av_st_r = Guns %>%
group_by(state) %>% # average all state
summarise(mean(violent),median(violent), sd(violent), max(violent), min(violent))%>%
mutate_if(is.numeric,funs(round(.,2)))
av_st_r
av_yr_r = Guns %>%
group_by(year) %>% # average all state
summarise(mean(violent),median(violent), sd(violent), max(violent), min(violent))%>%
mutate_if(is.numeric,funs(round(.,2)))
av_yr_r
av_styr = Guns %>% # average each year each state
group_by(year,state) %>%
summarise_all(mean)%>%
mutate_if(is.numeric,funs(round(.,2)))
av_styr
av_yr = Guns %>% # average all states each years
group_by(year)%>%
summarise_all(mean)%>%
mutate_if(is.numeric,funs(round(.,2)))
#delete null cols
av_yr1=av_yr
av_yr=dplyr::select(av_yr,-c(12,13))
av_yr1=dplyr::select(av_yr1,-c(7,8,12,13))
av_st=dplyr::select(av_st,-c(2,7,8,13))
all_dat=datatable(Guns,
class = 'cell-border stripe',
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;','Table 1: ',
htmltools::em(
'Full, Raw Dataset')))
names_yr=c("Year","Violent Crime Rate","Murder Rate","Robbery Rate","Prisoners","Percent Male","Population", "Income","Density")
dt_av_yr=datatable(av_yr1,
class = 'cell-border stripe',
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;','Table 2: ',
htmltools::em(
'Numeric Variable Averages By Year')),
colnames = c(
names_yr))
names_st=c("State","Violent Crime Rate","Murder Rate","Robbery Rate","Prisoners","Percent Male","Population", "Income","Density")
dt_av_st=datatable(av_st,
class = 'cell-border stripe',
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;','Table 3: ',
htmltools::em(
'Numeric Variable Averages By State')),
colnames = c(
names_st))
#RESPONSE
names_av_st_r=c("State","Mean Violence Rate", "Median Violence Rate", "Standard Deviation Violence Rate","Max Violence Rate","Min Violence Rate" )
dt_av_st_r=datatable(av_st_r,
class = 'cell-border stripe',
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;','Table 4: ',
htmltools::em(
'Violent Crime Summary Statistics by State')),
colnames = c(
names_av_st_r))
names_av_yr_r=c("Year","Mean Violence Rate", "Median Violence Rate", "Standard Deviation Violence Rate","Max Violence Rate","Min Violence Rate" )
dt_av_yr_r=datatable(av_yr_r,
class = 'cell-border stripe',
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;','Table 5: ',
htmltools::em(
'Violent Crime Summary Statistics by Year')),
colnames = c(
names_av_yr_r))
all_dat
dt_av_yr
dt_av_st
dt_av_st_r
dt_av_yr_r
#sapply(Guns, class) ensure all variables are as should be
panel_cor <- function(x, y){ #used to generate lower panel of pearson correlation coefficents for pairs plot
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- round(cor(x, y, use="complete.obs"), digits=2)
txt <- paste0("R = ", r)
cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor )
}
pairs(Guns[sapply(Guns, is.numeric)],
lower.panel = panel_cor)
#level plot
boxplot(Guns$vio~Guns$law,
xlab = "Is Shall-Carry Law in Effect?",
ylab = "Violent Crime Rate",
main="Level Comparision: Violent Crime Rate vs Law",
col=c('red','lightblue'))
v=Guns %>%
ggplot(
aes(x=year,
y=violent,
group=state,
color=state)) +
geom_line() +
labs(
title="Yearly Trends in Violent Crime by State",
x="School ID",
y = "Violent Crime Rate per 100K Poeple") +
theme(
text = element_text(
size=8,
face = 'bold'),
axis.text.x = element_text(angle=75,hjust=1),
axis.ticks.x = element_line(size=0.5)) +
scale_x_discrete(
guide = guide_axis(n.dodge =1),
expand=c(0, 0)) +
scale_fill_discrete(name = "State")
ggplotly(v)
m = Guns %>%
ggplot(
aes(x=year,
y=(violent),
group=year,
col=year)) +
geom_boxplot()+
labs(
title="USA Violent Crime by Year",
x="Year",
y = "Violent Crime Rate per 100K People") +
theme(
text = element_text(
size=8,
face = 'bold'),
axis.text.x = element_text(angle=75,hjust=1),
axis.ticks.x = element_line(size=0.5)) +
scale_x_discrete(
guide = guide_axis(n.dodge =1),
expand=c(0, 0)) +
scale_fill_discrete(name = "year")
ggplotly(m) # DC is a gross outlier, consider removal
ggplotly(m)
crime_fig <- plot_ly(av_yr, x= ~year) %>%
add_lines(
y= ~(violent),
name = "Violent Crime",
type = 'scatter',
mode = 'lines+markers') %>%
add_lines(
y= ~(robbery),
name = "Robberies",
visible = T) %>%
add_lines(
y= ~(murder),
name = "Murder",
visible = T) %>%
add_lines(
y= ~(prisoners),
name = "Incarceration",
visible = T) %>%
layout(
title = "Average Crime and Incarceration Rates in the USA",
xaxis = list(title = "Year",
tickangle = 75),
yaxis = list(title = "Average per 100k People",
margin = margin))
crime_fig
av_yr$other= (100-(av_yr$afam+av_yr$cauc))
demo_fig <- plot_ly(av_yr, x= ~year) %>%
add_lines(
y= ~(male),
name = "Males: Ages 10-29",
type = 'scatter',
mode = 'lines+markers') %>%
add_lines(
y= ~(afam),
name = "African Americans: Ages 10-64",
visible = T) %>%
add_lines(
y= ~(cauc),
name = "Caucation Americans: Ages 10-64",
visible = T) %>%
add_lines(
y= ~(other),
name = "Other American Ethnicities: Ages 10-64",
visible = T) %>%
layout(
title = "Average Yearly Population Demographics",
xaxis = list(title = "Year",
tickangle = 75),
yaxis = list(title = "Average percent of Population",
margin = margin),
legend)
demo_fig
#distribution of response
par(mfrow=c(1,2))
hist(Guns$violent,
xlab="Violent Crime Rate",
main= "Raw",
col="yellow") #TODO, ADD AXIS
hist(log(Guns$violent),
xlab="log(Violent Crime Rate)",
main="Log-Transformed",
col="lightgreen")
#before weights
mod_pre = lmer(log(violent)~law+male+robbery+murder+density+(1|state)+(1|year), data=Guns)
summary(mod_pre)
bal1= lm(male~law, Guns)
kpre1=kable(summary(bal1)[4],
caption = "lm(male ~ law)")
bal2= lm(robbery~law, Guns)
kpre2=kable(summary(bal2)[4],
caption = "lm(robbery ~ law)")
bal3= lm(murder~law, Guns)
kpre3=kable(summary(bal3)[4],
caption = "lm(murder ~ law)")
bal4= lm(density~law, Guns)
kpre4=kable(summary(bal4)[4],
caption = "lm(density ~ law)")
bal_pre = list(kpre1,kpre2, kpre3, kpre4)
bal_pre
#obtain weights using logistic regression and inverse weighted propensity scores
mod_log = glm(law~male+robbery+murder+density,family=binomial, Guns)
summary(mod_log)
prob= mod_log$fitted #probabilities
pscore = ifelse(Guns$law=='yes', prob, 1-prob)
Guns$pscore=pscore
weight=1/pscore
#distribution left skewed, not terrible extremes, consider trimming [0.01,0.99]
Guns$pscore=pscore
Guns$pscore=ifelse(pscore>0.99, 0.99, Guns$pscore)
Guns$pscore=ifelse(pscore<0.01, 0.01, Guns$pscore)
summary(Guns$pscore)
col.alpha<-function(color, alpha){ #changes opacity of color by factor of alpha
code=col2rgb(color)/256
rgb(code[1],code[2],code[3],alpha)
}
hist(pscore,
breaks=50,
col="violet",
freq=F,
ylim=c(0,5),
xlab="Propensity Score",
ylab="",
main="Propensity Score Distribution: Ignoring Level")
#better hist
hist(Guns$pscore[Guns$law=='yes'],
breaks=30,
col=col.alpha("red",.6),
freq=F,
ylim=c(0,5),
xlab="Propensity Score",
ylab="",
main="Propensity Score Distribution: Treatment vs Control")
hist((Guns$pscore[Guns$law=='no']),
breaks=30,
col=col.alpha("lightblue",.6),
freq=F, ylim=c(0,5),
xlab="Propensity Score",
ylab="",
main="",
add=T)
legend(.4,5,
legend=c("Law: Yes", "Law: No"),
col=c("red", "blue"),
fill = c("red", "lightblue") )
bal1= lm(male~law, Guns, weights=weight)
kpost1=kable(summary(bal1)[5],
caption = "lm(male ~ law)")
bal2= lm(robbery~law, Guns, weights=weight)
kpost2=kable(summary(bal2)[5],
caption = "lm(robbery ~ law)")
bal3= lm(murder~law, Guns, weights=weight)
kpost3=kable(summary(bal3)[5],
caption = "lm(murder ~ law)")
bal4= lm(density~law, Guns, weights=weight)
kpost4=kable(summary(bal4)[5],
caption = "lm(density ~ law)")
bal_post = list(kpost1,kpost2, kpost3, kpost4)
bal_post
mod_post = lmer(log(violent)~law+male+robbery+murder+density+(1|state)+(1|year), data=Guns, weights=weight)
summary(mod_post)
sm=summary(mod_post)
plot_redres(mod_post, type = "std_cond")
plot_resqq(mod_post)
hist(sm$residuals, breaks=30, main="Residuals Distribution", xlab="Residuals")
plot_ranef(mod_post)
vif(mod_post)
sessionInfo()sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19041)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DT_0.17 AER_1.2-9 survival_3.2-7 sandwich_3.0-0
## [5] lmtest_0.9-38 zoo_1.8-8 car_3.0-10 carData_3.0-4
## [9] knitr_1.31 stargazer_5.2.2 gridExtra_2.3 scales_1.1.1
## [13] lmerTest_3.1-3 lme4_1.1-26 Matrix_1.3-2 forcats_0.5.1
## [17] stringr_1.4.0 dplyr_1.0.4 purrr_0.3.4 readr_1.4.0
## [21] tidyr_1.1.2 tibble_3.0.6 tidyverse_1.3.0 MASS_7.3-53
## [25] plotly_4.9.3 ggplot2_3.3.3 redres_0.0.0.9
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_2.0-0 ellipsis_0.3.1
## [4] rio_0.5.26 fs_1.5.0 rstudioapi_0.13
## [7] farver_2.1.0 waldo_0.2.5 fansi_0.4.2
## [10] lubridate_1.7.10 xml2_1.3.2 splines_4.0.4
## [13] robustbase_0.93-7 Formula_1.2-4 jsonlite_1.7.2
## [16] nloptr_1.2.2.2 broom_0.7.5 dbplyr_2.1.0
## [19] shiny_1.6.0 compiler_4.0.4 httr_1.4.2
## [22] backports_1.2.1 assertthat_0.2.1 fastmap_1.1.0
## [25] lazyeval_0.2.2 cli_2.3.1 later_1.1.0.1
## [28] htmltools_0.5.1.1 tools_4.0.4 gtable_0.3.0
## [31] glue_1.4.2 Rcpp_1.0.6 cellranger_1.1.0
## [34] jquerylib_0.1.3 vctrs_0.3.6 nlme_3.1-152
## [37] crosstalk_1.1.1 xfun_0.21 testthat_3.0.2
## [40] openxlsx_4.2.3 rvest_0.3.6 mime_0.10
## [43] lifecycle_1.0.0 statmod_1.4.35 DEoptimR_1.0-8
## [46] hms_1.0.0 promises_1.2.0.1 qqplotr_0.0.4
## [49] yaml_2.2.1 curl_4.3 sass_0.3.1
## [52] stringi_1.5.3 highr_0.8 checkmate_2.0.0
## [55] boot_1.3-26 zip_2.1.1 rlang_0.4.10
## [58] pkgconfig_2.0.3 evaluate_0.14 lattice_0.20-41
## [61] labeling_0.4.2 htmlwidgets_1.5.3 cowplot_1.1.1
## [64] tidyselect_1.1.0 magrittr_2.0.1 R6_2.5.0
## [67] generics_0.1.0 DBI_1.1.1 pillar_1.5.0
## [70] haven_2.3.1 foreign_0.8-81 withr_2.4.1
## [73] abind_1.4-5 modelr_0.1.8 crayon_1.4.1
## [76] utf8_1.1.4 rmarkdown_2.7 readxl_1.3.1
## [79] data.table_1.14.0 reprex_1.0.0 digest_0.6.27
## [82] xtable_1.8-4 httpuv_1.5.5 numDeriv_2016.8-1.1
## [85] munsell_0.5.0 viridisLite_0.3.0 bslib_0.2.4