---
title: "ssmodels: A package for fit the sample selection models"
output: 
  rmarkdown::html_vignette:
    fig_caption: yes
    force_captions: yes
    #highlight: pygments
    toc: yes
    #Sumário flutuante
    #toc_float: true
    #numerar seções
    number_sections: true
    #Mostrar ou esconder os códigos (show ou hide)
    #code_folding: hide
    #Diversos modelos de documentos ver outros em http://bootswatch.com/
    #theme: united
vignette: >
  %\VignetteIndexEntry{ssmodels: A package for fit the sample selection models}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
  \usepackage{amsmath}
bibliography: "../inst/REFERENCES.bib"
biblio-style: "apalike"
---

```{r setup, include=FALSE}
require(knitr)
require(kfigr)
library(kableExtra)
options(knitr.table.format = "latex")
knitr::opts_chunk$set(echo = TRUE, fig.align = "center", message=FALSE, 
                      warning=FALSE, fig.height=5, fig.width=7.2)
```

# Introduction


<p style="text-align: justify;">
In order to facilitate the adjustment of the sample selection models existing in the literature, we created the ssmodels package. Our package allows the adjustment of the classic Heckman model (@heckman1976common, @heckman1979sample), and the estimation of the parameters of this model via the maximum likelihood method and two-step method, in addition to the adjustment of the Heckman-t models, introduced in the literature by @marchenko2012heckman and the Heckman-Skew model introduced in the literature by @ogundimu2016sample. We also implemented functions to adjust the generalized version of the Heckman model, introduced by @bastosBarreto, that allows the inclusion of covariables to the dispersion and correlation parameters and a function to adjust the Heckman-BS model introduced by @bastos that uses the Birnbaum-Saunders distribution as a joint distribution of the selection and primary regression variables.
</p>

<!--

<p style="text-align: justify;">
Com o objetivo de facilitar o ajuste dos modelos de seleção amostral existentes na literatura, 
criamos o pacote ssmodels. Nosso pacote permite o ajuste do modelo de Heckman clássico, e a estimação dos parâmetros desse modelo via método de máxima verossimilhança e método de dois passos, além do ajuste dos modelos Heckman-t, introduzido na literatura em 2012 por @marchenko2012heckman e do modelo Heckman-Skew introduzido na literatura por @ogundimu2016sample em 2016. Implementamos, ainda, funções para o ajuste da versão generalizada do modelo de Heckman que permite a inclusão de covariáveis aos parâmetros de dispersão e correlação e uma função para o ajuste do modelo Heckman-BS que utiliza a distribuição Birnbaum-Saunders como distribuição conjunta das variáveis de seleção e regressão primária. 
</p>



<p style="text-align: justify;">
Comparamos nesse texto os resultados obtidos com o ajuste das funções do nosso pacote e funções do pacote sampleSelection de @sampleSelection a dados simulados e também a banco de dados reais que foram incluídos no pacote.
</p>

# Modelo de Heckman Generalizado

## Simulated data with fixed dispersion and correlation 

<p style="text-align: justify;">
Inicialmente simulamos variáveis explicativas com distribuição uniforme padrão para as equações de seleção e regressão. Fixamos valores de alguns parâmetros e utilizamos a função ${\it mvtnorm()}$ do pacote mvtnorm de @mvtnorm para simular variáveis $y_{1}$ e $y_{2}$ bivariadas com distribuição normal. A variável $y_{1}$ foi utilizada como resposta de interesse e $y_{2}$ como equação de seleção de tal forma que observamos $yo=y_{1}*1(y_{2}>0).$ Abaixo segue os resultados da simulação.
</p>

```{r}
library("mvtnorm") #Pacote para geração dos dados
set.seed(0) 
#Parâmetros utilizados:
sigma <- 2.9
rho   <- -0.5
n     <- 100
gamma0<- 2.2
gamma1<- -3
beta0 <- -1
beta1 <- 2.7
#Covariáveis
xs    <- runif(n)
xo    <- runif(n)
#Vetor de parâmetros
b     <- rbind(beta0,beta1)
g     <- rbind(gamma0,gamma1)
mu1   <- (as.numeric(model.matrix(~xo) %*% b))
mu2   <- (as.numeric(model.matrix(~xs) %*% g))
#Simulação das variáveis bivariadas com distribuição normal
eps   <- matrix(NA,n,2)
for(k in 1:n){
  eps[k,] <- mvtnorm::rmvnorm(1,
                            c(mu1[k], mu2[k]),
                            matrix(c(sigma, rho*sigma, rho*sigma, 1), 2, 2))
}
#eps <- rmvnorm(1000, c(0,0), matrix(c(sigma, rho*sigma, rho*sigma, 1), 2, 2))
y1 <- cbind(eps[,1])
y2 <- cbind(eps[,2])>0
#Variável de seleção
ys<-1*y2
#Variável de interesse primário
yo<-y1*ys
```

