Ricardo Lowe Course: Mortality
Performing an APC model using Uniform Crime Report from the FBI (1976-2018)
library(readxl)
library(knitr)
library(kableExtra)
library(readxl)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazer#IPUMS_NHGIS 2010
ustract <- read.csv(file="~/Desktop/Data Center/IPUMS_ NHGIS/Population2010/nhgis0029_ds172_2010_tract.csv", header=TRUE, sep=",")
usplce <- read.csv(file="~/Desktop/Data Center/IPUMS_ NHGIS/Population2010/nhgis0029_ds172_2010_place.csv", header=TRUE, sep=",")
uscnty <- read.csv(file="~/Desktop/Data Center/IPUMS_ NHGIS/Population2010/nhgis0029_ds172_2010_county.csv", header=TRUE, sep=",")
bses18 <- read.csv(file="~/Desktop/Data Center/IPUMS_ NHGIS/bses/nhgis0039_ds240_20185_2018_tract.csv", header=TRUE, sep=",")
beduc18 <- read.csv(file="~/Desktop/Data Center/IPUMS_ NHGIS/nhgis0041_ds240_20185_2018_tract.csv", header=TRUE, sep=",")
#UCR Data
shr_1976_2018 <- read.csv(file="~/Desktop/Data Center/UCR_ICPSR/100699-V8/shr_1976_2018_csv/shr_1976_2018.csv", header=TRUE, sep=",")#Identifying Blck Homicides in UCR data
shr_1976_2018$black_vhomicides <- ifelse(shr_1976_2018$victim_1_race == "black", 1, 0)
shr_1976_2018$black_ohomicides <- ifelse(shr_1976_2018$offender_1_race == "black", 1, 0)
#Categorizing Age in UCR data to create Cohorts
#Victim
shr_1976_2018$victim_age <- as.numeric(shr_1976_2018$victim_1_age)
shr_1976_2018$victim_age_cohort <- ifelse(shr_1976_2018$victim_age<=4, "00,04",
ifelse(shr_1976_2018$victim_age>= 5 & shr_1976_2018$victim_age<= 9, "05,09",
ifelse(shr_1976_2018$victim_age>=10 & shr_1976_2018$victim_age<=14, "10,14",
ifelse(shr_1976_2018$victim_age>=15 & shr_1976_2018$victim_age<=19, "15,19",
ifelse(shr_1976_2018$victim_age>=20 & shr_1976_2018$victim_age<=24, "20,24",
ifelse(shr_1976_2018$victim_age>=25 & shr_1976_2018$victim_age<=29, "25,29",
ifelse(shr_1976_2018$victim_age>=30 & shr_1976_2018$victim_age<=34, "30,34",
ifelse(shr_1976_2018$victim_age>=35 & shr_1976_2018$victim_age<=39, "35,39",
ifelse(shr_1976_2018$victim_age>=40 & shr_1976_2018$victim_age<=44, "40,44",
ifelse(shr_1976_2018$victim_age>=45 & shr_1976_2018$victim_age<=49, "45,49",
ifelse(shr_1976_2018$victim_age>=50 & shr_1976_2018$victim_age<=54, "50,54",
ifelse(shr_1976_2018$victim_age>=55 & shr_1976_2018$victim_age<=59, "55,59",
ifelse(shr_1976_2018$victim_age>=60 & shr_1976_2018$victim_age<=64, "60,64",
ifelse(shr_1976_2018$victim_age>=65 & shr_1976_2018$victim_age<=69, "65,69",
ifelse(shr_1976_2018$victim_age>=70 & shr_1976_2018$victim_age<=74, "70,74",
ifelse(shr_1976_2018$victim_age>=75 & shr_1976_2018$victim_age<=79, "75,79",
ifelse(shr_1976_2018$victim_age>=80 & shr_1976_2018$victim_age<=84, "80,84",
ifelse(shr_1976_2018$victim_age>=85, "Over 85", ""))))))))))))))))))
#Offender
shr_1976_2018$offender_age <- as.numeric(shr_1976_2018$offender_1_age)
shr_1976_2018$offender_age_cohort <- ifelse(shr_1976_2018$offender_age<=4, "00,04",
ifelse(shr_1976_2018$offender_age>= 5 & shr_1976_2018$offender_age<= 9, "05,09",
ifelse(shr_1976_2018$offender_age>=10 & shr_1976_2018$offender_age<=14, "10,14",
ifelse(shr_1976_2018$offender_age>=15 & shr_1976_2018$offender_age<=19, "15,19",
ifelse(shr_1976_2018$offender_age>=20 & shr_1976_2018$offender_age<=24, "20,24",
ifelse(shr_1976_2018$offender_age>=25 & shr_1976_2018$offender_age<=29, "25,29",
ifelse(shr_1976_2018$offender_age>=30 & shr_1976_2018$offender_age<=34, "30,34",
ifelse(shr_1976_2018$offender_age>=35 & shr_1976_2018$offender_age<=39, "35,39",
ifelse(shr_1976_2018$offender_age>=40 & shr_1976_2018$offender_age<=44, "40,44",
ifelse(shr_1976_2018$offender_age>=45 & shr_1976_2018$offender_age<=49, "45,49",
ifelse(shr_1976_2018$offender_age>=50 & shr_1976_2018$offender_age<=54, "50,54",
ifelse(shr_1976_2018$offender_age>=55 & shr_1976_2018$offender_age<=59, "55,59",
ifelse(shr_1976_2018$offender_age>=60 & shr_1976_2018$offender_age<=64, "60,64",
ifelse(shr_1976_2018$offender_age>=65 & shr_1976_2018$offender_age<=69, "65,69",
ifelse(shr_1976_2018$offender_age>=70 & shr_1976_2018$offender_age<=74, "70,74",
ifelse(shr_1976_2018$offender_age>=75 & shr_1976_2018$offender_age<=79, "75,79",
ifelse(shr_1976_2018$offender_age>=80 & shr_1976_2018$offender_age<=84, "80,84",
ifelse(shr_1976_2018$offender_age>=85, "Over 85", ""))))))))))))))))))
#Categorizing Periods in UCR data
shr_1976_2018$period <- as.numeric(shr_1976_2018$year)
shr_1976_2018$period_cohort <-ifelse(shr_1976_2018$period>=1976 & shr_1976_2018$period<=1980, "1976 to 1980",
ifelse(shr_1976_2018$period>=1981 & shr_1976_2018$period<=1985, "1981 to 1985",
ifelse(shr_1976_2018$period>=1986 & shr_1976_2018$period<=1990, "1986 to 1990",
ifelse(shr_1976_2018$period>=1991 & shr_1976_2018$period<=1995, "1991 to 1995",
ifelse(shr_1976_2018$period>=1996 & shr_1976_2018$period<=2000, "1996 to 2000",
ifelse(shr_1976_2018$period>=2001 & shr_1976_2018$period<=2005, "2001 to 2005",
ifelse(shr_1976_2018$period>=2006 & shr_1976_2018$period<=2010, "2006 to 2010",
ifelse(shr_1976_2018$period>=2011 & shr_1976_2018$period<=2015, "2011 to 2015",
ifelse(shr_1976_2018$period>=2016 & shr_1976_2018$period<=2018, "2016 to 2018",
"NA")))))))))
#Categorizing Births in UCR data
#Victim
shr_1976_2018$victim_birth_cohort <- ifelse(shr_1976_2018$year == 1976, (1976- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1977, (1977- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1978, (1978- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1979, (1979- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1980, (1980- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1981, (1981- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1982, (1982- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1983, (1983- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1984, (1984- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1985, (1985- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1986, (1986- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1987, (1987- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1988, (1988- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1989, (1989- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1990, (1990- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1991, (1991- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1992, (1992- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1993, (1993- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1994, (1994- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1995, (1995- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1996, (1996- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1997, (1997- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1998, (1998- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 1999, (1999- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2000, (2000- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2001, (2001- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2002, (2002- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2003, (2003- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2004, (2004- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2005, (2005- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2006, (2006- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2007, (2007- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2008, (2008- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2009, (2009- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2010, (2010- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2011, (2011- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2012, (2012- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2013, (2013- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2014, (2014- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2015, (2015- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2016, (2016- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2017, (2017- shr_1976_2018$victim_age),
ifelse(shr_1976_2018$year == 2018, (2018- shr_1976_2018$victim_age), " ")))))))))))))))))))))))))))))))))))))))))))
shr_1976_2018$cohort_test<-as.numeric(shr_1976_2018$victim_birth_cohort)
shr_1976_2018$victim_cohort<- cut(shr_1976_2018$cohort_test, breaks=seq(1865, 2018, 5))
#Offender
shr_1976_2018$offender_birth_cohort <-ifelse(shr_1976_2018$year == 1976, (1976- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1977, (1977- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1978, (1978- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1979, (1979- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1980, (1980- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1981, (1981- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1982, (1982- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1983, (1983- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1984, (1984- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1985, (1985- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1986, (1986- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1987, (1987- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1988, (1988- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1989, (1989- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1990, (1990- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1991, (1991- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1992, (1992- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1993, (1993- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1994, (1994- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1995, (1995- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1996, (1996- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1997, (1997- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1998, (1998- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 1999, (1999- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2000, (2000- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2001, (2001- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2002, (2002- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2003, (2003- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2004, (2004- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2005, (2005- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2006, (2006- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2007, (2007- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2008, (2008- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2009, (2009- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2010, (2010- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2011, (2011- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2012, (2012- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2013, (2013- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2014, (2014- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2015, (2015- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2016, (2016- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2017, (2017- shr_1976_2018$offender_age),
ifelse(shr_1976_2018$year == 2018, (2018- shr_1976_2018$offender_age), " ")))))))))))))))))))))))))))))))))))))))))))
shr_1976_2018$cohort_testo<-as.numeric(shr_1976_2018$offender_birth_cohort)
shr_1976_2018$offender_cohort<- cut(shr_1976_2018$cohort_testo, breaks=seq(1865, 2018, 5))APC Model
#https://rpubs.com/corey_sparks/436812
#http://www.r-inla.org/download
library(haven)
library(car)
require(INLA)## Loading required package: INLA## Loading required package: Matrix## Loading required package: sp## Loading required package: parallel## This is INLA_19.09.03 built 2019-09-03 09:07:31 UTC.
## See www.r-inla.org/contact-us for how to get help.
## To enable PARDISO sparse library; see inla.pardiso()#install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
#options(repos = c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/testing"))
#update.packages("INLA", dep=TRUE)
library(INLA)
library(dplyr)
library(ggplot2)library(ggplot2)
shr_1976_2018%>%
filter(victim_age<85)%>%
group_by(victim_cohort)%>%
summarise(meanvhomi=mean(black_vhomicides, na.rm=T))%>%
ggplot(aes(y=meanvhomi, x=victim_cohort))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk (Victims)", subtitle = "UCR by FBI")+theme(axis.text.x = element_text(angle = 90, hjust = 1))## Warning: Factor `victim_cohort` contains implicit NA, consider using
## `forcats::fct_explicit_na`
shr_1976_2018%>%
filter(offender_age<85)%>%
group_by(offender_cohort)%>%
summarise(meanohomi=mean(black_ohomicides, na.rm=T))%>%
ggplot(aes(y=meanohomi, x=offender_cohort))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk (Offenders)", subtitle = "UCR by FBI")+theme(axis.text.x = element_text(angle = 90, hjust = 1))## Warning: Factor `offender_cohort` contains implicit NA, consider using
## `forcats::fct_explicit_na`
shr_1976_2018%>%
filter(victim_age<85)%>%
group_by(year)%>%
summarise(meanvhomi=mean(black_vhomicides, na.rm=T))%>%
ggplot(aes(y=meanvhomi, x=year))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk by year (Victims)", subtitle = "UCR by FBI")
shr_1976_2018%>%
filter(offender_age<85)%>%
group_by(year)%>%
summarise(meanohomi=mean(black_ohomicides, na.rm=T))%>%
ggplot(aes(y=meanohomi, x=year))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk by year (Offenders)", subtitle = "UCR by FBI")
shr_1976_2018%>%
filter(victim_age<85)%>%
group_by(victim_age_cohort)%>%
summarise(meanvhomi=mean(black_vhomicides, na.rm=T))%>%
ggplot(aes(y=meanvhomi, x=victim_age_cohort))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk by age (Victims)", subtitle = "UCR by FBI")
shr_1976_2018%>%
filter(offender_age<85)%>%
group_by(offender_age_cohort)%>%
summarise(meanohomi=mean(black_ohomicides, na.rm=T))%>%
ggplot(aes(y=meanohomi, x=offender_age_cohort))+geom_line(aes(group=1))+ggtitle(label = "Mean Black Homicide Risk by age (Offenders)", subtitle = "UCR by FBI")
shr_1976_2018%>%
filter(victim_age<85)%>%
group_by(year, victim_cohort)%>%
summarise(meanvhomi=mean(black_vhomicides, na.rm=T))%>%
ggplot(aes(y=meanvhomi, x=year))+geom_line(aes(group=victim_cohort, color=victim_cohort))+ggtitle(label = "Mean Black Homicide Risk by year and cohort (Victims)", subtitle = "UCR by FBI")## Warning: Factor `victim_cohort` contains implicit NA, consider using
## `forcats::fct_explicit_na`
shr_1976_2018%>%
filter(offender_age<85)%>%
group_by(year, offender_cohort)%>%
summarise(meanohomi=mean(black_ohomicides, na.rm=T))%>%
ggplot(aes(y=meanohomi, x=year))+geom_line(aes(group=offender_cohort, color=offender_cohort))+ggtitle(label = "Mean Black Homicide Risk by year and cohort (Offenders)", subtitle = "UCR by FBI")## Warning: Factor `offender_cohort` contains implicit NA, consider using
## `forcats::fct_explicit_na`
#Creating a JOIN variable to merge US Cen data with UCR data (Pasting ST FIPS and PL FIPS together in both datasets)
shr_1976_2018$plst <- as.numeric(paste(shr_1976_2018$fips_state_code,shr_1976_2018$fips_place_code,sep = ""))## Warning: NAs introduced by coercionusplce$plst <- as.numeric(paste(usplce$STATEA,usplce$PLACEA,sep = ""))
shr_1976_2018$bhomis<- scale(shr_1976_2018$black_vhomicides)
UCR_ay = merge(shr_1976_2018[,c("plst","year","black_vhomicides","victim_1_sex","bhomis","victim_cohort","victim_age","victim_age_cohort")], usplce[,c("plst","PLACE","STATE","PLACEA","GISJOIN","STATE")])
library(lme4)## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car#user lmer to find good starting values for INLA
fit_in1.mer0<- lmer(bhomis~ victim_1_sex+(1|year)+(1|victim_cohort), data=UCR_ay)
#INLA parameterizes variances on the log-precision scale, so if we take log(1/variance), we can get that number and feed it to INLA
vc0<-as.data.frame(VarCorr(fit_in1.mer0))
vc0$init<-log(1/vc0$vcov); vc0## grp var1 var2 vcov sdcor init
## 1 year (Intercept) <NA> 0.004511923 0.06717085 5.40103177
## 2 victim_cohort (Intercept) <NA> 0.040933160 0.20231945 3.19581480
## 3 Residual <NA> <NA> 0.968884698 0.98431941 0.03160966#Feed the init values to the hyperparameter inits in each f() #IID random effects for period and cohort fit_in0<- inla(bhomis~ victim_1_sex+f(year, model=“iid”,hyper = list(prec=list(initial=vc0\(init[2])) )+f(victim_cohort, model="iid",hyper = list(prec=list(initial=vc0\)init[1]))), family = “gaussian”, data=UCR_ay, num.threads = 2, verbose = F, control.family = list(hyper = list(prec=list(initial=vc0$init[3]))), control.compute = list(waic=T))
cohort as smooth RW1 random effect
fit_in2<- inla(bhomis~ victim_1_sex+f(year, model=“iid”,hyper = list(prec=list(initial=vc0\(init[2])) )+f(victim_cohort, model="rw1",hyper = list(prec=list(initial=vc0\)init[1]))), family = “gaussian”, data=UCR_ay, num.threads = 2, verbose = F, control.family = list(hyper = list(prec=list(initial=vc0$init[3]))), control.compute = list(waic=T)) gc()
summary(fit_in0)
```