Introduction to rifexpectile

Motivation

Recentered influence function (RIF) regressions extend the Oaxaca-Blinder decomposition from the mean to other distributional statistics. Applied to quantiles, this requires estimating an inverse density term at each quantile, which can be noisy and irregular, especially in the tails. rifexpectile implements a density-free alternative based on unconditional expectiles: the expectile RIF depends only on primitive moments of the outcome and requires no density estimation, no bandwidth choice, and no kernel smoothing.

The underlying methodology (the closed-form influence function, the two-sample asymptotic theory, and an empirical application) is described in a companion manuscript; this vignette focuses on how to use the package, not on re-deriving that theory.

library(rifexpectile)

A first example with simulated data

The package ships with a small simulated data set, wage_gap_sim, with a synthetic outcome logy observed in two waves (“2000” and “2010”), along with covariates x1 (continuous) and x2 (binary), and a stylized sampling weight. No real or proprietary data are used anywhere in this package.

data(wage_gap_sim)
head(wage_gap_sim)
#>   group        y     logy         x1 x2    weight
#> 1  2000 5661.229 8.641396  0.7677843  1 1.1484742
#> 2  2000 3417.577 8.136687  1.0126557  0 0.9333160
#> 3  2000 2850.788 7.955351  1.1053997  0 1.1191789
#> 4  2000 2283.712 7.733557 -0.1492884  0 1.1719493
#> 5  2000 3539.620 8.171775  0.9296988  1 1.6370367
#> 6  2000 1907.758 7.553684 -2.0784772  0 0.7371935

The main function, rif_decompose(), takes a formula, a data frame, and the name of the column identifying the two groups to compare:

fit <- rif_decompose(
  logy ~ x1 + x2,
  data = wage_gap_sim,
  group = "group",
  ref_group = "2000",
  alpha = seq(0.1, 0.9, by = 0.2),
  boot_reps = 100,
  seed = 1
)
fit
#> Density-free RIF-Oaxaca-Blinder decomposition of unconditional expectiles
#> Groups: '2000' (reference) vs. '2010'  (n = 2000 vs. 2000)
#> Bootstrap replications: 100
#> 
#>  alpha  total composition structure
#>    0.1 0.2241      0.0042    0.2199
#>    0.3 0.2831      0.0043    0.2788
#>    0.5 0.3211      0.0043    0.3168
#>    0.7 0.3563      0.0043    0.3520
#>    0.9 0.4145      0.0044    0.4102

The total column is always exactly the difference of the two groups’ sample expectiles at that level (this is a finite-sample identity, not an approximation, whenever an intercept is included – see ?rif_decompose for the algebra). The composition and structure columns are the RIF-regression-based two-fold decomposition of that total into an endowment effect and a returns effect.

plot(fit, which = "structure")

fit$detail_structure
#>                      0.1         0.3         0.5         0.7         0.9
#> (Intercept)  0.218822969 0.274841531 0.311810926 0.346230726 0.395926268
#> x1           0.002500413 0.002945101 0.002913782 0.003024582 0.003303249
#> x2          -0.001450154 0.001010056 0.002026771 0.002750263 0.010938976
colSums(fit$detail_structure)  # matches fit$structure
#>       0.1       0.3       0.5       0.7       0.9 
#> 0.2198732 0.2787967 0.3167515 0.3520056 0.4101685

fit_w <- rif_decompose(
  logy ~ x1 + x2,
  data = wage_gap_sim,
  group = "group",
  ref_group = "2000",
  alpha = c(0.1, 0.5, 0.9),
  weights = "weight",
  boot_reps = 100,
  seed = 1
)
fit_w
#> Density-free RIF-Oaxaca-Blinder decomposition of unconditional expectiles
#> Groups: '2000' (reference) vs. '2010'  (n = 2000 vs. 2000)
#> Bootstrap replications: 100
#> 
#>  alpha  total composition structure
#>    0.1 0.2193      0.0011    0.2182
#>    0.5 0.3158      0.0011    0.3146
#>    0.9 0.4131      0.0011    0.4119

A second example with public data

To illustrate the package on a widely used, publicly available data set rather than only a simulated one, this section uses wage1 from the wooldridge package (Wooldridge, Introductory Econometrics), a standard cross-section of US wages with a female indicator – the classic setting for an Oaxaca-Blinder-type decomposition of the gender wage gap.

data(wage1, package = "wooldridge")
wage1$group <- ifelse(wage1$female == 1, "female", "male")

fit_gap <- rif_decompose(
  lwage ~ educ + exper + tenure,
  data = wage1,
  group = "group",
  ref_group = "male",
  alpha = c(0.1, 0.25, 0.5, 0.75, 0.9),
  boot_reps = 200,
  seed = 1
)
fit_gap
#> Density-free RIF-Oaxaca-Blinder decomposition of unconditional expectiles
#> Groups: 'male' (reference) vs. 'female'  (n = 274 vs. 252)
#> Bootstrap replications: 200
#> 
#>  alpha   total composition structure
#>   0.10 -0.2585     -0.0757   -0.1828
#>   0.25 -0.3339     -0.0928   -0.2411
#>   0.50 -0.3972     -0.1066   -0.2906
#>   0.75 -0.4384     -0.1142   -0.3242
#>   0.90 -0.4647     -0.1195   -0.3452

plot(fit_gap, which = "structure")

Here ref_group = "male" means the composition term uses male coefficients as the reference and asks how much of the (female - male) gap in log-wage expectiles is attributable to differences in observed characteristics (education, experience, tenure) versus differences in their estimated returns, at each point of the log-wage distribution rather than only at the mean.

Reproducibility of the underlying data

simulate_income_data() is the function used to generate wage_gap_sim; the generating script is in data-raw/wage_gap_sim.R in the package source, so the bundled example data set can be regenerated or modified directly.

str(simulate_income_data(n = 5, seed = 42))
#> 'data.frame':    10 obs. of  6 variables:
#>  $ group : Factor w/ 2 levels "A","B": 1 1 1 1 1 2 2 2 2 2
#>  $ y     : num  7085 2246 1403 1659 5706 ...
#>  $ logy  : num  8.87 7.72 7.25 7.41 8.65 ...
#>  $ x1    : num  1.371 -0.565 0.363 0.633 0.404 ...
#>  $ x2    : int  1 0 1 1 0 0 0 1 0 1
#>  $ weight: num  1.051 0.806 0.805 1.256 0.747 ...