<p style="text-align: justify;">
Com os dados acima, obtemos um total de `r sum(1*yo==0)` observações censuradas. Ou seja, cerca de $30\%$ dos dados. É importante observar que as covariáveis são independentes e portanto a restrição de exclusão é garantida. Abaixo comparamos os resultados obtidos com o ajuste da função ${\it HeckmanCL()}$ do nosso pacote com a função ${\it selection()}$ do pacote sampleSelection. Observamos, como esperado, os mesmos resultados. Apresentamos na `r figr("fig1", TRUE, type="Figura")` o diagrama de dispersão dos dados, círculos negros são resultados observados e círculos vazios representam os dados não observados. A linha sólida representa o ajuste do modelo linear aos dados completos, ou seja, essa é a reta que melhor estima a variável de interesse. A reta tracejada, cor vermelha, representa o modelo de Heckman e a reta tracejada na cor azul representa o modelo linear simples ajustado aos dados observados. É possível inferir que na presença de dados censurados, o modelo de Heckman apresenta estimativas mais realistas do que o modelo linear simples ajustado somente aos dados observados.
</p>

```{r, message=FALSE}
library(ssmodels)
m0 <- HeckmanCL(ys~xs,yo~xo)
summary(m0)
library("sampleSelection")
m1 <- selection(ys~xs, yo ~xo, method = "ml")
summary(m1)
```


```{r fig1, fig.align= "center", fig.width = 7.2, fig.height= 5, anchor="Figura"}
opar <- par(mar=c(2,2,2,0) + 0.9, mgp=c(2,1,0))
pch <- c(1, 16)
plot(xo, y1, pch=pch[1 + ys], cex=0.5, lwd=0.5, main = "Figura 1: Ajuste dos modelos de Heckman Clássico e modelo linear \n simples a dados simulados com censura.")
# True dependence
abline(a=beta0, b=beta1, lty=1, lwd=2)
# Heckman's model
abline(a=coef(m0)[3], b=coef(m0)[4], lty=2, col="red", lwd=2)
# linear model
cf <- coef(lm(yo ~ xo, subset=ys==1))
abline(a=cf[1], b=cf[2], lty=3, col="blue", lwd=3)
par(opar)
```

## Data with variable dispersion and correlation
<p style="text-align: justify;">
Simulamos dados normais bivariados com dispersão e correlação variáveis. Nesse caso, comparamos o modelo de Heckman generalizado contra o modelo de Heckman clássico e o modelo linear simples. Observamos que o modelo de Heckman clássico apresentou estimativas dos parâmetros da equação principal semelhantes as estimativas do modelo linear simples, já nosso modelo, apresentou resultados próximos aos valores verdadeiros como podemos observar na `r figr("fig2", TRUE, type="Figura")`. 
</p>
```{r}
library("mvtnorm")
set.seed(0)
gamma0 <- 0.5
gamma1 <- -1
beta0 <- -2
beta1 <- 2
lambda0 <- 0.1
lambda1<- 0.5
kappa0 <- 0.3
kappa1 <- -0.5
a <- rbind(lambda0,lambda1)
b <- rbind(beta0,beta1)
g <- rbind(gamma0,gamma1)
d <- rbind(kappa0,kappa1)
n <- 200
xs <- rnorm(n)
xo <- rnorm(n)
mu1 <- (as.numeric(model.matrix(~xo) %*% b))
mu2 <- (as.numeric(model.matrix(~xs) %*% g))
sigma <- exp(as.numeric(model.matrix(~xo) %*% a))
rho   <- tanh(as.numeric(model.matrix(~xo) %*% d))
eps <- matrix(NA,n,2)
for(k in 1:n){
  eps[k,] <- mvtnorm::rmvnorm(1,c(mu1[k],mu2[k]),matrix(c((sigma[k])^2,rho[k]*sigma[k],rho[k]*sigma[k],1),2,2))
}
y1 <- cbind(eps[,1])
y2 <- cbind(eps[,2])>0
ys<-1*y2
yo<-y1*ys
```

No caso dos dados simulados acima, observamos `r sum(1*(yo==0))` valores censurados.

```{r, warning=FALSE}
m0 <- HeckmanCL(ys~1+xs, yo ~1+xo)
summary(m0)
m1 <- HeckmanGe(ys~xs, yo ~xo, cbind(xo), cbind(xo))
summary(m1)
```

Gráfico

