# Heckman with second step probit in R [closed]

The functions selection and heckit (package sampleSelection) support a binary dependent variable in the outcome equation:

The dependent variable of of the selection equation (specified by argument selection) must have exactly two levels (e.g., 'FALSE' and 'TRUE', or '0' and '1'). By default the levels are sorted in increasing order ('FALSE' is before 'TRUE', and '0' is before '1'). If the dependent variable of the outcome equation (specified by argument outcome) has exactly two levels, this variable is modelled as a binary variable.

Nevertheless, I can't find the possibility to choose the binary model to use in the second step. In particular I would like to use a probit model instead of a linear probability model.
Example using binary variable in the outcome equation:

## We estimate the probability of womens' education on their chances to get high wage (> \$5/hr in 1975 USD), using PSID data.
## We use education as explanatory variable and add age, kids, and non-work income as exclusion restrictions.

library(stargazer)
library(sampleSelection)

data(Mroz87)

ols=lm(wage >= 5 ~ educ, data=subset(Mroz87, lfp==1) )

m <- heckit(lfp ~ educ + age + kids5 + kids618 + nwifeinc,
wage >= 5 ~ educ, data = Mroz87 )

stargazer(ols, m, title="Married women's wage regressions", type="text",df=FALSE, digits=4)


As you can see by the result of the code in the example, the function uses the lpm and I can't find a way to modify it.

I do not know any package that does this. However, you can manually create a two-step probit. In Heckman's model, you take the inverse mill's ratio(imr) and pass it to the second step, which is a linear model. You can manually create the first step as a probit (confirm your findings and make sure you are on the right track), derive the imr, and write the second step. I will give it a try:

selection_fml <- lfp ~ educ + age + kids5 + kids618 + nwifeinc
probit_selection <- glm(selection_fml, data=Mroz87, family= binomial(link='probit'))

probit_lp <- -predict(probit_selection )
inv_mills <- dnorm(probit_lp)/(1-pnorm(probit_lp))

probit_outcome <- glm(wage >= 5 ~ educ + inv_mills, data = Mroz87, family= binomial(link='probit'))
summary(probit_outcome)


This should work.

I used a two-step probit before. I think it makes sense, especially if the second step is not an OLS-type situation. However, comparison can get complicated since it is hard to compare it with a single-step probit or Heckman's two-step model. AIC seems the only good way to do this. Using packages that help you print a decent table like stargazer is another difficulty if you do this manually.