# What does it mean that the decomposition is based on the linear systematic component? And how can I interpret my result?

I'm using the oaxaca package to implement a Blinder-Oaxaca decomposition on a logistic model with binary outcome.

The vignette says that:

Note that, if a non-linear function such as glm() is chosen, the decomposition will be based on the linear systematic component – usually associated with the estimation of the corresponding latent variable – of the regression method.

What does it mean that the decomposition is based on the linear systematic component?

I'm decomposing with:

results <- oaxaca(f,
data = df,
R = R,
reg.fun = glm,


The results is as follows:

> results$$n$$n.A
[1] 188751

$n.B [1] 87856$n.pooled
[1] 276607

> results$$y$$y.A
[1] 0.005260899

$y.B [1] 0.002663449$y.diff
[1] 0.00259745


So results$y looks like it contains probabilities. The difference in probability between the two group is 0.25%. But I don't understand the three-fold decomposition, which is as follows: > results$threefold$overall coef(endowments) se(endowments) coef(coefficients) se(coefficients) 0.59955787 0.12425443 0.23692063 0.13821627 coef(interaction) se(interaction) -0.04663127 0.13934151  If I sum the 3 components of the decomposition, I obtain: > te <- results$threefold$overall[["coef(endowments)"]] + + results$threefold$overall[["coef(coefficients)"]] + + results$threefold$overall[["coef(interaction)"]] > print(te) [1] 0.7898472  which should equal results$y$y.diff, but obviously isn't. ## 1 Answer To your first question - What does it mean that the decomposition is based on the linear systematic component? This means that the decomposition implemented here is on the log-odds (logit) scale. For a logistic regression, the $$\beta$$ coefficients are interpreted in terms of log-odds, which are linear and additive. The endowments, coefficients, and interaction outputs for your Blinder-Oaxaca decomposition are therefore all in log-odds terms. Addressing your second question about summing the decomposition estimates - this calculation is impacted by use of the log-odds scale, and the raw sum of these components will not be equal the simple difference in probabilities given by the results$y output.

To better understand what's happening, let's simulate a decomposition

# For reproducibility
set.seed(6000)

# Data generating process
n <- 10000
mu_group_A <- c(1, 2.5, 4)
mu_group_B <- c(0.5, 2.5, 3.5)
sigma <- rWishart(1,3,diag(3))[,,1]

# Simulate RHS observed data
x1 <- MASS::mvrnorm(n, mu_group_A,sigma)
x2 <- MASS::mvrnorm(n, mu_group_B,sigma)

# Regression formulas for two groups
intercept <- 1
betas_A <- c(0.25, 0.1, 0)
betas_B <- c(0.5, 0.75, 0.15)

# Calculate log-odds outcome for each
z1 <- intercept + x1 %*% betas_A
z2 <- intercept + x2 %*% betas_B

# Convert from log-odds to observed binary outcome
pr <- plogis(c(z1, z2))
y <- rbinom(n*2, 1, pr)

# Join into a data.frame
group_B <- rep(0:1, each = n)
df <- data.frame(y, rbind(x1, x2), group_B)


Now before we fit any models to our simulated data, we already know what the true decomposition components should be, based on the means of each variable by group ($$\bar{X}_A$$ and $$\bar{X}_B$$) and the log-odds scale $$\beta$$ coefficients we used to generate the outcome $$y$$.

\begin{align*} \text{Endowments} && (\bar{X}_A - \bar{X}_B)'\beta_B = 0.325 \\ \text{Coefficients} && \bar{X}_B' (\beta_A - \beta_B) = -2.275 \\ \text{Interaction} && (\bar{X}_A - \bar{X}_B)' (\beta_A - \beta_B) = -0.2 \\ \end{align*}

Which we can calculate in R by:

# True decomposition estimates
endowments <- t(mu_group_A - mu_group_B) %*% betas_B
coeffs <- t(mu_group_B) %*% (betas_A - betas_B)
interact <- t(mu_group_A - mu_group_B) %*% (betas_A - betas_B)

> (true_BO_par <- c(endowments, coeffs, interact))
[1]  0.325 -2.275 -0.200


Now let's run the Blinder-Oaxaca on our simulated data

library(oaxaca)
fit <- oaxaca(y ~ X1 + X2 + X3 | group_B,
data = df,
reg.fun = glm,
R = 100,


This takes a little while to run with the bootstrap procedure (used to estimate uncertainty on the decomposition components), but we eventually find

> fit$threefold$overall
coef(endowments)     se(endowments) coef(coefficients)
0.37101832         0.03174424        -2.48317386
se(coefficients)  coef(interaction)    se(interaction)
0.08839086        -0.23463968         0.02981978


Which agrees fairly well with our known decomposition parameters.

If we now take the sum of both the true and estimated decomposition components, you'll notice something interesting.

> sum(true_BO_par)
[1] -2.15
> sum(fit$threefold$overall[c(1, 3, 5)])
[1] -2.346795


Although they agree somewhat well with each other, they are definitely not on the probability scale (i.e. they don't equal fit$y$diff). But more importantly, transforming the summed decomposition components back to the probability scale doesn't actually give us the raw difference in probabilities.

This is because transforming back to the probability scale makes the previously independent linear effects of each term on the model (i.e., on the log-odds scale) dependent on the values of the other variables in the model. We need instead to use marginal effects.

To verify our results, let's estimate the separate marginal effects for each group, and find the difference

library(marginaleffects)

# Marginal effects needs a factor variable
df$$group_B <- as.factor(df$$group_B)

# Refit the models separately
m1 <- glm(y ~ X1 + X2 + X3,
data = df,
subset = group_B == 0,
m2 <- glm(y ~ X1 + X2 + X3,
data = df,
subset = group_B == 1,

> (mm1 <- marginalmeans(m1, type = "link"))
term value marginalmean  std.error conf.low conf.high p.value statistic
1 group_B     0     1.522338 0.02714284 1.469139  1.575537       0  56.08618
> (mm2 <- marginalmeans(m2, type = "link"))
term value marginalmean  std.error conf.low conf.high p.value statistic
1 group_B     1     3.869133 0.08093886 3.710496  4.027771       0  47.80316


And taking the difference in estimated marginal means for the separate models, we can see it matches our summed decomposition components (again, on the log-odds scale) perfectly:

> mm1$marginalmean - mm2$marginalmean
[1] -2.346795


To sum up, the decomposition components you're estimating by running a Blinder-Oaxaca model with a logistic regression have a marginal interpretation on the log-odds scale.