```{r fig2, fig.align= "center", fig.width = 7.2, fig.height= 5, anchor="Figura"}
opar <- par(mar=c(2,2,3,0) + 0.9, mgp=c(2,1,0))
pch <- c(1, 16)
plot(xo, y1, pch=pch[1 + ys], cex=0.5, lwd=0.3, ylim = c(-8,8), main = "Figura 2: Ajuste dos modelos de Heckman Generalizados e Clássico \n e modelo linear simples a dados simulados com censura.")
# True dependence
abline(a=beta0, b=beta1, lty=1, lwd=2)
# Heckman's model
abline(a=coef(m0)[3], b=coef(m0)[4], lty=2, col="blue", lwd=2)
# linear model
m2 <- coef(lm(y1 ~ xo, subset=ys==1))
abline(a=m2[1], b=m2[2], lty=3, col="green", lwd=2)
# Heckman Generalized
abline(a=coef(m1)[3], b=coef(m1)[4], lty=5, col="red")
par(opar)
```

# Comparison of Parametric Sample Selection Models
<p style="text-align: justify;">
@heckman1976common introduziu o modelo de seleção amostral clássico. A maior crítica em relação ao modelo foi a sensibilidade da estimação dos parâmetros à suposição de normailidade dos erros, frequentemente violada na prática. Devido a isso, @heckman1979sample introduziu um método mais prático e eficiente para estimação de parâmetros, tal método ficou conhecido como método de dois passos e é utilizado principalmente como ferramenta para encontrar valores iniciais para algoritmos iterativos de estimação por máxima verossimilhança. Mesmo após os dois trabalhos de Heckman, houve ainda muita discussão em relação a distribuição bivariada dos termos de erro das equações de seleção e regressão. Diversos trabalhos sobre problemas de multicolinearidade, comparações do método de dois passos com outros procedimentos, problema de heterocedásticidade e outros foram publicados. Somente no trabalho de @marchenko2012heckman, 36 (trinta e seis) anos depois, conseguiram substituir a distribuição normal pela distribuição t de Student, que é uma distribuição mais robusta, no sentido de se ajustar melhor a distribuição de dados discrepantes (distantes da média). Após mais mais quatro anos, no trabalho de @ogundimu2016sample, foi introduzido o modelo Heckman-Skew, em que a distribuição normal foi substítuida pela distribuição Skew-Normal e o modelo foi sugerido para ser ajustado a dados assimétricos. 
</p>
<p style="text-align: justify;">
O maior problema dos dois modelos, skew e t de Student, é a necessidade da estimação dos parâmetros de assimetria e de graus de liberdade, respectivamente. Ambos, reconhecidamente difíceis de serem estimados. Além disso, mesmo o modelo Skew sendo para dados assimétricos, em seu ajuste, a distribuição da variável de interesse não pode ser altamente assimétrica, quando este é o caso, devemos transformar a variável de interesse na tentiva de torna-la mais simétrica, o que afeta a interpretação dos parâmetros. Devido a isso, em 2018, introduzimos o modelo Heckman-BS, em que utilizamos a distribuição Birnbaum-Saunders em substituição a distribuição normal. Nosso modelo apresenta várias vantagens quando comparado aos modelos existentes na literatura, além de ser parcimonioso, pois possui a mesma quantidade de parâmetros do modelo clássico de Heckman, ele também é para dados assimétricos e positivos. E a maioria dos bancos de dados em que o modelo de Heckman é ajustado satisfaz tais premissas. Além disso, não necessita da transformação da variável de interesse para seu ajuste, em nenhum caso, o que torna a interpretação dos parâmetros mais realista e direta. Vamos comparar daqui em diante o ajuste e a estimação dos parâmetros dos modelos clássico de Heckman e os modelos Heckman-Skew, Heckman-t e Heckman-BS. Vamos simular dados assimétricos e positivos e também ajustar tais modelos a banco de dados com problemas de viés de seleção existentes na literatura.
</p>
## Simulation of asymmetric and positive data
<p style="text-align: justify;">
Quando desejamos ajustar um modelo de regressão a uma variável resposta que não é observada em sua totalidade e desconfiamos que a não observação é devido a alguma outra variável não observada, podemos estar diante de um problema de viés de seleção. Neste caso, o ajuste do modelo de Heckman é a principal indicação. No entanto, tais dados, podem ainda ser assimétricos, ou possuir muitas observações distantes da média, ou ainda apresentar heterocedasticidade, o que torna suspeita qualquer estimação obtida com a modelagem clássica de Heckman. Veremos adiante a comparação dos resultados obtidos com o ajuste dos modelos Clássico, Heckman-Skew, Heckman-t e Heckman-BS a dados simulados positivos e assimétricos.  
</p>
```{r, warning=FALSE}
set.seed(0)
n=200 #Tamanho das amostras
#Valores iniciais dos parametros usados para gerar as variaveis. 
#Ou seja, Valor verdadeiro dos parametros
#Para manter 30% de censura, manter os seguintes parametros:
gamma0 <- 1.6
gamma1 <- 0.8
gamma2 <- 0.2
gamma3 <- 0.7
beta0  <- 1
beta1  <- 0.7
beta2  <- 1.1
phi1   <- 1.2
phi2   <- 1
rho    <- -0.5
lambda <- 1
nu <- 10
#Matriz de covariaveis para gerar mu1
X1 <- rep(1,n)
X2 <- rnorm(n,0,1)
X3 <- rnorm(n,0,1)
X4 <- rnorm(n,0,1)
XO <- cbind(X2,X3)
#Matriz de covariaveis para gerar mu2, sem restriao de exclusao
XS <- cbind(X2,X3,X4)
#Vetor de valores verdadeiros dos parametros
b <- rbind(beta0,beta1,beta2)
g <- rbind(gamma0,gamma1,gamma2,gamma3)
#Vetor de medias 1
mu1 <- exp(as.numeric(model.matrix(~XO) %*% b))
#Vetor de medias 2
mu2 <- exp(as.numeric(model.matrix(~XS) %*% g))

u1  <-rnorm(n)
u2  <-rnorm(n)
z1  <-(((sqrt(1+rho))+(sqrt(1-rho)))/2)*u1+(((sqrt(1+rho))-(sqrt(1-rho)))/2)*u2
z2  <-(((sqrt(1+rho))-(sqrt(1-rho)))/2)*u1+(((sqrt(1+rho))+(sqrt(1-rho)))/2)*u2
#################################################################
#Variaveis (T1,T2)~BSB(mu1,phi1,mu2,1,rho)
##########################################################
T1<-(mu1/(1+(1/phi1)))*((1/2)*(sqrt(2/phi1))*z1+sqrt(1+((1/2)*(sqrt(2/phi1))*z1)^2))^2
T2<-(mu2/(1+(1/phi2)))*((1/2)*(sqrt(2/phi2))*z2+sqrt(1+((1/2)*(sqrt(2/phi2))*z2)^2))^2
################################################################
#Variavel indicadora
######################################
YS<-1*(T2>1)
#######################################################
#Variavel de interesse cuja densidade eh obtida apartir da distribuicao
#Birnbaum Saunders bivariada, observe que esta eh uma variavel com censura
#e de acordo com a minha especificacao dos parametros a censura eh de
#aproximadamente 30%.
############################################################
YO<-T1*YS
```

