Practical 3

Choose your own adventure

  1. Find some data relevant to your research interests and apply standardization and decomposition to them. Ask us for help if you get stuck!
  2. If you are still working on exercises from the previous practicals, feel free to go back and work on them
  3. Work through the bootstrapping examples below in your preferred language. There is a lot of detail here (particularly in the R example), so don’t worry if this feels like a lot. As a result we have provided all the code for you to run and suggestions for settings you might want to change to see how it affects the results.

Bootstrapping standardization and decomposition

These are the data we will be using. I don’t know anything about them, but this is the description from (Li 2017, p496)

[These data] compare desire for more children in two groups of women: those with 4 or more children (represented by parity 4+) and those with 1 child (represented as parity 1). Given that age is an important determinant of fertility and that most parity 1 women are likely younger than women with 4 or more children, the question is how to isolate the effect of age composition differences in the two parity groups. rdecompose can be applied to isolate the effect of age composition from actual differences in rates between the two groups of women

Table 2. Population Size and Percent Desiring More Children (Rate) by Age

Age group Parity 4+ Size (Nᵢ) Parity 4+ Rate (Tᵢ) Parity 1 Size (Nᵢ) Parity 1 Rate (Tᵢ)
20 to 24 27 37.037 363 90.083
25 to 29 152 19.079 208 76.923
30 to 34 224 15.179 96 56.250
35 to 39 239 5.021 59 20.339
40 to 44 211 6.161 48 10.417
All ages 853 11.489 774 72.093

Source: Clogg and Eliason (1988) from (Li 2017, p497)

R

Load packages and set the seed for the random number generator

library(tidyverse)
library(rsample)
library(DasGuptR)

set.seed(346)

Read in the data.

clogg_long <- read_csv("https://raw.githubusercontent.com/benmatthewsed/sgsss_std_decomp/refs/heads/main/data/example2-clogg-long.csv")

Running the original analysis

original_res <- 
DasGuptR::dgnpop(
  clogg_long,
  pop = "parity",
       factors = c("rate"),
      id_vars = "age_groups",
      crossclassified = "size") |> 
  DasGuptR::dg_table()

original_res
                         1        4      diff decomp
age_groups_struct 48.61896 25.54727 -23.07169  38.07
rate              55.84938 18.31685 -37.53253  61.93
crude             72.09319 11.48897 -60.60422 100.00

Exploding the data

Before resampling we need to turn the aggregate data into person-level data.

clogg_uncount <- 
clogg_long |> 
  uncount(size, .remove = FALSE) |> # we need to keep the size variable so set .remove as FALSE
  group_by(age_groups, parity) |> 
  mutate(
    d = as.integer(row_number() <= round(rate * size / 100))
  ) |> 
  ungroup()

clogg_uncount
# A tibble: 1,627 × 5
   age_groups parity  size  rate     d
   <chr>       <dbl> <dbl> <dbl> <int>
 1 20 to 24        4    27  37.0     1
 2 20 to 24        4    27  37.0     1
 3 20 to 24        4    27  37.0     1
 4 20 to 24        4    27  37.0     1
 5 20 to 24        4    27  37.0     1
 6 20 to 24        4    27  37.0     1
 7 20 to 24        4    27  37.0     1
 8 20 to 24        4    27  37.0     1
 9 20 to 24        4    27  37.0     1
10 20 to 24        4    27  37.0     1
# ℹ 1,617 more rows

Create the decomposition wrapper function

We need a shorthand way to perform the standardization and decomposition for each of our eventual bootstrap replications. We can create an R function to do this.

dg_clogg <- function(clogg_df){

DasGuptR::dgnpop(
  clogg_df,
  pop = "parity",
       factors = c("rate"),
      id_vars = "age_groups",
      crossclassified = "size") |> 
  DasGuptR::dg_table() |> 
  rownames_to_column() # this will be helpful later!

}

Re-aggregating the data

We also need a way to reaggregate each resampled dataset before the standardization and decomposition. We can create another function to do this. Note the use of rsample::analysis(splits) - this lets us access each bootstrap draw in a memory efficient way.

agg_clogg <- function(splits){
  
  rsample::analysis(splits) |> 
      count(age_groups, parity, d, name = "size") |> # this is the opposite of uncount()
  pivot_wider(names_from = "d", # turns the dataset wider to make the next step easier
                values_from = "size") |> 
  mutate(size = `1` + `0`,
  rate = `1` / size * 100) |> 
  select(age_groups, parity, size, rate) # keep only the columns we need
  
}

