I am working with censored data and would like to employ a selection model. As far as I can see, the most frequently applied selection model is the Heckman selection model that assumes a two stage process.
First we have a partly latent variable $Y\in \mathbb{R}$, where we only observe $\Delta Y \in \mathbb{R_+}$, with $\Delta=I(Y>0)$ beeing an indicator function that is one if $Y$ is greater than zero which is seen as equivalent to $Y$ beeing observed. Furthermore we assume that the latent variable $Y$ can be modelled using a linear model, using some covariates $X$: $Y=X\beta + \epsilon$ and that we can use covariates $Z$ in order to employ a probit model for the probability of response, i.e. $P(\Delta=1|Z)=\Phi(Z\gamma)$, with $\Phi(\cdot)$ denoting a standard normal CDF.
Using this, plus the assumption that the errors from both stages are jointly normal, as well as $Z\subset X$, we can arrive at $$E(\Delta Y|X,\Delta=1)=E(Y|X,\Delta=1)=X\beta + E(\epsilon|X,\Delta=1)=X\beta + \rho \sigma_\epsilon\lambda(Z\gamma) $$ with $\lambda(\cdot)$ representing the inverse Mill Ratio evaluated at $Z\gamma$. So far so good. But for some theoretical reasons I cannot assume in my application, that $P(\Delta=1|Z)$ is modelled correctly by a standard normal, but rather $P(\Delta=1|Z)=\Lambda(Z\gamma)$ with $\Lambda(\cdot)$ beeing the logistic function. Now the crucial assumption of joint normality of the errors breaks down as the errors from the logistic regression come from the type I generalized extreme value distribution.
Is there a theoretic possibility to fix this problem or do I have to rely on some numerical approximations for $E(\epsilon|X,\Delta=1)$