Observamos `r sum(YO==0)` valores censurados, ou seja, cerca de `r (sum(YO==0)/length(YO))*100` porcento dos dados.

```{r, warning = FALSE}
#Data frame com os dados simulados
dt1=data.frame(YO,YS,XO,XS)
names(dt1) <- c("YO","YS","XO1","XO2","XS1","XS2","XS3")
selectionEq <- YS~XS1+XS2+XS3
outcomeEq <- YO~XO1+XO2
mBS <- HeckmanBS(selectionEq,outcomeEq, data=dt1)
#Transformacao de y para ajuste dos demais modelos 
l_YO <- ifelse(YS==1,log(YO),0)
#Variavel resposta para ajuste dos demais modelos
dt2=data.frame(l_YO,YS,XO,XS)
names(dt2) <- c("l_YO","YS","XO1","XO2","XS1","XS2","XS3")
selection <- YS~XS1+XS2+XS3
outcome <- l_YO~XO1+XO2
mCL <- HeckmanCL(selection, outcome, data=dt2)
mtS <- HeckmantS(selection, outcome, data=dt2, nu)
mSK <- HeckmanSK(selection, outcome, data=dt2, lambda)
```


Grafico

```{r fig3, fig.align= "center", fig.width = 7.2, fig.height= 5, anchor="Figura"}
opar <- par(mar=c(2,2,3,0) + 0.9, mgp=c(2,1,0))
library("ggplot2")
library("gridExtra")
barfill <- "grey"
barlines <- "black"
p1 <- ggplot(dt1, aes(YO)) + geom_histogram(aes(y = ..density..), bins=10, colour = barlines, fill = barfill)+
  scale_x_continuous(name = "(a) Values of the variable of interest")
  #                    breaks = seq(0, 125, 25),
  #                    limits=c(0, 125))+
  # scale_y_continuous(name = "Count",
  #                    breaks = seq(0, 0.01, 0.005),
  #                    limits=c(0,  0.01))
p2 <- ggplot(dt2, aes(l_YO)) + geom_histogram(aes(y = ..density..), bins=10, colour = barlines, fill = barfill) + scale_x_continuous(name = "(b) Log of values of variable of interest", limits=c(-5, 5))
  #                    breaks = seq(-5, 5, 1),
  #                    limits=c(-5, 5))+
  # scale_y_continuous(name = "Count",
  #                    breaks = seq(0, 1, 0.2),
  #                    limits=c(0, 1))
grid.arrange(p1, p2, ncol=2)
```
<p style="text-align: justify;">
O modelo Heckman-BS é ajustado aos dados simulados sem nenhuma transformação, os demais modelos são ajustados ao logaritmo dos dados, os resultados da estimação seguem abaixo. Observem que o modelo Heckman-BS apresenta os melhores resultados, além disso, o intercepto da equação de interesse primário é altamente viesado quando estimado pelos demais modelos. @marchenko2012heckman justifica a importância deste parâmetro na interpretação prática do ajuste do modelo, principalmente em aplicações Econômicas. 
</p>
```{r}
Parameters <- c("$\\gamma_{0}$", "$\\gamma_{1}$", "$\\gamma_{2}$", "$\\gamma_{3}$", "$\\beta_{0}$", "$\\beta_{1}$", "$\\beta_{2}$", "$\\phi$", "$\\rho$", "$\\lambda$", "$\\nu$")
truevalue <- c(gamma0, gamma1, gamma2, gamma3, beta0, beta1, beta2, phi1, rho, lambda, nu)
HBS <- round(mBS$coefficients, digits = 3)
HCL <- round(mCL$coefficients, digits = 3)
HSK <- round(mSK$coefficients, digits = 3)
HtS <- round(mtS$coefficients, digits = 3)
Results <- data.frame("Parameters"= Parameters,
                      "truevalue" = truevalue,
                      "HeckmanBS" = c(HBS, "NA", "NA"), 
                      "HeckmanCL" = c(HCL, "NA", "NA"),
                      "HeckmantS" = c(HtS[1:9], "NA", HtS[10]),
                      "HeckmanSK" = c(HSK[1:9], HSK[10], "NA"))
kable(Results, format = "html", align = c("c", "c", "c", "c", "c"))
#kable_styling(, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = TRUE)
```
<p style="text-align: justify;">
É importante ainda observar que os valores estimados do parâmetro de graus de liberdade do modelo Heckman-t e do parâmetro estimado pelo modelo Heckman-Skew indicam que os dados transformados possuem distribuição normal e, provavelmente por isso, apresentam resultados tão próximos do modelo clássico de Heckman. 
</p>
-->