Resampling the data

Now we can bootstrap the data. Try changing the times argument from 100 to 200 and then to 1000 to see what happens.

clogg_straps <- 
rsample::bootstraps(
  clogg_uncount, times = 100
)

clogg_straps
# Bootstrap sampling 
# A tibble: 100 × 2
   splits             id          
   <list>             <chr>       
 1 <split [1627/598]> Bootstrap001
 2 <split [1627/594]> Bootstrap002
 3 <split [1627/612]> Bootstrap003
 4 <split [1627/611]> Bootstrap004
 5 <split [1627/596]> Bootstrap005
 6 <split [1627/614]> Bootstrap006
 7 <split [1627/599]> Bootstrap007
 8 <split [1627/604]> Bootstrap008
 9 <split [1627/590]> Bootstrap009
10 <split [1627/606]> Bootstrap010
# ℹ 90 more rows

Now we use the map function to reaggregate each bootstrap replication. This iterates over each bootstrap draw and applies our agg_clogg function to each one in turn, and then saves the results in a new column.

clogg_straps <- 
clogg_straps |> 
  mutate(clogg_agg = map(splits, agg_clogg))

clogg_straps
# Bootstrap sampling 
# A tibble: 100 × 3
   splits             id           clogg_agg        
   <list>             <chr>        <list>           
 1 <split [1627/598]> Bootstrap001 <tibble [10 × 4]>
 2 <split [1627/594]> Bootstrap002 <tibble [10 × 4]>
 3 <split [1627/612]> Bootstrap003 <tibble [10 × 4]>
 4 <split [1627/611]> Bootstrap004 <tibble [10 × 4]>
 5 <split [1627/596]> Bootstrap005 <tibble [10 × 4]>
 6 <split [1627/614]> Bootstrap006 <tibble [10 × 4]>
 7 <split [1627/599]> Bootstrap007 <tibble [10 × 4]>
 8 <split [1627/604]> Bootstrap008 <tibble [10 × 4]>
 9 <split [1627/590]> Bootstrap009 <tibble [10 × 4]>
10 <split [1627/606]> Bootstrap010 <tibble [10 × 4]>
# ℹ 90 more rows

Iterating the standardization and decomposition

With our aggregated bootstrap draws, we can use the map function to apply our custom dg_clogg function to run the standardization and decomposition to each draw.

clogg_straps <- 
clogg_straps |> 
  mutate(results = map(clogg_agg, dg_clogg))

clogg_straps
# Bootstrap sampling 
# A tibble: 100 × 4
   splits             id           clogg_agg         results     
   <list>             <chr>        <list>            <list>      
 1 <split [1627/598]> Bootstrap001 <tibble [10 × 4]> <df [3 × 5]>
 2 <split [1627/594]> Bootstrap002 <tibble [10 × 4]> <df [3 × 5]>
 3 <split [1627/612]> Bootstrap003 <tibble [10 × 4]> <df [3 × 5]>
 4 <split [1627/611]> Bootstrap004 <tibble [10 × 4]> <df [3 × 5]>
 5 <split [1627/596]> Bootstrap005 <tibble [10 × 4]> <df [3 × 5]>
 6 <split [1627/614]> Bootstrap006 <tibble [10 × 4]> <df [3 × 5]>
 7 <split [1627/599]> Bootstrap007 <tibble [10 × 4]> <df [3 × 5]>
 8 <split [1627/604]> Bootstrap008 <tibble [10 × 4]> <df [3 × 5]>
 9 <split [1627/590]> Bootstrap009 <tibble [10 × 4]> <df [3 × 5]>
10 <split [1627/606]> Bootstrap010 <tibble [10 × 4]> <df [3 × 5]>
# ℹ 90 more rows

We can now just select the results and the bootstrap id columns, then use unnest to take a look at our 100 decompositions. The columns 1 and 4 represent the one and four parity groups.

clogg_straps |> 
  mutate(results = map(clogg_agg, dg_clogg)) |> 
  select(id, results) |> 
  unnest(results) 
