Multi-stage selection model with panel data in R My data consists of individual level observations nested within countries over time. I would like to use multilevel models along with some sort of selection model. 
I have three related questions. 
1) Are there any issues or concerns with using a Heckman selection model twice? 
I have a model with two selection stages. While I can't conceive of any reason why this shouldn't work, I also haven't seen any examples of a heckman selection model being used twice. Assume I also have the necessary exclusion restrictions. 
2) My primary tool is R. There exists a package, sampleSelection, but my understanding is that algorithm it uses may not be appropriate for panel data / multilevel models. 
In stata, my understanding is the aforementioned is possible, but I'm not sure a package in R exists for doing heckman selection with multilevel models. 
3) Are there better ways for accomplishing what I am trying to do? Perhaps using instruments alone would be better? If so, are there any R packages you could point me to? If there are no R packages, I am still interested in this answer since I would be happy to code something if necessary. 
Thank you. 
 A: Looking at the Heckman article I see no reason, why two selection rules cannot be applied. In the article Heckman postulates the model as:
$Y_{1i}=X_{1i}\beta_1+U_{1i}$
$Y_{2i}=X_{2i}\beta_1+U_{2i}$
The goal is to estimate the first equation, but some data is missing. Which data is missing is controled by second equation. Then we have
$E(Y_{1i}| X_{1i}, \text{sample selection rule})=X_{1i}\beta_1+E(U_{1i}|\text{sample selection rule})$
If sample selection rule is $Y_{2i}\ge 0$ we have
$E(Y_{1i}|X_{1i},Y_{2i}\ge 0)=X_{1i}\beta_1+E(U_{1i}|U_{2i}\ge -X_{2i}\beta_2)$
And the article goes into details how this is estimated. From this I see no reason why we cannot add another sample selection rule $Y_{3i}\ge 0$, where:
$Y_{3i}=X_{3i}\beta_3+U_{3i}$
Then we get
$E(Y_{1i}|X_{1i},Y_{2i}\ge 0,Y_{3i}\ge 0)=X_{1i}\beta_1+E(U_{1i}|U_{2i}\ge -X_{2i}\beta_2,U_{3i}\ge X_{3i}\beta_3)$
Everything depends what assumptions are laid on $U_{1i},U_{2i},U_{3i}$. For estimation we can form maximum likelihood, or look for analog of formulas with Mills ration for bivariate normal density for trivariate normal density. In the latter case proceed as in article.  I would be quite optimistic of this idea succeeding. 
For the panel data framework look at Wooldridge book "Econometric analysis of cross-section and panel data". It has whole chapter dedicated for sample selection problems for panel data. The way I see it, it is straightforward application of ideas in Heckman's article.
As for the code in R, tough luck. If it is not implemented in plm package, it probably does not have readily available implementation.