# Real Databases

## MEPS 2001
<p style="text-align: justify;">
The MEPS is a set of large-scale surveys of families, individuals and their medical providers (doctors, hospitals, pharmacies, etc.) in the United States. It has data on the health services Americans use, how often they use them, the cost of these services and how they are paid, as well as data on the cost and reach of health insurance available to American workers. The sample is restricted to persons aged between 21 and 64 years and contains a variable response with 3328 observations of outpatient costs, of which 526 (15.8%) correspond to unobserved expenditure values and identified as zero expenditure for adjustment of the models. It also includes the following explanatory variables:
</p>
* educ: education status
* age: Age
* income: income
* female: gender
* vgood: a numeric vector
* good: a numeric vector
* hospexp: a numeric vector
* totchr: number of chronic diseases
* ffs: a numeric vector
* dhospexp: a numeric vector
* age2: a numeric vector
* agefem: a numeric vector
* fairpoor: a numeric vector
* year01: a numeric vector
* instype: a numeric vector
* ambexp: a numeric vector
* lambexp: log ambulatory expenditures
* blhisp: ethnicity
* instype_s1: a numeric vector
* dambexp: dummy variable, ambulatory expenditures
* lnambx: a numeric vector
* ins: insurance status
<p style="text-align: justify;">
These data were also used by @colin2009microeconometrics, @marchenko2012heckman and @zhelonkin2016robust to fit the classical Heckman and Heckman-t models and the robust version of the two-step method, respectively. We use the variable of interest $Y_{1}^{*}={\rm ambexp}$, which represents expenditure on medical services, on the logarithmic scale, since it is highly asymmetric, see `r figr("fig4", TRUE, type="Figura")`. The variable $Y_{2}^{*}={\rm dambexp},$ representing the willingness to spend, is not observed. We observe $U=1\{Y_{2}^{*}>0\},$ which represents the decision to spend on medical care.
</p>
```{r fig4, fig.align= "center", fig.width = 7.2, fig.height= 5, anchor="Figura"}
library(ssmodels)
#Leitura do dados MEPS2001
data(MEPS2001)
#tornando visiveis as colunas do data-frame
attach(MEPS2001)
barfill <- "grey"
barlines <- "black"
p1 <- ggplot(MEPS2001,aes(ambexp))+geom_histogram(colour = barlines, fill = barfill)+
  scale_x_continuous(name = "(a) Expenditures Medical",
                     breaks = seq(0, 15000, 2500),
                     limits=c(0, 15000))+
  scale_y_continuous(name = "Count",
                     breaks = seq(0, 800, 100),
                     limits=c(0, 800))
p2 <- ggplot(MEPS2001,aes(lambexp))+geom_histogram(colour = barlines, fill = barfill)+
  scale_x_continuous(name = "(b) Log of Expenditures Medical",
                     breaks = seq(0, 11, 1),
                     limits=c(0, 11))+
  scale_y_continuous(name = "Count",
                     breaks = seq(0, 300, 100),
                     limits=c(0, 300))
grid.arrange(p1, p2, ncol=2)
```
<p style="text-align: justify;">
According @marchenko2012heckman and @zhelonkin2016robust it is natural to fit a sample selection model to such data, since the willingness to spend $(Y_{2}^{*})$ is likely to be related to the expense amount $(Y_{1}^{*}).$ However, after fitting the classic Heckman model and using the Wald, likelihood ratio, or gradient tests for 
$H_{0}:\rho=0$ against $H_{1}:\rho\neq 0,$ the conclusion is that there is not sufficient evidence $(p\textrm{-valor}>0.1)$ to reject $H_{0},$ that is, there is no bias of selection. @colin2009microeconometrics suspected this conclusion and @marchenko2012heckman argued that a more robust model could evidence the 
selection bias present in the data and rejected the normality hypothesis of the data. However, we understand that the adjustment of the dispersion and/or the correlation can be an alternative to study these data without the need to modify the assumption of normality of the errors, because as our simulations indicate, the estimations of the parameters obtained by the adjustment of the model Heckman can be more severely affected by heteroscedasticity than by incorrect distribution of error terms.
</p>
<p style="text-align: justify;">
However, due to data asymmetry, we also adjusted the Heckman-BS model. Thus, we consider the following system of equations for classic and generalized Heckman models, Heckman-t and Heckman-Skew:
</p>
\begin{align}
lnambx &= \beta_{1}+ \beta_{2}age + \beta_{3}female + \beta_{4}educ + \beta_{5}blhisp + \beta_{6}totchr + \beta_{7}ins + \epsilon_{1i}, \label{eq_reg_aplic1} \\
dambexp&= \gamma_{1}+\gamma_{2}age + \gamma_{3}female + \gamma_{4}educ + \gamma_{5}blhisp + \gamma_{6}totchr + \gamma_{7}ins + \gamma_{8}income + \epsilon_{2i}, \label{eq_sel_aplic1} 
\end{align}
where errors are defined for classic Heckman models, Heckman-t and Heckman-Skew according to 
\begin{align}\label{2:disterro}
\begin{pmatrix}\epsilon_{1i}\\
\epsilon_{2i}
\end{pmatrix} &\overset{ind.}{\sim} \mathcal{N}
\begin{bmatrix}
\begin{pmatrix}
0\\
0
\end{pmatrix}\!\!,&
\begin{pmatrix}
\sigma^{2} & \rho\sigma \\
\sigma\rho & 1
\end{pmatrix}
\end{bmatrix}, i=1,\cdots,n,
\end{align}
and for generalized Heckman models how
\begin{align}\label{2:disterro1}
\begin{pmatrix}\epsilon_{1i}\\
\epsilon_{2i}
\end{pmatrix} &\overset{ind.}{\sim} \mathcal{N}
\begin{bmatrix}
\begin{pmatrix}
0\\
0
\end{pmatrix}\!\!,&
\begin{pmatrix}
\sigma_{i}^{2} & \rho_{i}\sigma_{i} \\
\sigma_{i}\rho_{i} & 1
\end{pmatrix}
\end{bmatrix}, i=1,\cdots,n,
\end{align}
in such a way that
\begin{align*}
\log{(\sigma_{i})} &= \lambda_{1}+\lambda_{2}age + \lambda_{3}female + \lambda_{4}totchr + \lambda_{5}ins\quad \textrm{and}\\ 
\rho_{i}&=\tanh{(\kappa_{0})}, i=1,\cdots, 3328.
\end{align*}