# A tibble: 300 × 6
   id           rowname             `1`   `4`  diff decomp
   <chr>        <chr>             <dbl> <dbl> <dbl>  <dbl>
 1 Bootstrap001 age_groups_struct  47.3  24.3 -23.0   39.3
 2 Bootstrap001 rate               53.6  18.1 -35.5   60.7
 3 Bootstrap001 crude              70.1  11.5 -58.5  100  
 4 Bootstrap002 age_groups_struct  51.2  26.8 -24.3   40.7
 5 Bootstrap002 rate               56.8  21.3 -35.5   59.3
 6 Bootstrap002 crude              71.9  12.1 -59.8  100  
 7 Bootstrap003 age_groups_struct  49.1  23.2 -25.9   41.0
 8 Bootstrap003 rate               54.7  17.5 -37.2   59.0
 9 Bootstrap003 crude              73.2  10.1 -63.1  100  
10 Bootstrap004 age_groups_struct  46.7  23.9 -22.8   38.2
# ℹ 290 more rows

We can then calculate the standard error of the standardization and decomposition for each of the factors in the decomposition.

clogg_boot_res <- 
clogg_straps |> 
 select(id, results) |> 
  unnest(results) |> 
 # filter(is.na(diff))
  group_by(rowname) |> 
  summarise(se_diff = sd(diff)) 

clogg_boot_res
# A tibble: 3 × 2
  rowname           se_diff
  <chr>               <dbl>
1 age_groups_struct    2.46
2 crude                1.94
3 rate                 3.28

And now we can attach these to the results from the original data, and calculate confidence intervals using a normal approximation.

original_res |> 
  rownames_to_column() |> 
  left_join(clogg_boot_res) |> 
  mutate(conf_low = diff - 1.96 * se_diff,
  conf_upp = diff + 1.96 * se_diff)
            rowname        1        4      diff decomp  se_diff  conf_low
1 age_groups_struct 48.61896 25.54727 -23.07169  38.07 2.462872 -27.89892
2              rate 55.84938 18.31685 -37.53253  61.93 3.278694 -43.95877
3             crude 72.09319 11.48897 -60.60422 100.00 1.944092 -64.41464
   conf_upp
1 -18.24446
2 -31.10629
3 -56.79380

We don’t have to use a normal approximation - we can use the bootstrap draws themselves to calculate confidence intervals (although be careful - you need a lot of bootstrap draws to get reliable estimates of 95% confidence intervals)

clogg_straps |> 
  select(id, results) |> 
  unnest(results) |> 
 # filter(is.na(diff))
  group_by(rowname) |> 
  reframe(
      estimate = quantile(diff, c(0.025, 0.5, 0.975), na.rm = TRUE),
      ci = c(0.025, 0.5, 0.975)
    ) |> 
  pivot_wider(id_cols = rowname,
  names_from = ci,
values_from = estimate)
# A tibble: 3 × 4
  rowname           `0.025` `0.5` `0.975`
  <chr>               <dbl> <dbl>   <dbl>
1 age_groups_struct   -27.0 -23.0   -18.1
2 crude               -64.2 -60.5   -57.3
3 rate                -44.5 -37.9   -31.6

Compare these to the values calculated with the Normal approximation above. How similar are they?

Stata

First read in the data

import delimited "https://raw.githubusercontent.com/benmatthewsed/sgsss_std_decomp/refs/heads/main/data/example2-clogg-long.csv", varnames(1) clear

Run the original decomposition from (Li 2017) Table 2.

rdecompose size rate, group(parity) transform(size) sum(age_groups)

Set up the decomposition program

Code from (Li 2017, p500) with variable names slightly amended.

program mydecompose, eclass
preserve
collapse (count) size= d (mean) rate = d, by(age_groups parity)
quietly rdecompose size rate, group(parity) transform(size) sum(age_groups)
matrix b = e(b) * 100
ereturn post b
restore
end

Explode the original data

This code expands out the dataset based on the size variable and then creates a new 0/1 column d based on the rate variable.

expand size

by age_groups parity, sort: generate d = _n<=round(rate*size/100)

We can check what we’ve done by summing up again by age and parity

by age_groups parity, sort: tab d

Bootstrapping the standardization and decomposition

Now we can resample the data and run the standardization and decomposition. Try changing the reps() value from 1000 to 100 and see what happens. The rseed command makes the results reproducible.

bootstrap, nowarn nodots reps(1000): mydecompose rseed 356

References

Li, Jinjing. 2017. “Rate Decomposition for Aggregate Data Using Das Gupta’s Method.” The Stata Journal 17 (2): 490–502.