Selection Models — MNAR (R)

sampleSelection — Heckman's two-step and maximum likelihood

The R counterpart to selection_python.ipynb. Where the Python notebook builds the Bayesian selection sampler from scratch, this notebook uses sampleSelection (Toomet & Henningsen) — the package that implements Heckman's classic estimators: the two-step procedure (probit selection, then an outcome regression augmented with the inverse-Mills ratio) and the full-information maximum likelihood fit. Same benchmark: the Mroz wage data, where the log wage is observed only for working women.

options(repr.plot.width=9, repr.plot.height=4.5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(sampleSelection))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; GREY<-"#718096"
data(Mroz87); Mroz87$lwage <- ifelse(Mroz87$lfp==1, log(Mroz87$wage), NA)
cat("Mroz87:", nrow(Mroz87), "women;", sum(Mroz87$lfp), "working (log wage observed), ",
    sum(Mroz87$lfp==0), "not (missing).\n")

Warning message:
"package 'sampleSelection' was built under R version 4.6.1"

Warning message:
"package 'maxLik' was built under R version 4.6.1"

Warning message:
"package 'miscTools' was built under R version 4.6.1"

Mroz87: 753 women; 428 working (log wage observed),  325 not (missing).

1. Heckman selection vs naive OLS

We estimate the wage equation $\log \text{wage}\sim\text{educ}+\text{exper}+\text{exper}^2$ with participation (the selection) depending additionally on age, young children and non-wife income — the exclusion restrictions. selection() fits it by maximum likelihood; lm() on the working women alone is the naive comparison. The estimated rho and the significance of the inverse-Mills term say whether selection matters.

ml <- selection(lfp ~ educ + exper + I(exper^2) + age + kids5 + kids618 + nwifeinc,
                lwage ~ educ + exper + I(exper^2), data=Mroz87, method="ml")
ts <- selection(lfp ~ educ + exper + I(exper^2) + age + kids5 + kids618 + nwifeinc,
                lwage ~ educ + exper + I(exper^2), data=Mroz87, method="2step")
ols <- lm(lwage ~ educ + exper + I(exper^2), data=Mroz87)
cat("returns to education (educ coefficient in the log-wage equation):\n")
cat(sprintf("  Heckman ML   : %.3f\n", coef(ml)["educ"]))
cat(sprintf("  Heckman 2step: %.3f\n", coef(ts)["educ"]))
cat(sprintf("  naive OLS    : %.3f\n", coef(ols)["educ"]))
cat(sprintf("\nestimated selection correlation rho = %.2f (ML)\n", ml$estimate["rho"]))
print(summary(ml)$estimate[c("rho","sigma"),])

returns to education (educ coefficient in the log-wage equation):

  Heckman ML   : 0.131

  Heckman 2step: 0.131

  naive OLS    : 0.107

estimated selection correlation rho = 0.03 (ML)

        Estimate Std. Error    t value      Pr(>|t|)
rho   0.02660697  0.1470779  0.1809039  8.564926e-01
sigma 0.66339757  0.0227075 29.2149123 2.494880e-125

The rho estimate and its standard error are the crux — and they expose a subtlety. Here rho is not distinguishable from zero, so there is little evidence of selection on unobservables; yet the Heckman returns to education still sit somewhat above the OLS value. The inverse-Mills-ratio term is a nonlinear function of the selection covariates — which include education — so it is collinear with education in the wage equation, and adding it shifts the coefficient even when the selection itself is weakly identified. This is the well-known fragility of Heckman corrections under weak identification: the point estimate can move without the data really supporting a selection effect.

est<-c("Heckman ML"=coef(ml)["educ"], "Heckman 2step"=coef(ts)["educ"], "naive OLS"=coef(ols)["educ"])
barplot(est, col=c(GREEN,BLUE,RED), ylab="returns to education", main="Mroz: Heckman shifts the estimate even though rho is insignificant")
rho_se <- summary(ml)$estimate["rho","Std. Error"]
cat(sprintf("rho = %.2f (SE %.2f) -> %s\n", ml$estimate["rho"], rho_se,
    ifelse(abs(ml$estimate["rho"]/rho_se) < 2, "NOT significant: little evidence of selection on unobservables", "significant")))
cat(sprintf("returns to education: OLS %.3f  ->  Heckman %.3f  (a %.0f%% shift from an insignificant rho)\n",
    coef(ols)["educ"], coef(ml)["educ"], 100*(coef(ml)["educ"]/coef(ols)["educ"]-1)))
cat("So the maximum-likelihood correction moves the estimate even though selection is weakly identified -- the\n")
cat("inverse-Mills term is collinear with education. The Bayesian sampler in the Python notebook, regularised by its\n")
cat("prior, instead stays near the OLS value (~0.108): same weak-identification problem, two different symptoms.\n")

rho = 0.03 (SE 0.15) -> NOT significant: little evidence of selection on unobservables

returns to education: OLS 0.107  ->  Heckman 0.131  (a 22% shift from an insignificant rho)

So the maximum-likelihood correction moves the estimate even though selection is weakly identified -- the

inverse-Mills term is collinear with education. The Bayesian sampler in the Python notebook, regularised by its

prior, instead stays near the OLS value (~0.108): same weak-identification problem, two different symptoms.

2. Summary

sampleSelection fits the field-standard Heckman estimators and, on Mroz, tells a subtler story than "the correction does nothing". The estimated selection correlation rho is not distinguishable from zero — the data cannot identify selection on unobservables — yet the maximum-likelihood and two-step returns to education (~0.13) sit noticeably above OLS (~0.11), because the inverse-Mills-ratio term is collinear with education. That gap, from an insignificant rho, is the classic fragility of Heckman corrections under weak identification. The Bayesian sampler in selection_python.ipynb faces the same weak identification but, regularised by its prior, keeps the estimate near OLS — the wide posterior for rho being its honest signal that the data are uninformative. Two estimators, one weakly-identified problem, two different symptoms.

sampleSelection (with the selection() function and its heckit alias) is the standard tool for Heckman models and the frequentist mirror of the Bayesian selection sampler. Selection models factor the missing-data problem as $p(y)\,p(r\mid y)$ and generalise the Tobit model (a selection model with $\rho=1$). The correction rests on an exclusion restriction and the joint-normality assumption — untestable choices that make MNAR a matter of judgement. The next project, pattern-mixture models, makes that judgement an explicit sensitivity parameter.