To adjust the Heckman-BS model we consider
\begin{align*}
\log\mu_{1i} &= \beta_{1} + \beta_{2}age + \beta_{3}female  + \beta_{4}educ + \beta_{5}blhisp + \beta_{6}totchr + \beta_{7}ins,\\
\log\mu_{2i} &= \gamma_{1}+\gamma_{2}age + \gamma_{3}female + \gamma_{4}educ + \gamma_{5}blhisp + 
\gamma_{6}totchr + \gamma_{7}ins + \gamma_{8}income,\quad i=1,2,\cdots,3328. 
\end{align*}
<p style="text-align: justify;">
and $(ambexp,dambexp)$ is independently distributed with Birnbaum-Saunders (BS) Bivariate parameter distribution $\mu_{1}>0, \mu_{2}>0, \phi_{1}=\phi>0, \phi_{2}=1$ and $-1<\rho<1,$ that is, 
$$(ambexp_{i},dambexp_{i})\stackrel{ind}{\sim}\mbox{ BS}(\mu_{1i},\mu_{2i}, 
\phi,1,\rho), \mu_{1i}>0, \mu_{2i}>0, \phi>0\quad \textrm{and}\quad -1<\rho<1, i=1,\cdots, n$$
</p>
<p style="text-align: justify;">
We observed similar values for the parameters estimated by all models, mainly the parameters of the variable of interest. However, the only model that allows a direct interpretation of such parameters is the Heckman-BS model. We can, for example, affirm, based on the result of this model, that keeping the other parameters fixed, changing a unit in age represents an increase of $24.6\%$ in ambulatory expenses. For the other models, the interpretation is related to the log of outpatient expenses. See the results of the estimates of the parameters below:
</p>
```{r, warning=FALSE}
selectEq <- dambexp ~ age + female + educ + blhisp + totchr + ins + income
outcomeEq <- lnambx ~ age + female + educ + blhisp + totchr + ins 
outcomeS <- cbind(age,female,totchr,ins)
outcomeC <- 1
outcomeBS <- ambexp ~ age + female + educ + blhisp + totchr + ins 
mCL <- HeckmanCL(selectEq, outcomeEq, data = MEPS2001)
mBS <- HeckmanBS(selectEq, outcomeBS, data = MEPS2001)
mSK <- HeckmanSK(selectEq, outcomeEq, data = MEPS2001,lambda = 1)
mtS <- HeckmantS(selectEq, outcomeEq, data = MEPS2001,df=12)
mGe <- HeckmanGe(selectEq, outcomeEq,outcomeS, outcomeC, data = MEPS2001)
Parameters <- c("Intercept", "age", "female", "educ", "blhisp", "totchr", "ins", "income",
   "Intercept", "age", "female", "educ", "blhisp", "totchr", "ins", "sigma", "age", "female", "totchr", "ins", "rho", "nu", "lambda")
HBS <- round(mBS$coefficients, digits = 3)
HCL <- round(mCL$coefficients, digits = 3)
HSK <- round(mSK$coefficients, digits = 3)
HtS <- round(mtS$coefficients, digits = 3)
HGe <- round(mGe$coefficients, digits = 3)
Results <- data.frame("Parameters"= Parameters,
                      "HeckmanGe" = c(HGe[1:21], "NA", "NA"), 
                      "HeckmanCL" = c(HCL[1:16], "NA", "NA", "NA", "NA", HCL[17], "NA", "NA"),
                      "HeckmanBS" = c(HBS[1:16], "NA", "NA", "NA", "NA", HBS[17], "NA", "NA"), 
                      "HeckmantS" = c(HtS[1:16], "NA", "NA", "NA", "NA", HtS[17:18], "NA" ),
                      "HeckmanSK" = c(HSK[1:16], "NA", "NA", "NA", "NA", HSK[17], "NA", HSK[18]))
kable(Results, format = "html", align = c("c", "c", "c", "c", "c"))
```
<p style="text-align: justify;">
In addition, there are also the results of the significance of the parameters. The classic Heckman model adjusted to such data indicates no correlation between the variables value spent and the decision to spend. According to @colin2009microeconometrics, such a result is suspect and should be better analyzed. 
</p>
```{r, warning=FALSE}
summary(mCL)
```
<p style="text-align: justify;">
The generalized Heckman model, on the other hand, indicates that there is the presence of heteroscedasticity in the data and, probably, due to this, the classic Heckman model is not the most suitable for modeling these data. Our model, in addition to indicating the presence of variable dispersion, also indicates a significant correlation between the variables amount spent and the decision to spend.
</p>
```{r, warning=FALSE}
summary(mGe)
```
<p style="text-align: justify;">
On the other hand, as the observed data show asymmetry, the Heckman-BS model can be a good indication to adjust such data. Note that the Heckman-BS model also indicates the presence of a significant correlation. It also presents, as previously mentioned, the advantage of parsimony and the direct interpretation of the relationship of the parameters with the response variable of interest.
</p>
```{r, warning=FALSE}
summary(mBS)
```
<p style="text-align: justify;">
The Heckman-t model, even adjusted to the expense log, indicates the presence of deviation from normality, since the degree of freedom parameter is approximately equal to $13$ and is significant. Such a model also observes the significant correlation parameter.
</p>
```{r, warning=FALSE}
summary(mtS)
```
<p style="text-align: justify;">
Finally, the Heckman-Skew model also indicates the presence of a significant correlation between the variables of interest and selection. In addition, it also indicates deviation from the normality of the transformed data. It is important to note that without the data transformation, the iterative method of estimating the parameters of the Heckman-Skew model does not converge.
</p>
```{r, warning=FALSE}
summary(mSK)
```


