library(DasGuptR)
library(tidyverse)
library(rsample)
data("reconv")
# seed required for bootstrapping
set.seed(as.numeric(as.Date("2026/05/21")))Standardization and decomposition with small differences in crude rates
Results from standardization and decomposition express change in crude rates as a percentage. When the absolute difference in crude rates is very small these percentages can become extremely large - and in the case where the crude rates are identical they are not defined. As such analysts may want to asses whether the difference in crude rates between comparators in a standardization and decomposition is ‘big enough’ to give reliable results. In the case where rates reflect aggregations of binary variables (as in the reconvictions data included with the package) one way to assess whether differences in crude rates are ‘big enough’ is via a statistical significance test for proportions. If the test is statistically significant then the results of the standardization and decomposition are likely to be robust to small variations in the input data. Conversely, if crude rates in comparator populations are not ‘statistically significant’ results from standardization and decomposition may be fragile to small changes in the input data. This vignette uses bootstrapping to demonstrate this phenomenon.
Load packages
Comparing 2004 and 2016
First we conduct a standardization and decomposition for the first and last years in the reconvictions data.
# original analysis
dat_2004_2016 <-
reconv |>
filter(year == "2004" | year == "2016")
res_2004_2016 <-
dgnpop(
dat_2004_2016,
pop="year",
factors=c("prev_rate"),id_vars=c("Sex","Age"),
crossclassified="offenders") |>
dg_table() |>
rownames_to_column()
res_2004_2016 rowname 2004 2016 diff decomp
1 Age_struct 0.31 0.28 -0.02147 41.3
2 prev_rate 0.31 0.28 -0.02975 57.2
3 Sex_struct 0.29 0.29 -0.00076 1.5
4 crude 0.32 0.27 -0.05198 100.0This works as anticipated.
Uncertainty in comparison between 2004 and 2016
We may want to calculate uncertainty intervals around our results. To do this we bootstrap the underlying reconvictions data we find that the results from the 2004 to 2016 comparison are robust sensitive to perturbations of the input data.
# bootstrap helper functions
uncount_split <- function(splits){
rsample::analysis(splits) |>
count(year, Age, Sex, reconvicted) |>
pivot_wider(names_from = "reconvicted",
values_from = "n") |>
mutate(offenders = n_not_reconvicted + reconvicted,
prev_rate = reconvicted / offenders)
}
calc_dg <- function(d){
d |>
dgnpop(pop="year",
factors=c("prev_rate"),id_vars=c("Sex","Age"),
crossclassified="offenders",
agg = TRUE) |>
dg_table() |>
rownames_to_column()
}
# bootstrap:
n_draws <- 500
dat_2004_2016_long <-
dat_2004_2016 |>
mutate(n_not_reconvicted = offenders - reconvicted) |>
select(year, Age, Sex, n_not_reconvicted, reconvicted) |>
pivot_longer(n_not_reconvicted:reconvicted,
names_to = "reconvicted",
values_to = "n") |>
mutate(n = as.integer(n)) |>
uncount(n)
draws_2004_2016 <-
dat_2004_2016_long |>
group_by(year) |>
rsample::bootstraps(times = n_draws) |>
mutate(dat = map(splits, uncount_split),
results = map(dat, calc_dg),
draw = 1:n_draws) |>
select(draw, results)Plotting the bootstrap draws against the observed data shows that the original results align with the median result of the bootstrap draws.
draws_2004_2016 |>
unnest(results) |>
pivot_longer(`2004`:`2016`,
names_to = "year",
values_to = "rate") |>
filter(year == "2016") |>
ggplot(aes(x = decomp, y = rowname)) +
ggdist::stat_halfeye() +
geom_point(data = res_2004_2016,
aes(x = decomp, y = rowname),
size = 4,colour = "red") +
ggtitle(label = "Red point is real data") +
facet_wrap(~ rowname, scales = "free_x")
draws_2004_2016 |>
unnest(results) |>
pivot_longer(`2004`:`2016`,
names_to = "year",
values_to = "rate") |>
filter(year == "2016",
rowname != "crude") |>
group_by(rowname, year) |>
summarise(confs = list(quantile(decomp, c(0.025, 0.5, 0.975))),
cis = list(c("lower_ci", "median", "upper_ci")),
.groups = "drop") |>
unnest(c(confs, cis)) |>
pivot_wider(names_from = cis, values_from = confs)# A tibble: 3 × 5
rowname year lower_ci median upper_ci
<chr> <chr> <dbl> <dbl> <dbl>
1 Age_struct 2016 36.66 41.38 47.41
2 Sex_struct 2016 0.9548 1.485 2.05
3 prev_rate 2016 51.14 57.15 62.18res_2004_2016 rowname 2004 2016 diff decomp
1 Age_struct 0.31 0.28 -0.02147 41.3
2 prev_rate 0.31 0.28 -0.02975 57.2
3 Sex_struct 0.29 0.29 -0.00076 1.5
4 crude 0.32 0.27 -0.05198 100.0However, the crude difference in reconviction rates in 2004 and 2016 is quite large (2016 rate is 0.84 the rate in 2004) conditions favourable to standardization and decomposition. Plotting the crude rates alongside their confidence intervals using prop.test() gives a sense of how different the crude rates in pairs of comparison years are whilst adjusting for the size of the reconviction cohorts.
crude_rates <-
reconv |>
group_by(year) |>
summarise(crude_rate = weighted.mean(prev_rate, offenders),
offenders = sum(offenders),
reconvicted = sum(reconvicted))
ref <-
crude_rates |>
filter(year == "2016")
diff_cis <-
crude_rates |>
filter_out(year == "2016") |>
mutate(
res = map2(
reconvicted, offenders,
~ prop.test(
c(..1, ref$reconvicted),
c(..2, ref$offenders)
) |> broom::tidy()
)
)
diff_cis |>
unnest(res) |>
mutate(diff = estimate1 - estimate2) |>
ggplot(aes(x = year, y = diff)) +
geom_hline(yintercept = 0) +
geom_point() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high)) +
labs(y = "Difference in crude reconviction rate compared to 2016",
x = "Reconviction cohort")
The figure above shows that confidence intervals for the difference in crude rates between 2004 and 2016 does not include zero - but the interval for the 2015-2016 differences does (and we can see that the p-value for this comparison is not significant at the .05 level):
diff_cis |>
unnest(res) |>
filter(p.value >= 0.05) |>
select(year, estimate1, estimate2, p.value)# A tibble: 1 × 4
year estimate1 estimate2 p.value
<int> <dbl> <dbl> <dbl>
1 2015 0.2725 0.2718 0.8196Indeed, the reconviction rates in 2015 and 2016 only differ at the 3rd decimal.
Comparing 2015 and 2016
If we attempt to do the same decomposition as we did for 2004 and 2016, but this time using 2015 and 2016, then Das Gupta’s maths appears to work as before, and the output tells us that 114% of the difference between these two years is due to changes in age structure:
dat_2015_2016 <-
reconv |>
filter(year == "2015" | year == "2016")
res_2015_2016 <-
dgnpop(
dat_2015_2016,
pop="year",
factors=c("prev_rate"),id_vars=c("Sex","Age"),
crossclassified="offenders") |>
dg_table() |>
rownames_to_column()
res_2015_2016 rowname 2015 2016 diff decomp
1 Age_struct 0.2725 0.2717 -8.246e-04 114.25
2 prev_rate 0.2721 0.2722 1.044e-04 -14.47
3 Sex_struct 0.2721 0.2721 -1.569e-06 0.22
4 crude 0.2725 0.2718 -7.218e-04 100.00However, we are decomposing a very small crude rate difference. These decomposition effects expressed as percentages are extremely unreliable. In short, this is because we are dividing by such a small value of -7.2^{-4}.
Uncertainty in comparison between 2015 and 2016
If we repeat the process for the 2015 to 2016 comparison, we can see that not only are the bootstrap distribution of estimates extremely wide, but the observed results do not align with the center of these distributions:
# bootstrapping 2015 & 2016
dat_2015_2016_long <-
dat_2015_2016 |>
mutate(n_not_reconvicted = offenders - reconvicted) |>
select(year, Age, Sex, n_not_reconvicted, reconvicted) |>
pivot_longer(n_not_reconvicted:reconvicted,
names_to = "reconvicted",
values_to = "n") |>
mutate(n = as.integer(n)) |>
uncount(n)
draws_2015_2016 <-
dat_2015_2016_long |>
group_by(year) |>
rsample::bootstraps(times = n_draws) |>
mutate(dat = map(splits, uncount_split),
results = map(dat, calc_dg),
draw = 1:n_draws) |>
select(draw, results)draws_2015_2016 |>
unnest(results) |>
pivot_longer(`2015`:`2016`,
names_to = "year",
values_to = "rate") |>
filter(year == "2016",
rowname != "crude") |>
ggplot(aes(x = decomp, y = rowname)) +
ggdist::stat_halfeye() +
geom_point(data = res_2015_2016 |> filter(rowname != "crude"),
aes(x = decomp, y = rowname),
size = 4, colour = "red"
) +
ggtitle(label = "Red point is real data") +
facet_wrap(~ rowname, scales = "free_x") +
scale_x_continuous(limits = c(-250, 250))
The red points (the observed results) do not align with the original results, as the results table confirms:
draws_2015_2016 |>
unnest(results) |>
pivot_longer(`2015`:`2016`,
names_to = "year",
values_to = "rate") |>
filter(year == "2016") |>
group_by(rowname, year) |>
summarise(confs = list(quantile(decomp, c(0.025, 0.5, 0.975))),
cis = list(c("lower_ci", "median", "upper_ci")),
.groups = "drop") |>
unnest(c(confs, cis)) |>
pivot_wider(names_from = cis, values_from = confs)# A tibble: 4 × 5
rowname year lower_ci median upper_ci
<chr> <chr> <dbl> <dbl> <dbl>
1 Age_struct 2016 -221.0 17.63 418.8
2 Sex_struct 2016 -41.82 0.365 32.14
3 crude 2016 100 100 100
4 prev_rate 2016 -303.2 83.02 327.9 Summary
Calculating appropriate confidence intervals for ratio data is a well studied problem in the statistics literature but can be challenging in applied settings1. The process described here - calculating a statistical significance test for the difference in crude rates and/or bootstrapping the underlying data before performing the decomposition - can identify situations where researchers may want to treat the outputs of a standardization and decomposition with caution.
Footnotes
See Franz, V. H. (2007). Ratios: A short guide to confidence limits and proper use. arXiv preprint arXiv:0710.2024.↩︎