#### CONTEXTUAL INFLUENCES ON THE SENTENCING OF INDIVIDUALS CONVICTED OF SEXUAL CRIMES #### Replication script #### Setup Workspace #### install.packages(c("tidyverse", "haven", "magrittr", "Amelia", "MatchIt", "brms", "car", "lme4", "gridExtra")) library(tidyverse) # for data manipulation library(haven) # for reading in files library(Amelia) # for missing data library(MatchIt) # for counterfactuals library(magrittr) # for more data manipulation library(brms) # for model estimation library(sjstats) # for post-estimation quantities library(gridExtra) # for graphics #### Read in Data #### d = read_csv("https://www.dropbox.com/s/y777vc4d1mhhwg9/SO%20Sentencing%20Context%20Replication%20Data.csv?dl=1") # Relevel Factor Variables d$incarcerate = factor(d$incarcerate, levels = c('Non-Custodial Sanction', 'Incarcerated')) d$sex_offender = factor(d$sex_offender, levels = c('Matched Violent Offender', 'Sexual Offender')) d$female = factor(d$female, levels = c('Male', 'Female')) d$race_eth = factor(d$race_eth, levels = c('White', 'Black', 'Hispanic')) d$guideline_edition = factor(d$guideline_edition, levels = c('6th Ed', '5th Ed', '6th Ed, Revised', '7th Ed')) d$presump_inc = factor(d$presump_inc, levels = c('Presumptive Non-Custodial', 'Presumptive Incarceration')) d$mand_minimum = factor(d$mand_minimum, levels = c('No Mandatory Minimum', 'Mandatory Minimum Eligible')) d$disposition = factor(d$disposition, levels = c('Plea', 'Trial')) d$vic_under_13 = factor(d$vic_under_13, levels = c('Victim 13 or Older', 'Victim Under 13')) d$multi_convict = factor(d$multi_convict, levels = c('Single Conviction Offense', 'Multiple Conviction Offenses')) d$recidivist = factor(d$recidivist, levels = c('Non-Recidivist Offender', 'Recidivist Offender')) d$post_sandusky = factor(d$post_sandusky, levels = c('Pre-Sandusky Sentence', 'Post-Sandusky Sentence')) d$year = factor(d$year) d$court_size = factor(d$court_size, levels = c('Medium Court', 'Large Court', 'Small Court')) summary(d) jdist = read_csv("https://www.dropbox.com/s/e627hfrvnpamwpw/SO%20Sentencing%20Context%20JDistrict.csv?dl=1") summary(jdist) #### Missing Data Imputation #### # Check Missingness apply(d, 2, function(x) length(which(is.na(x)))) # Setup Amelia specification noms = c("sex_offender", "female", "race_eth", "guideline_edition", "presump_inc", "mand_minimum", "disposition", "vic_under_13", "multi_convict", "recidivist", "post_sandusky", "year", "court_size") ords = c("prs", "ogs") bds = matrix(c(6, 15, 100), nrow = 1, ncol = 3) ncpus = parallel::detectCores() m = 20 # Perform Imputation set.seed(938568) # for reproducibility d.am = amelia(d, m = m, noms = noms, ords = ords, bounds = bds, cs = "jdistrict", idvars = c("defendant_id", "incarcerate", "min_sentence", "LABEL", "OFN_LABEL", "OFN_GRADE"), empri = .01*nrow(d), parallel = "multicore", ncpus = ncpus) rm(noms, ords, bds, ncpus, m) # Selected at random - 6th imputation reported in manuscript d_complete = d.am$imputations[[6]] #### Precision Matching #### # Create integer distribution for age d_complete$age_f = factor(round(d_complete$age)) # Perform matching - Note: This will take some time set.seed(768947) # for reproducibility m_out = matchit(sex_offender == 'Sexual Offender' ~ factor(prs) + age_f + female + race_eth + recidivist, data = d_complete, method = "nearest", exact = c("prs", "age_f", "female", "race_eth", "recidivist")) summary(m_out) d_match = match.data(m_out) save(d_match, file = "Matched Data.RData") # Remove Judicial District Measures for Scaling d_match %<>% select(defendant_id, jdistrict, sex_offender, min_sentence, incarcerate, age, female, race_eth, prs, ogs, presump_inc, presump_sent, guideline_edition, mand_minimum, disposition, multi_convict, vic_under_13, recidivist, year, post_sandusky) #### Table 1 - Judicial District Descriptive Statistics #### # Means jdist %>% summarize_if(is.numeric, mean) # Standard Deviations jdist %>% summarize_if(is.numeric, sd) # % for Court Size prop.table(table(jdist$court_size)) * 100 #### Table 2 - Offender Level Descriptive Statistics #### # Overall Sample with(d_match, prop.table(table(incarcerate))) * 100 d_match %>% summarize(mean = mean(min_sentence), sd = sd(min_sentence)) d_match %>% summarize(mean = mean(age), sd = sd(age)) with(d_match, prop.table(table(female))) * 100 with(d_match, prop.table(table(race_eth))) * 100 d_match %>% summarize(mean = mean(prs), sd = sd(prs)) with(d_match, prop.table(table(recidivist))) * 100 d_match %>% summarize(mean = mean(ogs), sd = sd(ogs)) with(d_match, prop.table(table(vic_under_13))) * 100 with(d_match, prop.table(table(disposition))) * 100 with(d_match, prop.table(table(multi_convict))) * 100 with(d_match, prop.table(table(post_sandusky))) * 100 with(d_match, prop.table(table(presump_inc))) * 100 d_match %>% summarize(mean = mean(presump_sent), sd = sd(presump_sent)) with(d_match, prop.table(table(mand_minimum))) * 100 with(d_match, prop.table(table(guideline_edition))) * 100 # Sexual and Matched Violent Convictions options(pillar.sigfig = 4) # for decimals with(d_match, prop.table(table(sex_offender, incarcerate), 1)) * 100 d_match %>% group_by(sex_offender) %>% summarize(mean = mean(min_sentence), sd = sd(min_sentence)) d_match %>% group_by(sex_offender) %>% summarize(mean = mean(age), sd = sd(age)) with(d_match, prop.table(table(sex_offender, female), 1)) * 100 with(d_match, prop.table(table(sex_offender, race_eth), 1)) * 100 d_match %>% group_by(sex_offender) %>% summarize(mean = mean(prs), sd = sd(prs)) with(d_match, prop.table(table(sex_offender, recidivist), 1)) * 100 d_match %>% group_by(sex_offender) %>% summarize(mean = mean(ogs), sd = sd(ogs)) with(d_match, prop.table(table(sex_offender, vic_under_13), 1)) * 100 with(d_match, prop.table(table(sex_offender, disposition), 1)) * 100 with(d_match, prop.table(table(sex_offender, multi_convict), 1)) * 100 with(d_match, prop.table(table(sex_offender, post_sandusky), 1)) * 100 with(d_match, prop.table(table(sex_offender, presump_inc), 1)) * 100 d_match %>% group_by(sex_offender) %>% summarize(mean = mean(presump_sent), sd = sd(presump_sent)) with(d_match, prop.table(table(sex_offender, mand_minimum), 1)) * 100 with(d_match, prop.table(table(sex_offender, guideline_edition), 1)) * 100 # Scaling functions center_this = function(x){ (x - mean(as.numeric(x), na.rm = TRUE)) } scale_this = function(x){ (x - mean(as.numeric(x), na.rm = TRUE)) / sd(as.numeric(x), na.rm = TRUE) } jdist %<>% mutate(jail_capacityZ = scale_this(jail_capacity), court_caseloadZ = scale_this(court_caseload), pct_femaleZ = scale_this(pct_female), pct_under14Z = scale_this(pct_under14), political_competitionZ = scale_this(political_competition), religious_heteroZ = scale_this(religious_heterogeneity)) d_match = merge(d_match, jdist, by = "jdistrict", all.x = TRUE) d_match$court_size = factor(d_match$court_size, levels = c('Medium Court', 'Large Court', 'Small Court')) d_match %<>% group_by(sex_offender) %>% mutate(ageZ = scale_this(age), femaleC = center_this(female == 'Female'), blackC = center_this(race_eth == 'Black'), hispanicC = center_this(race_eth == 'Hispanic'), prsZ = scale_this(prs), ogsZ = scale_this(ogs), presump_sentZ = scale_this(presump_sent), presump_incC = center_this(presump_inc == 'Presumptive Incarceration'), recidivistC = center_this(recidivist == 'Recidivist Offender'), vic_under_13C = center_this(vic_under_13 == 'Victim Under 13'), trialC = center_this(disposition == 'Trial'), mand_minimumC = center_this(mand_minimum == 'Mandatory Minimum Eligible'), multi_convictC = center_this(multi_convict == 'Multiple Conviction Offenses'), post_sanduskyC = center_this(post_sandusky == 'Post-Sandusky Sentence'), ged_5thC = center_this(guideline_edition == '5th Ed'), ged_6thRC = center_this(guideline_edition == '6th Ed, Revised'), ged_7thC = center_this(guideline_edition == '7th Ed'), year_2005C = center_this(year == 2005), year_2006C = center_this(year == 2006), year_2007C = center_this(year == 2007), year_2008C = center_this(year == 2008), year_2009C = center_this(year == 2009), year_2010C = center_this(year == 2010), year_2011C = center_this(year == 2011), year_2012C = center_this(year == 2012), year_2013C = center_this(year == 2013), year_2014C = center_this(year == 2014), year_2015C = center_this(year == 2015)) #### Pre-Estimation Checks #### # Square Root of the VIF sqrt(car::vif(glm(min_sentence > 0 ~ sex_offender + ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + recidivistC + ogsZ + vic_under_13C + trialC + multi_convictC + post_sanduskyC + presump_sentZ + mand_minimumC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ +religious_heteroZ, data = d_match, family = binomial(link = "logit")))) # Null Model Intraclass Correlations library(lme4) so_inc_null = glmer(min_sentence > 0 ~ 1 + (1 | jdistrict), data = d_match, family = binomial(link = "logit"), subset = sex_offender == 'Sexual Offender') icc(so_inc_null) so_length_null = glmer.nb(round(min_sentence) ~ 1 + (1 | jdistrict), data = d_match, subset = sex_offender == 'Sexual Offender') icc(so_length_null) nso_inc_null = glmer(min_sentence > 0 ~ 1 + (1 | jdistrict), data = d_match, family = binomial(link = "logit"), subset = sex_offender == 'Matched Violent Offender') icc(nso_inc_null) nso_length_null = glmer.nb(round(min_sentence) ~ 1 + (1 | jdistrict), data = d_match, subset = sex_offender == 'Matched Violent Offender') icc(nso_length_null) #### Estimate Models #### # Set Weakly Informative Prior for Parameters hurdle_prior = c(prior(normal(0, 5), class = "b"), prior(student_t(3, 0, 1), class = "Intercept"), prior(student_t(3, 0, 1), class = "sd"), prior(normal(0, 5), class = "b", dpar = "hu"), prior(logistic(0, 1), class = "Intercept", dpar = "hu"), prior(student_t(3, 0, 1), class = "sd", dpar = "hu"), prior(gamma(0.01, 0.01), class = "shape")) # Sexual Offenders Model set.seed(200024) so_hurdle = brm(bf(round(min_sentence) ~ ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + recidivistC + ogsZ + vic_under_13C + trialC + multi_convictC + post_sanduskyC + presump_sentZ + mand_minimumC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ + (1 | jdistrict), hu ~ ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + ogsZ + vic_under_13C + trialC + multi_convictC + post_sanduskyC + presump_incC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ + (1 | jdistrict)), data = subset(d_match, sex_offender == 'Sexual Offender'), family = hurdle_negbinomial(), chains = 4, cores = 4, iter = 8000, warmup = 4000, thin = 8, control = list(adapt_delta = 0.80), prior = hurdle_prior, inits = 0) print(so_hurdle, prob = 0.95, digits = 3) save(so_hurdle, file = "SO Hurdle brms.RData") # Highest Density Interval hdi(so_hurdle, prob = c(.95, .99, .999)) # Matched Violent Offenders Model set.seed(716642) nso_hurdle = brm(bf(round(min_sentence) ~ ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + recidivistC + ogsZ + vic_under_13C + trialC + multi_convictC + post_sanduskyC + presump_sentZ + mand_minimumC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ + (1 | jdistrict), hu ~ ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + ogsZ + vic_under_13C + trialC + multi_convictC + post_sanduskyC + presump_incC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ + (1 | jdistrict)), data = subset(d_match, sex_offender == 'Matched Violent Offender'), family = hurdle_negbinomial(), chains = 4, cores = 4, iter = 8000, warmup = 4000, thin = 8, control = list(adapt_delta = 0.80), prior = hurdle_prior, inits = 0) print(nso_hurdle, prob = 0.95, digits = 3) save(nso_hurdle, file = "NSO Hurdle brms.RData") # Highest Density Intervals hdi(nso_hurdle, prob = c(.95, .99, .999)) # Interaction Model to Determine Group Differences set.seed(303069) int_hurdle = brm(bf(round(min_sentence) ~ (ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + recidivistC + ogsZ + vic_under_13C + trialC + multi_convictC + post_sandusky + presump_sentZ + mand_minimumC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ) * sex_offender + (1 | jdistrict), hu ~ (ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + ogsZ + vic_under_13C + trialC + multi_convictC + post_sandusky + presump_incC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ) * sex_offender + (1 | jdistrict)), data = d_match, family = hurdle_negbinomial(), chains = 4, cores = 4, iter = 8000, warmup = 4000, thin = 8, control = list(adapt_delta = 0.80), prior = hurdle_prior, inits = 0) summary(int_hurdle, prob = 0.95, digits = 6) save(int_hurdle, file = "Int Hurdle brms.RData") # Highest Density Interval hdi(int_hurdle, prob = c(.95, .99, .999)) #### Marginal Effects #### # Need to estimate separate logit model to extract incarceration probabilities d_match$incarcerate = as.numeric(d_match$incarcerate) - 1 d_match$trials = 1 logit_prior = c(prior(normal(0, 5), class = "b"), prior(student_t(3, 0, 1), class = "Intercept"), prior(student_t(3, 0, 1), class = "sd")) set.seed(303069) # Same as hurdle above int_logit = brm(incarcerate | trials(trials) ~ (ageZ + I(ageZ^2) + femaleC + blackC + hispanicC + prsZ + ogsZ + vic_under_13C + trialC + multi_convictC + post_sandusky + presump_incC + ged_5thC + ged_6thRC + ged_7thC + year_2005C + year_2006C + year_2007C + year_2008C + year_2009C + year_2010C + year_2011C + year_2012C + year_2013C + year_2014C + year_2015C + jail_capacityZ + court_caseloadZ + court_size + pct_femaleZ + pct_under14Z + political_competitionZ + religious_heteroZ) * sex_offender + (1 | jdistrict), data = d_match, family = binomial(link = "logit"), chains = 4, cores = 4, iter = 8000, warmup = 4000, thin = 8, control = list(adapt_delta = 0.80), prior = logit_prior, inits = 0) summary(int_logit, prob = 0.95, digits = 6) save(int_logit, file = "Int Logit brms.RData") ic = list(jail_capacityZ = seq(-2, 2, 0.25)) jail_cap = marginal_effects(int_logit, effects = "sex_offender:jail_capacityZ", int_conditions = ic) jci = ggplot(jail_cap[[1]], aes(x = parse_number(jail_capacityZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Jail Capacity (Z Score)") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) jci jail_cap = marginal_effects(int_hurdle, effects = "sex_offender:jail_capacityZ", int_conditions = ic) jcl = ggplot(jail_cap[[1]], aes(x = parse_number(jail_capacityZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Jail Capacity (Z Score)") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) jcl ic = list(court_caseloadZ = seq(-2, 2, 0.25)) caseload = marginal_effects(int_logit, effects = "sex_offender:court_caseloadZ", int_conditions = ic) cli = ggplot(caseload[[1]], aes(x = parse_number(court_caseloadZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Court Caseload (Z Score)") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) cli caseload = marginal_effects(int_hurdle, effects = "sex_offender:court_caseloadZ", int_conditions = ic) cll = ggplot(caseload[[1]], aes(x = parse_number(court_caseloadZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Court Caseload (Z Score)") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) cll court_size = marginal_effects(int_logit, effects = "sex_offender:court_size")[[1]] court_size$court_size = factor(court_size$court_size, levels = c("Large Court", "Medium Court", "Small Court")) csi = ggplot(court_size, aes(x = court_size, y = estimate__, group = sex_offender)) + geom_pointrange(aes(ymin = lower__, ymax = upper__, shape = sex_offender), color = "black", position = position_dodge(width = 0.5), size = 0.5) + xlab("Court Size") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_shape_manual(name = NULL, values = c(16, 17), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) csi court_size = marginal_effects(int_hurdle, effects = "sex_offender:court_size")[[1]] court_size$court_size = factor(court_size$court_size, levels = c("Large Court", "Medium Court", "Small Court")) csl = ggplot(court_size, aes(x = court_size, y = estimate__, group = sex_offender)) + geom_pointrange(aes(ymin = lower__, ymax = upper__, shape = sex_offender), color = "black", position = position_dodge(width = 0.5), size = 0.5) + xlab("Court Size") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_shape_manual(name = NULL, values = c(16, 17), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) csl ic = list(political_competitionZ = seq(-2, 2, 0.25)) polcomp = marginal_effects(int_logit, effects = "sex_offender:political_competitionZ", int_conditions = ic) pci = ggplot(polcomp[[1]], aes(x = parse_number(political_competitionZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Political Competition (Z Score)") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) pci polcomp = marginal_effects(int_hurdle, effects = "sex_offender:political_competitionZ", int_conditions = ic) pcl = ggplot(polcomp[[1]], aes(x = parse_number(political_competitionZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Political Competition (Z Score)") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) pcl ic = list(religious_heteroZ = seq(-2, 2, 0.25)) relhet = marginal_effects(int_logit, effects = "sex_offender:religious_heteroZ", int_conditions = ic) rhi = ggplot(relhet[[1]], aes(x = parse_number(religious_heteroZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Religious Heterogeneity (Z Score)") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) rhi relhet = marginal_effects(int_hurdle, effects = "sex_offender:religious_heteroZ", int_conditions = ic) rhl = ggplot(relhet[[1]], aes(x = parse_number(religious_heteroZ), y = estimate__, group = sex_offender)) + geom_line(aes(linetype = sex_offender), size = 1.0, color = "black") + geom_ribbon(aes(ymin = lower__, ymax = upper__, fill = sex_offender), alpha = 0.3, color = NA) + xlab("Religious Heterogeneity (Z Score)") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_linetype_manual(name = NULL, values = c(3, 1), guide = FALSE) + scale_fill_manual(name = NULL, values = c("#bdbdbd", "#636363"), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) rhl post_sandusky = marginal_effects(int_logit, effects = "sex_offender:post_sandusky")[[1]] post_sandusky$post_sandusky = recode_factor(post_sandusky$post_sandusky, `Pre-Sandusky Sentence` = "Pre-Sandusky", `Post-Sandusky Sentence` = "Post-Sandusky") psi = ggplot(post_sandusky, aes(x = post_sandusky, y = estimate__, group = sex_offender)) + geom_pointrange(aes(ymin = lower__, ymax = upper__, shape = sex_offender), color = "black", position = position_dodge(width = 0.5), size = 0.5) + xlab("Sentencing Relative to Sandusky Case") + ylab("Probabilty of Incarceration") + scale_y_continuous(limits = c(0.70, 1)) + scale_shape_manual(name = NULL, values = c(16, 17), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) psi post_sandusky = marginal_effects(int_hurdle, effects = "sex_offender:post_sandusky")[[1]] post_sandusky$post_sandusky = recode_factor(post_sandusky$post_sandusky, `Pre-Sandusky Sentence` = "Pre-Sandusky", `Post-Sandusky Sentence` = "Post-Sandusky") psl = ggplot(post_sandusky, aes(x = post_sandusky, y = estimate__, group = sex_offender)) + geom_pointrange(aes(ymin = lower__, ymax = upper__, shape = sex_offender), color = "black", position = position_dodge(width = 0.5), size = 0.5) + xlab("Sentencing Relative to Sandusky Case") + ylab("Sentence Length (Months)") + scale_y_continuous(limits = c(0, 40)) + scale_shape_manual(name = NULL, values = c(16, 17), guide = FALSE) + theme_classic() + theme(axis.text = element_text(color = "black")) psl # Arrange Grids p1 = grid.arrange(jci, jcl, cli, cll, csi, csl, nrow = 3) p2 = grid.arrange(pci, pcl, rhi, rhl, psi, psl, nrow = 3) ggsave(file = "Sentencing Context Figure 1.png", p1, width = 6.5, height = 7, units = "in", dpi = 800) ggsave(file = "Sentencing Context Figure 2.png", p2, width = 6.5, height = 7, units = "in", dpi = 800)