## Nhanes 2003-2004
<p style="text-align: justify;">
The US National Health and Nutrition Examination Study (NHANES) is a survey data collected by the US National Center for Health Statistics. The survey data dates back to 1999, where individuals of all ages are interviewed in their home annually and complete the health examination component of the survey. The study variables include demographic variables (e.g. age and annual household income), physical measurements (e.g. BMI – body mass index), health variables (e.g. diabetes status), and lifestyle variables (e.g. smoking status). This data frame contains the following columns:
</p>
* id: Individual identifier
* age: Age
* gender: Sex 1=male, 0=female
* educ: Education is dichotomized into high school and above versus less than high school
* race: categorical variable with five levels
* income: Household income ($1000 per year) was reported as a range of values in dollar (e.g. 0–4999, 5000–9999, etc.)
* and had 10 interval categories. 
*Income: Household income ($1000 per year) was reported as a range of values in dollar (e.g. 0–4999, 5000–9999, etc.)
* and had 10 interval categories. 
* bmi: body mass index
* sbp: systolic blood pressure

```{r, warning=FALSE}
library(ssmodels)
data(nhanes)
attach(nhanes)
perc <- function(x,data){
    nna <- ifelse(sum(is.na(x))!=0,summary(x)[[7]],0)
    perc <- ifelse(sum(is.na(x))!=0,(nna/length(data$id))*100,0)
   return(perc)
}
Variables <- c("SBP (mm Hg)", "Age (year)", "Gender", "BMI (Kg/$m^{2}$)", "Education (years)", "Race", "Income ($\\$1000$ per year)", "Numbers Obs.")
perc1 <- round(perc(sbp,nhanes), digits = 2)
perc2 <- round(perc(age,nhanes), digits = 2)
perc3 <- round(perc(gender,nhanes), digits = 2)
perc4 <- round(perc(bmi,nhanes), digits = 2)
perc5 <- round(perc(educ,nhanes), digits = 2)
perc6 <- round(perc(race,nhanes), digits = 2)
perc7 <- round(perc(Income,nhanes), digits = 2)
nObs <- length(Income)
Percentage <- c(perc1, perc2, perc3, perc4, perc5, perc6, perc7, nObs)
df <- subset(nhanes, !is.na(sbp))
df <- subset(df, !is.na(bmi))
attach(df)
perc11 <- round(perc(sbp,df), digits = 2)
perc12 <- round(perc(age,df), digits = 2)
perc13 <- round(perc(gender,df), digits = 2)
perc14 <- round(perc(bmi,df), digits = 2)
perc15 <- round(perc(educ,df), digits = 2)
perc16 <- round(perc(race,df), digits = 2)
perc17 <- round(perc(Income,df), digits = 2)
nObs1 <- length(Income)
Percentage1 <- c(perc11, perc12, perc13, perc14, perc15, perc16, perc17, nObs1)
table <- data.frame("Variables" = Variables, "Percentage of Missing"= Percentage, "Without Missing"= Percentage1)
kable(table, format = "html", align = c("c", "c", "c"))
```

