##' Run Statistical Analysis of Monthly Background Checks of Gun Purchase
##'
##' @param debug Optional boolean switch to indicate whether interim data is displayed;
##' default is \sQuote{FALSE}
##'
##' @return A \code{data.frame} is returned, contained all different prepared columns.
##'
##' @author Gregor Aisch and Josh Keller wrote the R code; Dirk Eddelbuettel created
##' and maintains the package.
##' @seealso The NY Times article presenting this analsysi undertaken by this package is
##' at \url{http://www.nytimes.com/interactive/2015/12/10/us/gun-sales-terrorism-obama-restrictions.html?}
##'
##' @examples
##' \dontrun{
##' gs <- analysis()
##' plot_gunsales(gs)
##' ggplot_gunsales(gs)
##' }
analysis <- function(debug=FALSE) {
## estimate gun sales using formula by Jurgen Brauer, published here
## http://www.smallarmssurvey.org/fileadmin/docs/F-Working-papers/SAS-WP14-US-Firearms-Industry.pdf
##
## note: the column `multiple_corrected` is a copy of `multiple` in which
## we set the checks in the "multiple" category to 0 for California
alldata <- alldata %>% mutate(guns_sold=(handgun + longgun) * 1.1 + multiple_corrected * 2)
#alldata <- mutate(alldata, guns_sold=(handgun + longgun) * 1.1 + multiple_corrected * 2)
## let's look at the total numbers; state_ts() is a helper function
total <- alldata %>% state_ts('Totals', 'guns_sold')
#total <- state_ts(all, 'Totals', 'guns_sold')
## compute seasonally adjusted gun sales (using final() and seas() from seasonal)
totalSeas <- total %>% seas %>% final
#totalSeas <- final(seas(total))
poptotal <- poptotal %>%
filter(year >= 2000) %>%
#filter(year < 2015 | month <= 11) %>%
with(ts(res_pop, start=c(2000,1), frequency = 12))
#poptotal <- ts(select(filter(poptotal, year >= 2000), "res_pop"),
# start=c(2000,1),
# freqency=12)
## normalize gun sales by population
totalSeasPop <- totalSeas / poptotal * 1000
totalSeasScaled <- totalSeas / 280726
## create a new data frame that eventually stores all the
## data we need in the final piece
out_data <- ts_to_dataframe(total, 'guns_total') %>%
mutate(guns_total=round(guns_total, 3))
## expand the data.frame, adding more volumns
out_data <- data.frame(out_data,
guns_total_seas=as.matrix(totalSeas),
guns_total_per_1000=round(as.matrix(totalSeasPop), digits=3),
guns_total_per_1000_scaled=round(as.matrix(totalSeasScaled), digits=3))
if (debug) {
print(head(out_data))
print(tail(out_data))
}
## create a temporary matrix for computing the
## handgun_share and longgun_share columns
## cbind works correctly here as it operates on timeseries object
tmp <- cbind(final(seas(state_ts(alldata, 'Totals', 'handgun'))),
final(seas(state_ts(alldata, 'Totals', 'longgun'))),
final(seas(state_ts(alldata, 'Totals', 'other'))),
final(seas(state_ts(alldata, 'Totals', 'multiple_corrected'))))
colnames(tmp) <- c('handgun', 'longgun', 'other', 'multiple')
out_data <- data.frame(out_data, tmp)
## convert NAs to 0 in column other
out_data$other[is.na(out_data$other)] <- 0
## compute the handgun/longgun share
out_data <- within(out_data, {
handgun_share=round(handgun / (handgun+longgun+other+multiple*0.5), 4)
longgun_share=round(longgun / (handgun+longgun+other+multiple*0.5), 4)
})
## plot percent of national for selected states
show_states <- c('New Jersey', 'Maryland', 'Georgia',
'Louisiana', 'Mississippi', 'Missouri')
for (s in show_states) {
s.ts <- state_data(alldata, s, total, totalSeas)
## merge with out_data
temp <- mutate(ts_to_dataframe(s.ts), value=round(value,3))
colnames(temp) <- c('year', 'month', gsub(' ', '_', tolower(s)))
out_data <- data.frame(out_data, temp[,3,drop=FALSE])
}
if (debug) {
print(head(out_data))
print(tail(out_data))
}
## compute handgun sales for DC: handung * 1.1 + multiple
dchandgun_checks <- state_ts(alldata, 'District of Columbia', 'handgun', outer_zeros_to_na=F)
dcmultiple <- state_ts(alldata, 'District of Columbia', 'multiple', outer_zeros_to_na=F)
dchandgun <- (dchandgun_checks * 1.1 + dcmultiple + 1) %>% seas %>% final - 1
totalHandgun <- (state_ts(alldata, 'Totals', 'handgun') * 1.1 +
state_ts(alldata, 'Totals', 'multiple')) %>% seas %>% final
dchandgunPct <- dchandgun / totalHandgun * 100000
## merge with out_data
temp <- ts_to_dataframe(round(dchandgunPct, 1), 'dc_handguns_per_100k_national_sales')
out_data <- data.frame(out_data, temp[,3,drop=FALSE])
## estimate how much more guns are sold missouri after law change
missouri <- state_data(alldata, 'Missouri', normalize = F, adj_seasonal = F)
missouri.avg_pre_2007 <- mean(missouri[73:84])
missouri.avg_post_2008 <- mean(missouri[97:108])
print(paste('Increase in monthly gun sales in Missouri =', round(missouri.avg_post_2008 - missouri.avg_pre_2007, digits=2)))
invisible(out_data)
}
## R-devel CMD check still whines about these:
Date <- dc_handguns_per_100k_national_sales <- guns_total <- guns_total_per_1000 <- NULL
guns_total_per_1000_scaled <- guns_total_seas <- handgun <- handgun_share <- longgun <- NULL
longgun_share <- month.num <- multiple <- multiple_corrected <- other <- res_pop <- NULL
state <- value <- year <- NULL
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.