```{r, warning=FALSE}
library("ggplot2")
library("gridExtra")
barfill <- "grey"
barlines <- "black"
p1 <- ggplot(df, aes(Income)) + geom_histogram( breaks = seq(0, 10, 0.5), aes(y = ..density..), colour = barlines, fill = barfill)+
  scale_x_continuous(name = "Income",
                   breaks = seq(0, 10, 2),
                   limits=c(0, 10))

p2 <- ggplot(df, aes(log(sbp))) + geom_histogram( colour = barlines, fill = barfill) +
  scale_x_continuous(name = "Log Systolic blood pressure")
grid.arrange(p1, p2, ncol=2)
```

```{r, warning = FALSE}
df$YS <- ifelse(is.na(df$Income),0,1)
df$educ <- ifelse(df$educ<=2,0,1)
df$Income <- ifelse(is.na(df$Income),0,df$Income)
attach(df)
selectionEq <- YS~age+gender+educ+race
outcomeEq   <- log(sbp)~age+gender+educ+bmi+Income
outcomeBS   <- sbp~age+gender+educ+bmi+Income
mCL <- HeckmanCL(selectionEq, outcomeEq, data = df)
mBS <- HeckmanBS(selectionEq, outcomeBS, data = df)
mSK <- HeckmanSK(selectionEq, outcomeEq, data = df, lambda = 0)
mtS <- HeckmantS(selectionEq, outcomeEq, data = df, df = 15)

    Parameters <- c("Intercept", "age", "gender", "educ", "race", "Intercept", "age", "gender", "educ", "bmi", "income", "sigma", "rho", "nu", "lambda")
HBS <- round(mBS$coefficients, digits = 5)
HCL <- round(mCL$coefficients, digits = 5)
HSK <- round(mSK$coefficients, digits = 5)
HtS <- round(mtS$coefficients, digits = 5)
Results <- data.frame("Parameters"= Parameters, 
                      "HeckmanCL" = c(HCL[1:13], "NA", "NA"),
                      "HeckmanBS" = c(HBS[1:13], "NA", "NA"), 
                      "HeckmantS" = c(HtS[1:13], HtS[14], "NA"),
                      "HeckmanSK" = c(HSK[1:13], "NA", HSK[14]))
kable(Results, format = "html", align = c("c", "c", "c", "c", "c"))
```

```{r, warning = FALSE}
summary(mCL)
```

```{r, warning = FALSE}
summary(mtS)
```

```{r, warning = FALSE}
summary(mBS)
```

```{r, warning = FALSE}
summary(mSK)
```



# References
