Modelo de Seleção de Heckman

Fernando de Souza Bastos -Universidade Federal de Minas Gerais

21 de Outubro, 2016

Variáveis Dependentes Limitadas

Como sabemos, variável dependente é uma medida que resulta do valor de outra medida variável.

Numa pesquisa que relaciona o stress a frequência cardíaca. Considera-se como a variável independente o stress e a variável dependente a frequência cardíaca;

Se quisermos relacionar o efeito da educação sobre o rendimento para medir o efeito do nível de escolaridade sobre a renda anual, a variável independente é o nível de escolaridade e a variável dependente é a renda anual;

Numa pesquisa sobre rendimento e horas de trabalho, nossa variável dependente é também o rendimento e a variável independente as horas de trabalho;

etc.

Variáveis dependentes ás vezes são limitadas em seu domínio e nem sempre podemos observar os dados sobre as variáveis dependentes e independentes ao longo de toda uma população. Uma variável limitada pode ser classificada em duas categorias, variável truncada ou variável com censura. Também definimos amostras como truncadas ou com censura dependendo da natureza da limitação da variável dependente pesquisada.

Amostra com Censura

Dizemos que uma amostra é censurada quando algumas observações sobre a variável dependente, correspondentes a valores conhecidos da(s) variável(is) independente(s), não são observados. Ou seja, quando dados sobre a variável resposta não estão completamente disponíveis para algumas unidades da amostra. No entanto, para estas unidades, os dados sobre as variáveis regressoras são totalmente conhecidos. Assim, temos:

C1: $\ y=y^{*}.1 (y^{*}>c)+c.1(y^{*}\leq c);$
C2: $\ y=y^{*}.1 (y^{*}<d)+d.1(y^{*}\geq d);$
C3: $\ y=y^{*}.1(c<y^{*}<d)+c.1(y^{*}\leq c)+d.1(y^{*}\geq d).$

Amostra Truncada

Neste caso, os valores das variáveis independentes são conhecidos somente quando a variável dependente é observada. De outra forma, a amostragem truncada geralmente surge quando uma pesquisa tem como alvo um subconjunto específico da população e ignora inteiramente a outra parte da população.

T1: $\ y=y^{*},\ \textrm{se}\ y^{*}>c,$ não observada caso contrário;
T2: $\ y=y^{*},\ \textrm{se}\ y^{*}<d,$ não observada caso contrário;
T3: $\ y=y^{*},\ \textrm{se}\ c<y^{*}<d,$ não observada caso contrário;

Densidade de uma v.a. truncada

Proposição: Seja $Y$ uma variável aleatória continua com função densidade de probabilidade $f(y)$ e $``a"$ constante real. A densidade condicional de $Y$ dado que $Y>a,$ fica dada por: $\begin{eqnarray}\label{prop1} g(y|Y>a)=\frac{f(y)}{1-F(a)}, \end{eqnarray}$

$F(.),$ função de distribuição acumulada de $Y.$

Se $Y\sim N(\mu,\sigma^2),$ então: $P(Y>a)=1-\Phi\left(\frac{a-\mu}{\sigma}\right)=1-\Phi(b)$ onde $b=\dfrac{a-\mu}{\sigma}$ e $\Phi(.)$ é a função de distribuição acumulada da normal padrão.
Segue neste caso que, $\begin{eqnarray} g(y|Y>a)=\dfrac{f(y)}{1-\Phi(b)}=\dfrac{\frac{1}{\sigma}\phi(\frac{y-\mu}{\sigma})}{1-\Phi(b)}, \end{eqnarray}$ onde $\phi(.)$ é a densidade da normal padrão.

Momentos e Razão Inversa de Mills

$\begin{eqnarray} \nonumber E(Y|Y>a)&=&\displaystyle \mu+\sigma\lambda(b)\\ \nonumber Var(Y|Y>a)&=&\sigma^{2}[1-\delta(b)], \end{eqnarray}$

em que $b=\dfrac{a-\mu}{\sigma}.$

$\begin{eqnarray} \nonumber \lambda(b)&=&\dfrac{\phi(b)}{1-\Phi(b)}\ \textrm{se o truncamento é}\ y>a, \\ \nonumber \lambda(b)&=&-\dfrac{\phi(b)}{\Phi(b)}\ \textrm{se o truncamento é}\ y<a, \end{eqnarray}$

e $\delta(b)=\lambda(b)[\lambda(b)-b].$

$\lambda(.)$ é conhecida como razão inversa de Mills.

Truncamento Incidental

Proposição: Sejam $Y$ e $Z$ duas variáveis aleatórias continuas com distribuição bivariada, correlação $\rho,$ e função densidade conjunta denotada por $f(y,z).$ A densidade de $Y,$ dado que $Z>a$ é:

$\begin{eqnarray} g(y|Z>a)=\frac{f(y,z)}{P(Z>a)}, a\in \mathcal{R} \ \textrm{qualquer.} \end{eqnarray}$

Momentos

Se $Y$ e $Z$ possuem distribuição normal bivariada com médias $\mu_{y}, \mu_{z},$ desvios-padrão $\sigma_{y}, \sigma_{z},$ respectivamente, e coeficiente de correlação $\rho,$ então,

$\begin{eqnarray} \displaystyle E(Y|Z>a)&=&\mu_{y}+\rho\sigma_{y}\lambda(b_{z})\\ \nonumber \displaystyle Var(Y|Z>a)&=&\sigma^{2}_{y}[1-\rho^{2}\delta(b_{z})], \end{eqnarray}$

com $b_{z}=\dfrac{a-\mu_{z}}{\sigma_{z}},\ \lambda(b_{z})=\dfrac{\phi(b_{z})}{1-\Phi(b_{z})}$ e $\delta(b_{z})=\lambda(b_{z})[\lambda(b_{z})-b_{z}]$

Modelo de Regressão Truncado

Considere uma variável aleatória $y_{i}^{*}$ linearmente dependente de $\mathbf{x}_{i},$ isto é:

$y_{i}^{*}=\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\epsilon_{i},\ \textrm{onde,}$

$\epsilon_{i}|\mathbf{x}_{i}\overset{iid}{\sim} N(0,\sigma^2)\Rightarrow y_{i}^{*}|x_{i}\overset{iid}{\sim} N(\mathbf{x}_{i}^{'}\boldsymbol{\beta},\sigma^{2})\ \textrm{e}\ E(y_{i}^{*})=\mathbf{x}_{i}^{'}\boldsymbol{\beta}.$

A observação

$i$ é observada somente se

$y_{i}^{*}$ está acima de um certo limite conhecido, isto é,

$\begin{eqnarray} y_{i}=\left\{ \begin{array}{cc} y_{i}^{*},& \textrm{se}\ y_{i}^{*}> a,\\ &\\ n.a,& \textrm{se}\ y_{i}^{*}\leq a \end{array} \right. \end{eqnarray}$

A função densidade de $y_{i}$ é dada por:

$\begin{eqnarray} \displaystyle f(y_{i}|x_{i})=f(y_{i}^{*}|y_{i}^{*}>a,x_{i})=\frac{1}{\sigma}\dfrac{\phi\left(\dfrac{y_{i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)}{\Phi\left(\frac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}-a}{\sigma}\right)} \end{eqnarray}$ Note que o valor esperado da variável observada é não linear em

$\mathbf{x}_{i}.$

$\begin{eqnarray} \nonumber E(y_{i}|\mathbf{x}_{i})&=&E(y_{i}^{*}|y_{i}^{*}>a,\mathbf{x}_{i})\\ \nonumber &=&\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\sigma\dfrac{\phi\left[(a-\mathbf{x}_{i}^{'}\boldsymbol{\beta})/\sigma\right]}{1-\Phi\left[(a-\mathbf{x}_{i}^{'}\boldsymbol{\beta})/\sigma\right]}\\ \nonumber &=&\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\sigma\lambda(b_{i}) \end{eqnarray}$

Estimação

A regressão linear simples da variável observada $y_{i}$ em $x_{i},$

$y_{i}=\mathbf{x}_{i}^{'}\boldsymbol{\beta}+u_{i},$

irá produzir estimativas viesadas de $\boldsymbol{\beta},$ pois o erro $u_{i}=(\epsilon_{i}|y_{i}^{*}>a)$ é correlacionado com $\mathbf{x}_{i}$ e $E(u_{i})=E(\epsilon_{i}|y_{i}^{*}>a)=\sigma\lambda(b_{i})>0.$

Assim, a regressão truncada é, em geral, estimada por máxima verossimilhança.

A função log de verossimilhança é

$\begin{eqnarray} \nonumber \mathcal{L}(\theta|\mathbf{x}_{i})=\displaystyle \sum_{i=1}^{n}\ln \left[\sigma^{-1}\phi\left(\dfrac{y_{i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right]-\sum_{i=1}^{n} \ln\left[1-\Phi\left(\dfrac{a-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right], \end{eqnarray}$

permitindo estimar tanto $\boldsymbol{\beta}$ quanto $\sigma$ por um procedimento numérico iterativo. Propriedades usuais de MV (consistência, eficiência assintótica, e etc) se aplicam.

Interpretação dos parâmetros

A interpretação dos parâmetros depende muito da questão de pesquisa. O efeito marginal sobre a variável dependente observada é dado por:

$\begin{eqnarray} \nonumber \dfrac{\partial E(y_{i}|x_{i})}{\partial \mathbf{x}_{ij}}&=&\dfrac{\partial E(y_{i}^{*}|y_{i}^{*}>a,x_{i})}{\partial \mathbf{x}_{ij}}\\ \nonumber &=&\boldsymbol{\beta}_{j}\left(1-\lambda_{i}^{2}+\alpha_{i}\lambda_{i}\right)\\ \nonumber &=& \boldsymbol{\beta}_{j}\underbrace{(1-\delta(\alpha_{i}))}_{0<\delta<1} \end{eqnarray}$

onde $\alpha_{i}=\left(\dfrac{a-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)$ e $\lambda_{i}=\lambda(\alpha_{i}).$

Modelo Tobit (Modelo de Regressão com Censura)

Considere uma variável aleatória $y_{i}^{*}$ linearmente dependente de $\mathbf{x}_{i},$ isto é:

$y_{i}^{*}=\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\epsilon_{i},\ \textrm{onde,}$

Suponha que o valor observado $y_{i}$ é censurado abaixo do zero, isto é,

$\begin{eqnarray} y_{i}=\left\{ \begin{array}{cc} y_{i}^{*},& \textrm{se}\ y_{i}^{*}> 0,\\ &\\ 0,& \textrm{se}\ y_{i}^{*}\leq 0 \end{array} \right. \end{eqnarray}$

A variável observada é uma variável aleatória mistura com massa de probabilidade em zero, $P(y_{i}=0|\mathbf{x}_{i})=P(y_{i}^{*}\leq 0|\mathbf{x}_{i})=\Phi\left(-\dfrac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)$ e para valores acima do zero com densidade: $f(y_{i}|x_{i})=\dfrac{1}{\sigma}\phi\left(\dfrac{y_{i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)$

O valor esperado da variável observada é,

$\begin{eqnarray} \nonumber E(y_{i}|\mathbf{x}_{i})&=&0.P(y_{i}^{*}\leq 0|x_{i})+E(y_{i}^{*}|y_{i}^{*}>0,\mathbf{x}_{i}).P(y_{i}^{*}> 0|x_{i})\\ %\nonumber &=&\left[\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\sigma\dfrac{\phi\left(\mathbf{x}_{i}^{'}\boldsymbol{\beta}/\sigma\right)}{\Phi\left(\mathbf{x}_{i}^{'}\boldsymbol{\beta}/\sigma\right)}\right]\Phi\left(\mathbf{x}_{i}^{'}\boldsymbol{\beta}/\sigma\right)\\ \nonumber &=&\Phi\left(\dfrac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\left[\mathbf{x}_{i}^{'}\boldsymbol{\beta}+\sigma\lambda\left(-\dfrac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right] \end{eqnarray}$

Estimação

Log de Verossimilhança:

$\begin{eqnarray} \nonumber \mathcal{L}(\theta|\mathbf{x}_{i})&=&\displaystyle \sum_{i|y_{i}>0} \ln\left[\sigma^{-1}\phi\left(\dfrac{y_{i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right]+\sum_{i|y_{i}=0}\ln\left[1-\Phi\left(\frac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right], \end{eqnarray}$

$\theta$ vetor de parâmetros.

Interpretação dos Parâmetros

A interpretação dos parâmetros depende muito da questão de pesquisa. Se o pesquisador está interessado no efeito marginal sobre a variável dependente latente

$(y_{i}^{*})$ ,

$\begin{eqnarray} \dfrac{\partial E(y_{i}^{*}|x_{i})}{\partial \mathbf{x}_{i,j}}=\boldsymbol{\beta_{j}}\ \textrm{ou}\ \dfrac{\partial E(y_{i}^{*}|y_{i}^{*}>0,x_{i})}{\partial \mathbf{x}_{i,j}}=\boldsymbol{\beta_{j}}[1-\lambda(b)^{2}-\lambda(b)] \end{eqnarray}$ No entanto, se o pesquisador está interessado no efeito marginal sobre o valor esperado dos valores observados, o efeito marginal é dado por,

$\begin{eqnarray} \dfrac{\partial E(y_{i}|x_{i})}{\partial \mathbf{x}_{i,j}}=\boldsymbol{\beta_{j}}\Phi\left(\dfrac{\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right) \end{eqnarray}$

Um dos problemas da metodologia Tobit é considerar que a equação que nos diz se uma observação está no limite é a mesma equação que nos diz o valor da variável dependente. Isso nem sempre é real. Temos então o modelo de Heckman, também chamado de modelo Tobit tipo II ou modelo de seleção amostral.

Modelo de Heckman ou Tobit tipo II

O modelo de seleção amostral é representado pelo seguinte sistema de regressão:

$\begin{eqnarray} \nonumber y_{1,i}^{*} &=& x_{i}^{T}\beta +u_{1i}\\ \nonumber y_{2,i}^{*} &=& z_{i}^{T}\gamma+u_{2i}, \end{eqnarray}$

onde as respostas $y_{1i}^{*}$ e $y_{2i}^{*}$ são variáveis latentes, $x_{i}$ e $z_{i}$ são as covariáveis, $\beta \in \mathcal{R}^{p+1}$ e $\gamma\in \mathcal{R}^{q+1}$ são vetores de parâmetros, e os termos de erro seguem uma distribuição normal bivariada,

$\left( \begin{array}{c} u_{1i} \\ u_{2i} \end{array} \right)\sim N\left\{\left( \begin{array}{c} 0 \\ 0 \end{array} \right),\left( \begin{array}{ll} \sigma^{2} & \rho\sigma \\ \rho\sigma & 1 \end{array} \right)\right\},$

com variância $\sigma_{2}^{2}=1, \sigma_{1}^{2}=\sigma^{2}$ e correlação $\rho.$ $\sigma_{2}^{2}=1$ para assegurar a identificabilidade.

Observa-se uma indicadora $\delta_{i}$ quando a variável latente $y_{2i}^{*}$ é positiva e definimos o valor da variável $y_{i}=y_{1,i}^{*}$ se a indicadora é 1. Ou seja, assim podemos fazer:

$\begin{eqnarray} \nonumber \delta_{i}=\left\{ \begin{array}{cc} 1,& \textrm{se}\ y_{2,i}^{*}>0,\\ &\\ 0,& \textrm{se}\ y_{2,i}^{*}\leq 0 \end{array} \right. \end{eqnarray}$

$\begin{eqnarray} \nonumber y_{i}=\left\{ \begin{array}{cc} y_{1,i}^{*}, & \textrm{se}\ \delta_{i}=1,\\ &\\ \textrm{0},& \textrm{se}\ \delta_{i}=0 \end{array} \right. \end{eqnarray}$ Considerando o modelo descrito acima temos:

$\begin{eqnarray} E(y_{i}|y_{1,i}^{*}\ \textrm{é observado},\mathbf{x}_{i},\mathbf{z}_{i})&=& \mathbf{x}_{i}^{'}\boldsymbol{\beta} + \boldsymbol{\beta_{\lambda}}\left[\frac{\phi(\mathbf{z}_{i}^{'}\boldsymbol{\gamma})}{\Phi(\mathbf{z}_{i}^{'}\boldsymbol{\gamma})}\right] \end{eqnarray}$ e,

$\begin{eqnarray} Var(y_{i}|y_{1,i}^{*}\ \textrm{é observado},\mathbf{x}_{i},\mathbf{z}_{i})&=&\sigma[1-\rho^{2}\delta(-\mathbf{z}_{i}^{'}\boldsymbol{\gamma})] \end{eqnarray}$

Se os erros são independentes, então $E(y_{i}|y_{1,i}^{*}\ \textrm{é observado},\mathbf{x}_{i},\mathbf{z}_{i})=\mathbf{x}_{i}^{'}\boldsymbol{\beta}$ e o método de mínimos quadrados ordinários (MQO) nos dão estimativas consistentes de $\boldsymbol{\beta}$ . No entanto, qualquer correlação entre os dois erros significa que a média truncada não é mais $\mathbf{x}_{i}^{'}\boldsymbol{\beta}$ e o MQO produz estimativas viesadas de $\boldsymbol{\beta},$ pois o fator $\rho\sigma\left[\frac{\phi(\mathbf{z}_{i}^{'}\boldsymbol{\gamma})}{\Phi(\mathbf{z}_{i}^{'}\boldsymbol{\gamma})}\right]$ é omitido e torna-se parte do termo de erro $u_{1,i}.$ O viés resultante é chamado de viés de seleção amostral ou, somente, viés de seleção.

Estimação

Método de Máxima Verossimilhança

Utilizando o pressuposto de normalidade bivariada, a função de verossimilhança é dada por:

$\begin{eqnarray} \mathcal{L}(\theta|\mathbf{x}_{i},\mathbf{z}_{i})&=&\displaystyle \sum_{\delta_{i}=0} \ln \left[ \Phi\left\{ \dfrac{\mathbf{z}_{i}^{'}\boldsymbol{\gamma}+\frac{\rho}{\sigma} (y_{1,i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta})}{\sqrt{1-\rho^{2}}}\right\}\right]\\ \nonumber &+&\sum_{\delta_{i}=1} \ln \left\{ \sigma^{-1}\phi\left(\dfrac{y_{1,i}-\mathbf{x}_{i}^{'}\boldsymbol{\beta}}{\sigma}\right)\right\}\\ \nonumber &+&\sum_{\delta_{i}=1} \ln \left(\Phi(-\mathbf{z}_{i}^{'}\boldsymbol{\gamma})\right) \end{eqnarray}$ onde

$\theta$ é um vetor de parâmetros.

O estimador obtido é consistente, assintoticamente normal e eficiente, mas tem vários inconvenientes. Primeiro de tudo é não-linear e, obviamente, requer métodos numéricos iterativos. É muito caro do ponto de vista computacional, apesar de ser possível implementa-lo, existe o problema da não linearidade, a função de verossimilhança também tem máximos locais (Olsen 1982), o que requer um bom ponto de partida para o algoritmo numérico.

Dois passos de Heckman

Heckman propôs um processo simples de duas etapas que envolve a estimativa de um probit padrão e um modelo de regressão linear. Os dois passos são:

Estima-se a equação probit por máxima verossimilhança para obter uma estimativa de $\boldsymbol{\gamma}.$ Para cada observação da amostra selecionada, computamos: $\begin{eqnarray} \hat{\lambda}=\frac{\phi(\mathbf{z}_{i}^{'}\hat{ \boldsymbol{\gamma}})}{ \Phi(\mathbf{z}_{i}^{'}\hat{\boldsymbol{\gamma}})}\ \textrm{e}\ \hat{\delta}_{i}=\hat{\lambda}_{i}(\hat{\lambda}_{i}-\mathbf{z}_{i}^{'} \hat{\boldsymbol{\gamma}}). \end{eqnarray}$

Estima-se $\boldsymbol{\beta},\ \boldsymbol{\beta}_{\lambda}=\rho\sigma, \rho$ e $\sigma$ por mínimos quadrados a partir da equação: $\begin{eqnarray} \displaystyle y_{i}=\mathbf{x}_{i}^{'}\boldsymbol{\beta} + \boldsymbol{\beta_{\lambda}}\hat{\lambda_{i}}+\varepsilon_{i} \end{eqnarray}$

Vejamos alguns Exemplos:

MROZ 1987

De acordo com Hill, Griffiths, and Lim (2008), em um problema de estimação de parâmetros, se os dados são obtidos por amostragem aleatória, métodos de estimação, clássicos, tais como mínimos quadrados, funcionam bem. No entanto, se os dados são obtidos por um procedimento de amostragem que não é aleatório, os procedimentos normais não funcionam adequadamente. Geralmente enfrentamos tais problemas na amostragem. Uma ilustração famosa vem de economia do trabalho. Se quisermos estudar os determinantes dos salários das mulheres casadas, nos deparamos com um problema de seleção da amostra. Pois, ao recolher dados sobre as mulheres casadas, e perguntar-lhes o salário que ganham, muitas vão responder que são donas de casa. Nós só observamos os dados sobre salários de mercado quando a mulher escolhe entrar na força de trabalho. Uma estratégia é ignorar as mulheres que são donas de casa, omiti-las a partir da amostra, em seguida, usar mínimos quadrados para estimar a equação de salários para aquelas que trabalham. Esta estratégia é falha, a razão para a falha é que a amostra não é uma amostra aleatória, uma vez que os dados observados, ao omitir as mulheres que são donas de casa, são selecionados por um processo sistemático.

O primeiro exemplo é retirado do livro de Hill, Griffiths, and Lim (2008), vamos reconsiderar a análise dos salários recebidos por mulheres casadas usando o conjunto de dados obtidos de Mroz (1987). Na amostra de 753 mulheres casadas, 428 são empregadas no mercado formal e recebem salários diferentes de zero. Vejam uma parte do conjunto de dados:

setwd("~/GitHub/Modelos de Heckman/Modelo-de-Heckman")
#install.packages("Rtools")
#devtools::install_github("hadley/readxl")
library(readxl)
dados <- read_excel("~/GitHub/Modelos de Heckman/Modelo-de-Heckman/mroz_z.xls",col_names = TRUE, col_types = NULL)

## DEFINEDNAME: 0b 00 00 11 02 00 00 00 00 00 00 00 00 00 00 5f 78 6c 66 6e 2e 4e 4f 52 4d 2e 53 2e 44 49 53 54 1c 1d 
## DEFINEDNAME: 0b 00 00 11 02 00 00 00 00 00 00 00 00 00 00 5f 78 6c 66 6e 2e 4e 4f 52 4d 2e 53 2e 44 49 53 54 1c 1d 
## DEFINEDNAME: 0b 00 00 11 02 00 00 00 00 00 00 00 00 00 00 5f 78 6c 66 6e 2e 4e 4f 52 4d 2e 53 2e 44 49 53 54 1c 1d 
## DEFINEDNAME: 0b 00 00 11 02 00 00 00 00 00 00 00 00 00 00 5f 78 6c 66 6e 2e 4e 4f 52 4d 2e 53 2e 44 49 53 54 1c 1d

dados[423:433,c(5,6,19,23,14,1,7,21)]

## # A tibble: 11 × 8
##      age  educ exper  kids    mtr  inlf   wage     lwage
##    <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>     <dbl>
## 1     32    17    14     1 0.7215     1 4.7826 1.5649840
## 2     36    10     2     1 0.7215     1 2.3118 0.8380265
## 3     40    12    21     1 0.6215     1 5.3061 1.6688570
## 4     43    13    22     1 0.5815     1 5.8675 1.7694290
## 5     33    12    14     1 0.5815     1 3.4091 1.2264480
## 6     30    12     7     1 0.6915     1 4.0816 1.4064890
## 7     49    12     2     1 0.6615     0     NA        NA
## 8     30    16     5     1 0.6615     0     NA        NA
## 9     30    12    12     1 0.6615     0     NA        NA
## 10    41    12     1     1 0.5815     0     NA        NA
## 11    45    12    12     1 0.6615     0     NA        NA

attach(dados)
#tail(MEPS2001)

par(mfrow=c(1,2))
#library(stringr)
#dados$wage<-as.numeric(dados$wage)
#wage2<- str_replace(dados$wage, pattern="NA", replacement= 0)
#wage2<-as.numeric(wage2)
#dados$wage2
#dados$lwage<-as.numeric(dados$lwage)
#lwage2<- str_replace(dados$lwage, pattern="NA", replacement= 0)
#lwage2<-as.numeric(lwage2)
#dados$lwage2
hist(wage,ylim=c(0,200),xlim = c(-1,30),xlab = "Salário de mulheres com emprego formal", ylab = "Frequência",main = "Dados MROZ")
hist(lwage,ylim=c(0,150),xlim = c(-3,5),xlab = "Log do Salário de mulheres com emprego formal", ylab = "Frequência",main = "Dados MROZ")

Primeiro, vamos estimar uma equação simples para o salário, explicando $ln(salário)$ como uma função da educação, EDUC, e anos de experiência no mercado de trabalho (EXPER), utilizando as 428 mulheres que têm salários positivos. O modelo é:

$ln(wage)_{i}=\beta_{0}+\beta_{1}Educ+\beta_{2}exper+\epsilon_{i},$

com, $\epsilon\sim N(0,1).$

as estimativas dos parâmetros são:

lreg<-lm(lwage~educ+exper,data = dados[dados$inlf==1,])
summary(lreg)

## 
## Call:
## lm(formula = lwage ~ educ + exper, data = dados[dados$inlf == 
##     1, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.05608 -0.30524  0.05599  0.38846  2.27384 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.400174   0.190368  -2.102 0.036132 *  
## educ         0.109489   0.014167   7.728 7.94e-14 ***
## exper        0.015674   0.004019   3.900 0.000112 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.669 on 425 degrees of freedom
## Multiple R-squared:  0.1484, Adjusted R-squared:  0.1444 
## F-statistic: 37.02 on 2 and 425 DF,  p-value: 1.512e-15

O retorno estimado para a educação é cerca de 11%, e os coeficientes estimados de educação e experiência são estatisticamente significativos.

Agora vamos ajustar o modelo:

$\begin{eqnarray*} lwage_{i} &=& \beta_{0} + \beta_{1}educ + \beta_{2}exper + u_{1} \end{eqnarray*}$

e assumimos que o log do salário $(lwage)$ é observado se:

$\begin{eqnarray*} \beta_{0}+\beta_{1}Age+\beta_{2}Educ_{i}+\beta_{3}kids_{i}+\beta_{4}mtr+u_{2}>0, \end{eqnarray*}$

onde $u_{1}$ e $u_{2}$ são correlacionados.

Ajuste com função Glm e lm do R

O procedimento Heckit (seleção amostral) começa por estimar um modelo probit de participação na força de trabalho. Como variáveis explicativas usamos a idade da mulher, seus anos de escolaridade, uma variável de indicador para saber se ela tem filhos, e a taxa de imposto marginal que ela iria pagar sobre rendimentos se empregada. O modelo probit é dado por:

$\begin{eqnarray*} inlf_{i} &\sim& Binomial (n, \pi_i)\\ g(\pi_i) &=& \beta_{0}+\beta_{1}Age+\beta_{2}Educ_{i}+\beta_{3}kids_{i}+\beta_{4}mtr \end{eqnarray*}$

Sendo $inlf$ a variável resposta binária, $Age_{i},Educ_{i},kids_{i}$ e $mtr_{i}$ as i-ésimas realizações das respectivas variáveis explicativas $Age,Educ,kids$ e $mtr,$ respectivamente, e $g(\pi_i)=\Phi^{-1}(\pi_{i})$ a função de ligação probit. Como houve indivíduos com o mesmo conjunto de covariáveis considera-se o número de repetições $n=428$ para a distribuição de $inlf_{i}.$

as estimativas dos parâmetros são:

#Modelo Probit
fit1<-glm(inlf~age+educ+kids+mtr,family=binomial(link=probit),data=dados)
summary(fit1)

## 
## Call:
## glm(formula = inlf ~ age + educ + kids + mtr, family = binomial(link = probit), 
##     data = dados)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.908  -1.208   0.816   1.060   1.588  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.192261   0.726243   1.642 0.100656    
## age         -0.020615   0.007031  -2.932 0.003369 ** 
## educ         0.083776   0.023120   3.623 0.000291 ***
## kids        -0.313885   0.122509  -2.562 0.010403 *  
## mtr         -1.393817   0.626166  -2.226 0.026017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1029.75  on 752  degrees of freedom
## Residual deviance:  988.29  on 748  degrees of freedom
## AIC: 998.29
## 
## Number of Fisher Scoring iterations: 4

Como esperado, os efeitos da idade, a presença de crianças, e as perspectivas de impostos mais altos reduzem significativamente a probabilidade de que uma mulher se junte à força de trabalho, enquanto a educação aumenta a probabilidade. Utilizando os coeficientes estimados calculamos a razão inversa de Mills para as 428 mulheres com salários de mercado.

library(sampleSelection)
dados$IMR <- invMillsRatio(fit1)$IMR1#Criar uma nova coluna com a covariável Razão Inversa de Mills 
dados[423:433,c(5,6,19,23,14,1,7,21,32)]

## # A tibble: 11 × 9
##      age  educ exper  kids    mtr  inlf   wage     lwage       IMR
##    <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>     <dbl>     <dbl>
## 1     32    17    14     1 0.7215     1 4.7826 1.5649840 0.4412366
## 2     36    10     2     1 0.7215     1 2.3118 0.8380265 0.8181494
## 3     40    12    21     1 0.6215     1 5.3061 1.6688570 0.6793267
## 4     43    13    22     1 0.5815     1 5.8675 1.7694290 0.6340350
## 5     33    12    14     1 0.5815     1 3.4091 1.2264480 0.5657426
## 6     30    12     7     1 0.6915     1 4.0816 1.4064890 0.6164304
## 7     49    12     2     1 0.6615     0     NA        NA 0.8290024
## 8     30    16     5     1 0.6615     0     NA        NA 0.4219251
## 9     30    12    12     1 0.6615     0     NA        NA 0.5929940
## 10    41    12     1     1 0.5815     0     NA        NA 0.6586579
## 11    45    12    12     1 0.6615     0     NA        NA 0.7763784

Este é então incluído na equação de regressão múltipla de salários e aplicado mínimos quadrados para obter as seguintes estimativas:

#Modelo de Regressão Linear Simples com wage>0
fit2<- lm(lwage~educ+exper+IMR, data = dados[dados$inlf==1,])
summary(fit2)

## 
## Call:
## lm(formula = lwage ~ educ + exper + IMR, data = dados[dados$inlf == 
##     1, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.02634 -0.28901  0.05418  0.37278  2.42408 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.810529   0.494476   1.639  0.10192    
## educ         0.058459   0.023850   2.451  0.01464 *  
## exper        0.016320   0.003998   4.082 5.34e-05 ***
## IMR         -0.866429   0.326988  -2.650  0.00836 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6643 on 424 degrees of freedom
## Multiple R-squared:  0.1622, Adjusted R-squared:  0.1563 
## F-statistic: 27.37 on 3 and 424 DF,  p-value: 3.381e-16

Notemos que o coeficiente estimado da razão inversa de Mills é estatisticamente significativo, o que implica que existe um viés de seleção presente nos resultados de quadrados mínimos da primeira equação de regressão sem a covariável IMR. Além disso, o retorno salárial estimado para a educação diminuiu de 11% para aproximadamente 6%.

Outra forma de proceder os calculos acima é utilizar o método de Máxima Verossimilhança ou o método de duas etapas de Heckman para estimar todos os parâmetros do modelo de regressão via pacote sampleSelection do R. Vejamos as saídas:

Ajuste com o pacote Sample Selection

Processo de duas etapas

library("sampleSelection")
# Two-step estimation
fit3<-heckit( inlf~age+educ+kids+mtr,
              lwage~educ+exper, dados)
summary(fit3)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 753 observations (325 censored and 428 observed)
## 11 free parameters (df = 743)
## Probit selection equation:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.192296   0.720544   1.655 0.098404 .  
## age         -0.020616   0.007045  -2.926 0.003534 ** 
## educ         0.083775   0.023205   3.610 0.000326 ***
## kids        -0.313885   0.123711  -2.537 0.011376 *  
## mtr         -1.393853   0.616575  -2.261 0.024070 *  
## Outcome equation:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.810542   0.610799   1.327 0.184910    
## educ        0.058458   0.029635   1.973 0.048915 *  
## exper       0.016320   0.004202   3.884 0.000112 ***
## Multiple R-Squared:0.1622,   Adjusted R-Squared:0.1563
##    Error terms:
##               Estimate Std. Error t value Pr(>|t|)  
## invMillsRatio  -0.8664     0.3993   -2.17   0.0303 *
## sigma           0.9326         NA      NA       NA  
## rho            -0.9291         NA      NA       NA  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Note neste caso a diferença entre as t-estatísticas. Os valores anteriores das t-estatísticas, calculados via função glm e lm são baseados em erros padrão como normalmente calculado pelo uso da regressão de mínimos quadrados. Os habituais erros padrão não levam em conta o fato de que a razão inversa de Mills é um valor estimado. Assim, o pacote sampleSelection corrige os erros padrão ao levar em conta a estimativa do primeiro estágio probit, estes são usados para construir as t-estatísticas ajustadas. Como podemos ver as t-estatísticas ajustadas são um pouco menores na saída do sampleSelection, indicando que os erros padrão ajustados são um pouco maiores do que os da saída do glm e lm.

É preferível estimar o modelo completo, tanto a equação de seleção e a equação de interesse, em conjunto por máxima verosimilhança. Pois, os erros padrão com base no procedimento de máxima verossimilhança são menores do que aqueles gerados pelo método de estimação de dois passos.

Método de Máxima Verossimilhança

# ML estimation
fit4<-selection(inlf~age+educ+kids+mtr,
                lwage~educ+exper, dados)
summary(fit4)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 6 iterations
## Return code 2: successive function values within tolerance limit
## Log-Likelihood: -913.561 
## 753 observations (325 censored and 428 observed)
## 10 free parameters (df = 743)
## Probit selection equation:
##              Estimate Std. error t value  Pr(> t)    
## (Intercept)  1.595958   0.623731   2.559  0.01051 *  
## age         -0.013262   0.005939  -2.233  0.02555 *  
## educ         0.063931   0.021745   2.940  0.00328 ** 
## kids        -0.152592   0.099587  -1.532  0.12546    
## mtr         -2.291885   0.537565  -4.263 2.01e-05 ***
## Outcome equation:
##             Estimate Std. error t value Pr(> t)    
## (Intercept) 0.668587   0.235006   2.845 0.00444 ** 
## educ        0.065816   0.016635   3.957 7.6e-05 ***
## exper       0.011767   0.004094   2.875 0.00404 ** 
##    Error terms:
##       Estimate Std. error t value Pr(> t)    
## sigma  0.84944    0.04248   20.00  <2e-16 ***
## rho   -0.83947    0.03490  -24.05  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Neste caso a estimativa do parâmetro da razão inversa de Mills é:

IMR<-coef(summary(fit4))["rho",1]*coef(summary(fit4))["sigma",1]
IMR

## [1] -0.7130812

#res<-fit4$estimate
#IMR<-res[9]*res[10]
#IMR

Ajuste com o pacote SSMROB

library(ssmrob)
#Equação de seleção que será ajustada via modelo probit
 selectEq <- inlf~age+educ+kids+mtr
#Equação de regressão que será ajustada via modelo linear simples
outcomeEq <- lwage~educ+exper
fit5<-heckitrob(outcomeEq,selectEq,control=heckitrob.control(tcc=3.2,weights.x1="robCov"))
summary(fit5)

## ------------------------------------------------------------- 
## Robust 2-step Heckman / heckit M-estimation 
## Probit selection equation: 
##                  Estimate   Std.Error t-value  p-value    
## XS(Intercept)  1.19229732 0.726234392   1.642 0.101000    
## XSage         -0.02061554 0.007031349  -2.932 0.003370  **
## XSeduc         0.08377532 0.023119989   3.624 0.000291 ***
## XSkids        -0.31388480 0.122508531  -2.562 0.010400   *
## XSmtr         -1.39385431 0.626155191  -2.226 0.026000   *
## Outcome equation: 
##                  Estimate   Std.Error t-value  p-value    
## XO(Intercept)  0.77909013 0.596402209   1.306 1.91e-01    
## XOeduc         0.06440744 0.028843743   2.233 2.56e-02   *
## XOexper        0.01585891 0.003503069   4.527 5.98e-06 ***
## imrData$IMR1  -0.86768171 0.393821292  -2.203 2.76e-02   *
## ---
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
## -------------------------------------------------------------

MEPS 2001 Sem a Covariável Rendimento

O artigo de Marchenko and Genton (2012) considerou os dados sobre os gastos ambulatoriais do Medical Expenditure Panel Survey 2001 (MEPS2001), analisadas por Cameron e Trivedi (2010). MEPS é a fonte de dados sobre o custo e a utilização de cuidados de saúde e cobertura de seguro de saúde mais completa dos Estados Unidos, segundo a Agência de Investigação de Saúde e Qualidade (AHRQ) dos EUA. A amostra é restrita a apenas aqueles indivíduos que estão cobertos por seguros privados, com idades entre 21 e 64 anos. Os dados consistem de 3328 observações, dos quais 526 (15,8%) correspondem aos valores de despesas zero. O conjunto de dados inclui diversas variáveis explicativas, tais como idade, sexo, anos de escolaridade, entre outros.

library(ssmrob)
data(MEPS2001)
attach(MEPS2001)
#head(MEPS2001)
MEPS2001[1:10,c(1,2,4,8,18,22,20,16,17,21)]

##    educ age female totchr blhisp ins dambexp ambexp   lambexp    lnambx
## 1    11 3.3      0      2      0   0       1    760 6.6333184 6.6333184
## 2    14 2.1      0      0      1   0       1    497 6.2085900 6.2085900
## 3    12 4.0      1      1      1   0       1   1002 6.9097533 6.9097533
## 4    14 5.2      1      0      1   1       1    745 6.6133842 6.6133842
## 5    16 5.0      1      0      0   0       1   2728 7.9113240 7.9113240
## 6    12 3.7      0      0      0   0       1    636 6.4551988 6.4551988
## 7    14 3.2      1      0      0   0       1      2 0.6931472 0.6931472
## 8    12 4.6      0      0      0   1       0      0        NA 0.0000000
## 9    12 4.0      1      0      1   1       0      0        NA 0.0000000
## 10   17 6.4      0      0      0   0       1   2826 7.9466176 7.9466176

#tail(MEPS2001)

Vejam a distribuição dos dados brutos de despesas ambulatoriais e dos dados de despesas ambulatoriais logaritmados:

library(ssmrob)
#Carregando o conjunto de dados MEPS2001 - dados de despesas ambulatoriais
data(MEPS2001)
attach(MEPS2001)
par(mfrow=c(1,2))
hist(ambexp,ylim = c(0,3500),xlim=c(0,20000) ,xlab = "Despesas Ambulotariais", ylab = "Frequência",main = "Dados do MEPS 2001")
hist(lnambx,ylim = c(0,800),xlim=c(0,12), xlab = "Log das Despesas Ambulotariais", ylab = "Frequência",main = "Dados do MEPS 2001")

Antes de dar sequência com a aplicação dos modelos de seleção amostral, considere uma regressão linear múltipla relacionando os gastos ambulotariais com as diversas outras variáveis disponíveis no banco de dados e supondo erro aleatório com distribuição normal, ou seja, considere o modelo:

$lnambx_{i}=\beta_{0}+\beta_{1}Age+\beta_{2}female+\beta_{3}Educ+\beta_{4}blhisp+\beta_{5}totchr+\beta_{6}ins+\epsilon_{i},$

com, $\epsilon\sim N(0,1).$

as estimativas dos parâmetros são:

lreg<-lm(lnambx~age+female+educ+blhisp+totchr+ins)
summary(lreg)

## 
## Call:
## lm(formula = lnambx ~ age + female + educ + blhisp + totchr + 
##     ins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8228 -0.9666  0.5640  1.5731  5.8482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.72876    0.28126   6.147 8.86e-10 ***
## age          0.32473    0.03835   8.468  < 2e-16 ***
## female       1.14470    0.08334  13.735  < 2e-16 ***
## educ         0.11411    0.01654   6.898 6.27e-12 ***
## blhisp      -0.73418    0.09289  -7.904 3.64e-15 ***
## totchr       1.05940    0.05537  19.133  < 2e-16 ***
## ins          0.20783    0.08691   2.391   0.0168 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.381 on 3321 degrees of freedom
## Multiple R-squared:  0.2346, Adjusted R-squared:  0.2332 
## F-statistic: 169.7 on 6 and 3321 DF,  p-value: < 2.2e-16

Podemos ver que todos os parâmetros são altamente significativos, esse ajuste foi realizado com os dados completos incluindo os valores zero. Poderíamos pensar que estes não devem fazer parte do ajuste e retirá-los, neste caso teríamos uma mudança significativa no ajuste:

lreg<-lm(lnambx~age+female+educ+blhisp+totchr+ins,data = MEPS2001[ MEPS2001$dambexp == 1, ] )
summary(lreg)

## 
## Call:
## lm(formula = lnambx ~ age + female + educ + blhisp + totchr + 
##     ins, data = MEPS2001[MEPS2001$dambexp == 1, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0471 -0.7959  0.0818  0.8644  4.6335 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.907825   0.168151  29.187  < 2e-16 ***
## age          0.217233   0.022223   9.775  < 2e-16 ***
## female       0.379376   0.048577   7.810 8.04e-15 ***
## educ         0.022239   0.009761   2.278   0.0228 *  
## blhisp      -0.238532   0.055195  -4.322 1.60e-05 ***
## totchr       0.561817   0.030508  18.416  < 2e-16 ***
## ins         -0.020827   0.050006  -0.416   0.6771    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.27 on 2795 degrees of freedom
## Multiple R-squared:  0.1918, Adjusted R-squared:  0.1901 
## F-statistic: 110.6 on 6 and 2795 DF,  p-value: < 2.2e-16

Mas o fato é que não podemos deixar de considerar as pessoas que possuem gasto zero com despesas ambulatoriais, pois do contrário, nossa amostra não seria obtida de forma aleatória, ou seja, caso não consideremos essa parte da amostragem teríamos um problema de seleção amostral. Assim, a disposição para gastar será relacionada com algumas covariáveis por meio de um modelo probit e posteriormente iremos criar uma nova covariável (razão inversa de Mills) e iremos analisar o modelo de regressão de interesse acrescido desta variável. Dessa forma, vamos utilizar o modelo de seleção amostral clássico de Heckman para analisar esses dados. Como a distribuição dos gastos é altamente viesada, a análise foi realizada utilizando a escala logarítmica.

Agora vamos ajustar o modelo:

$\begin{eqnarray*} lnambx_{i} &=& \beta_{0}+\beta_{1}Age+\beta_{2}female_{i}+\beta_{3}educ_{i}+\beta_{4}blhisp+\beta_{5}totchr+\beta_{6}ins + u_{1} \end{eqnarray*}$

e assumimos que o log das despesas ambulatoriais $(lnambx)$ é observado se:

$\begin{eqnarray*} \beta_{0}+\beta_{1}Age+\beta_{2}female_{i}+\beta_{3}educ_{i}+\beta_{4}blhisp+\beta_{5}totchr+\beta_{6}ins+u_{2}>0, \end{eqnarray*}$

onde $u_{1}$ e $u_{2}$ são correlacionados.

Ajuste com as funções glm e lm do R

Primeiro analisamos os dados utilizando as funções glm e lm do R.

O modelo probit utlizado foi:

$\begin{eqnarray*} dambexp_{i} &\sim& Binomial (n, \pi_i)\\ g(\pi_i) &=& \beta_{0}+\beta_{1}Age+\beta_{2}female_{i}+\beta_{3}educ_{i}+\beta_{4}blhisp+\beta_{5}totchr+\beta_{6}ins \end{eqnarray*}$

Os parâmetros estimados são:

library(aod)
#Modelo probit
fit1<-glm( dambexp ~ age+female+educ+blhisp+totchr+ins, family = binomial(link = "probit"),data=MEPS2001)
summary(fit1)

## 
## Call:
## glm(formula = dambexp ~ age + female + educ + blhisp + totchr + 
##     ins, family = binomial(link = "probit"), data = MEPS2001)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4495   0.1284   0.3873   0.6262   1.5804  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.71770    0.19086  -3.760 0.000170 ***
## age          0.09731    0.02703   3.601 0.000317 ***
## female       0.64421    0.06022  10.698  < 2e-16 ***
## educ         0.07017    0.01123   6.247 4.19e-10 ***
## blhisp      -0.37449    0.06156  -6.083 1.18e-09 ***
## totchr       0.79352    0.07140  11.113  < 2e-16 ***
## ins          0.18124    0.06227   2.911 0.003608 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2904.9  on 3327  degrees of freedom
## Residual deviance: 2395.3  on 3321  degrees of freedom
## AIC: 2409.3
## 
## Number of Fisher Scoring iterations: 6

Com os parâmetros do probit estimado foi possível encontrar a razão inversa de Mills.

MEPS2001$IMR <- invMillsRatio(fit1)$IMR1#Criar uma nova coluna com a covariável Razão Inversa de Mills 
MEPS2001[1:10,c(1,2,4,8,18,22,21,20,23)]

##    educ age female totchr blhisp ins    lnambx dambexp        IMR
## 1    11 3.3      0      2      0   0 6.6333184       1 0.05966001
## 2    14 2.1      0      0      1   0 6.2085900       1 0.73870733
## 3    12 4.0      1      1      1   0 6.9097533       1 0.12209510
## 4    14 5.2      1      0      1   1 6.6133842       1 0.21276541
## 5    16 5.0      1      0      0   0 7.9113240       1 0.13082628
## 6    12 3.7      0      0      0   0 6.4551988       1 0.51722607
## 7    14 3.2      1      0      0   0 0.6931472       1 0.21318779
## 8    12 4.6      0      0      0   1 0.0000000       0 0.38796767
## 9    12 4.0      1      0      1   1 0.0000000       0 0.30092333
## 10   17 6.4      0      0      0   0 7.9466176       1 0.25274431

Assim, a covariável $IMR$ foi utilizada no ajuste da regressão de interesse:

fit2 <- lm( lnambx ~ age+female+educ+blhisp+totchr+ins + IMR,
                data = MEPS2001[ MEPS2001$dambexp == 1, ] )
summary(fit2)

## 
## Call:
## lm(formula = lnambx ~ age + female + educ + blhisp + totchr + 
##     ins + IMR, data = MEPS2001[MEPS2001$dambexp == 1, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0498 -0.7840  0.0794  0.8689  4.6281 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.30257    0.29072  18.240  < 2e-16 ***
## age          0.20212    0.02400   8.422  < 2e-16 ***
## female       0.28916    0.07278   3.973 7.27e-05 ***
## educ         0.01199    0.01154   1.039  0.29871    
## blhisp      -0.18106    0.06509  -2.781  0.00545 ** 
## totchr       0.49833    0.04884  10.204  < 2e-16 ***
## ins         -0.04740    0.05248  -0.903  0.36647    
## IMR         -0.48017    0.28852  -1.664  0.09617 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.269 on 2794 degrees of freedom
## Multiple R-squared:  0.1926, Adjusted R-squared:  0.1906 
## F-statistic: 95.23 on 7 and 2794 DF,  p-value: < 2.2e-16

#wald.test(b = coef(fit3), Sigma = vcov(fit3), Terms = 8)

Notemos que o parâmetro que acompanha a covariável $IMR$ foi não significativo e poderíamos assumir que não há viés de seleção amostral e desconsiderar os valores iguais a zero para a variável gasto com despesas ambulotariais. O que implica que gastar com despesas ambulatoriais não está relacionado com a decisão de gastar e podem ser analisados separadamente por meio de Mínimos Quadrados Ordinários. Esta conclusão parece não plausível. Como observado por Cameron and Trivedi (2009), a suposição de normalidade dos erros é muito suspeita para esses dados. O que é possível de se verificar por meio da visualização do histograma do log das despesas apresentado anteriormente e devido a não significância da covariável $IMR$ .

Esse mesmo ajuste pode ser feito através do pacote sampleSelection usando o método de dois passos e máxima verossimilhança:

Ajuste com o pacote Sample Selection

No processo de duas etapas do pacote sample selection, manteve-se o resultado anterior com a não significância do parâmetro que acompanha a covariável $IMR.$

Processo de duas etapas

library("sampleSelection")
# Two-step estimation
fit3<-heckit( dambexp ~ age+female+educ+blhisp+totchr+ins,
                 lnambx ~ age+female+educ+blhisp+totchr+ins, MEPS2001 )
summary(fit3)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 3328 observations (526 censored and 2802 observed)
## 17 free parameters (df = 3312)
## Probit selection equation:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.71771    0.19247  -3.729 0.000195 ***
## age          0.09732    0.02702   3.602 0.000320 ***
## female       0.64421    0.06015  10.710  < 2e-16 ***
## educ         0.07017    0.01134   6.186 6.94e-10 ***
## blhisp      -0.37449    0.06175  -6.064 1.48e-09 ***
## totchr       0.79352    0.07112  11.158  < 2e-16 ***
## ins          0.18124    0.06259   2.896 0.003809 ** 
## Outcome equation:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.30257    0.29414  18.028  < 2e-16 ***
## age          0.20212    0.02430   8.319  < 2e-16 ***
## female       0.28916    0.07369   3.924 8.89e-05 ***
## educ         0.01199    0.01168   1.026    0.305    
## blhisp      -0.18106    0.06585  -2.749    0.006 ** 
## totchr       0.49833    0.04947  10.073  < 2e-16 ***
## ins         -0.04740    0.05315  -0.892    0.373    
## Multiple R-Squared:0.1926,   Adjusted R-Squared:0.1906
##    Error terms:
##               Estimate Std. Error t value Pr(>|t|)  
## invMillsRatio  -0.4802     0.2907  -1.652   0.0986 .
## sigma           1.2932         NA      NA       NA  
## rho            -0.3713         NA      NA       NA  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Mesmo resultado utilizando máxima verossimilhança.

Método de Máxima Verossimilhança

# ML estimation
fit4<-selection( dambexp ~ age+female+educ+blhisp+totchr+ins,
                    lnambx ~ age+female+educ+blhisp+totchr+ins, MEPS2001)
summary(fit4)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 4 iterations
## Return code 2: successive function values within tolerance limit
## Log-Likelihood: -5838.397 
## 3328 observations (526 censored and 2802 observed)
## 16 free parameters (df = 3312)
## Probit selection equation:
##             Estimate Std. error t value  Pr(> t)    
## (Intercept) -0.72444    0.19243  -3.765 0.000167 ***
## age          0.09845    0.02699   3.648 0.000264 ***
## female       0.64367    0.06014  10.703  < 2e-16 ***
## educ         0.07025    0.01134   6.195 5.85e-10 ***
## blhisp      -0.37263    0.06173  -6.036 1.58e-09 ***
## totchr       0.79467    0.07103  11.188  < 2e-16 ***
## ins          0.18212    0.06255   2.912 0.003595 ** 
## Outcome equation:
##             Estimate Std. error t value  Pr(> t)    
## (Intercept)  5.03743    0.22619  22.271  < 2e-16 ***
## age          0.21229    0.02296   9.247  < 2e-16 ***
## female       0.34973    0.05967   5.861 4.61e-09 ***
## educ         0.01887    0.01053   1.793 0.072971 .  
## blhisp      -0.21960    0.05948  -3.692 0.000222 ***
## totchr       0.54095    0.03906  13.848  < 2e-16 ***
## ins         -0.02954    0.05104  -0.579 0.562801    
##    Error terms:
##       Estimate Std. error t value Pr(> t)    
## sigma  1.27074    0.01821   69.77  <2e-16 ***
## rho   -0.12421    0.14438   -0.86    0.39    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Neste caso a estimativa do parâmetro da razão inversa de Mills é:

IMR<-coef(summary(fit4))["rho",1]*coef(summary(fit4))["sigma",1]
IMR

## [1] -0.1578385

#res<-fit4$estimate
#IMR<-res[15]*res[16]
#IMR

Poderíamos então assumir erroneamente que não há viés de seleção amostral, mesmo que essa suposição seja estranha dada a possivel relação das variáveis, mas foi o resultado apresentando pelos ajustes de dois passos de Heckman e máxima verossimilhança. Porém, outro pacote que ajusta modelos de seleção amostral utilizando os dois passos de Heckman corrige um pouco esse resultado, apesar de continuar assumindo distribuição normal bivariada para os erros do modelo o ajuste é mais sensível no sentido de detectar o viés de seleção amostral. Vejamos a análise com esse pacote.

Ajuste com o pacote SSMROB

#Equação de seleção que será ajustada via modelo probit
selectEq <- dambexp ~ age+female+educ+blhisp+totchr+ins
#Equação de regressão que será ajustada via modelo linear simples
outcomeEq <- lnambx ~ age+female+educ+blhisp+totchr+ins
fit5<-heckitrob(outcomeEq,selectEq,control=heckitrob.control(tcc=3.2,weights.x1="robCov"))
summary(fit5)

## ------------------------------------------------------------- 
## Robust 2-step Heckman / heckit M-estimation 
## Probit selection equation: 
##                  Estimate  Std.Error t-value  p-value    
## XS(Intercept) -0.74914476 0.19506999  -3.840 1.23e-04 ***
## XSage          0.10541500 0.02769588   3.806 1.41e-04 ***
## XSfemale       0.68740832 0.06225762  11.040 2.41e-28 ***
## XSeduc         0.07011568 0.01146521   6.116 9.62e-10 ***
## XSblhisp      -0.39774532 0.06264878  -6.349 2.17e-10 ***
## XStotchr       0.83283613 0.08027772  10.370 3.24e-25 ***
## XSins          0.18256005 0.06371471   2.865 4.17e-03  **
## Outcome equation: 
##                  Estimate  Std.Error t-value  p-value    
## XO(Intercept)  5.40154264 0.27672891  19.520 7.53e-85 ***
## XOage          0.20061658 0.02450765   8.186 2.70e-16 ***
## XOfemale       0.25501033 0.06992954   3.647 2.66e-04 ***
## XOeduc         0.01324867 0.01161609   1.141 2.54e-01    
## XOblhisp      -0.15508435 0.06506654  -2.383 1.72e-02   *
## XOtotchr       0.48115830 0.03822948  12.590 2.52e-36 ***
## XOins         -0.06706633 0.05159205  -1.300 1.94e-01    
## imrData$IMR1  -0.67676033 0.25927579  -2.610 9.05e-03  **
## ---
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
## -------------------------------------------------------------

Se considerarmos os erros padrão, vemos que o erro padrão da estimativa robusta $(0,2593)$ é menor do que a do estimador clássico $(0,2907).$ O valor de p do teste de viés de seleção robusto é $p=0.009$ , o que leva à conclusão da presença de viés de seleção amostral mesmo o estimador clássico sugerindo a ausência de viés.

Os mesmos dados foram utilizados agora acrescentando a covariável rendimento na equação de seleção, impondo a restrição de exclusão do modelo, embora o uso da renda para esta finalidade é discutível. Todos os fatores considerados são fortes preditores da decisão de gastar. Notemos que os resultados quanto a significância do parâmetro que acompanha a covariável $IMR$ foram os mesmos.

Meps 2001 Com a Covariável Rendimento

Ajuste com as funções glm e lm do R

#Modelo probit
fit1<-glm( dambexp ~ age+female+educ+blhisp+totchr+ins+income, family = binomial(link = "probit"),data=MEPS2001)
summary(fit1)

## 
## Call:
## glm(formula = dambexp ~ age + female + educ + blhisp + totchr + 
##     ins + income, family = binomial(link = "probit"), data = MEPS2001)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4420   0.1278   0.3859   0.6254   1.5329  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.668644   0.192364  -3.476 0.000509 ***
## age          0.086815   0.027495   3.157 0.001591 ** 
## female       0.663505   0.060952  10.886  < 2e-16 ***
## educ         0.061884   0.011891   5.204 1.95e-07 ***
## blhisp      -0.365784   0.061731  -5.925 3.11e-09 ***
## totchr       0.795747   0.071472  11.134  < 2e-16 ***
## ins          0.169107   0.062594   2.702 0.006900 ** 
## income       0.002677   0.001313   2.038 0.041514 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2904.9  on 3327  degrees of freedom
## Residual deviance: 2391.0  on 3320  degrees of freedom
## AIC: 2407
## 
## Number of Fisher Scoring iterations: 6

MEPS2001$IMR <- invMillsRatio(fit1)$IMR1#Criar uma nova coluna com a covariável Razão Inversa de Mills

#Modelo Linear Simples com razão de Mill's
fit2 <- lm( lnambx ~ age+female+educ+blhisp+totchr+ins + IMR,
                data = MEPS2001[ MEPS2001$dambexp == 1, ] )
summary(fit2)

## 
## Call:
## lm(formula = lnambx ~ age + female + educ + blhisp + totchr + 
##     ins + IMR, data = MEPS2001[MEPS2001$dambexp == 1, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1002 -0.7869  0.0810  0.8656  4.6380 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.28893    0.28541  18.531  < 2e-16 ***
## age          0.20247    0.02395   8.455  < 2e-16 ***
## female       0.29213    0.07174   4.072 4.79e-05 ***
## educ         0.01239    0.01144   1.083  0.27873    
## blhisp      -0.18287    0.06465  -2.829  0.00471 ** 
## totchr       0.50063    0.04797  10.436  < 2e-16 ***
## ins         -0.04651    0.05235  -0.888  0.37440    
## IMR         -0.46371    0.28065  -1.652  0.09859 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.269 on 2794 degrees of freedom
## Multiple R-squared:  0.1926, Adjusted R-squared:  0.1906 
## F-statistic: 95.23 on 7 and 2794 DF,  p-value: < 2.2e-16

#wald.test(b = coef(fit3), Sigma = vcov(fit3), Terms = 8)

Ajuste com o pacote Sample Selection

Processo de duas etapas

library("sampleSelection")
# Two-step estimation
fit3<-heckit( dambexp ~ age+female+educ+blhisp+totchr+ins+income,
               lnambx ~ age+female+educ+blhisp+totchr+ins, MEPS2001 )
summary(fit3)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 3328 observations (526 censored and 2802 observed)
## 18 free parameters (df = 3311)
## Probit selection equation:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.668647   0.194125  -3.444 0.000579 ***
## age          0.086815   0.027456   3.162 0.001581 ** 
## female       0.663505   0.060965  10.883  < 2e-16 ***
## educ         0.061884   0.012039   5.140 2.90e-07 ***
## blhisp      -0.365784   0.061909  -5.908 3.81e-09 ***
## totchr       0.795750   0.071217  11.174  < 2e-16 ***
## ins          0.169107   0.062930   2.687 0.007240 ** 
## income       0.002677   0.001310   2.043 0.041131 *  
## Outcome equation:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.28893    0.28852  18.331  < 2e-16 ***
## age          0.20247    0.02422   8.359  < 2e-16 ***
## female       0.29213    0.07258   4.025 5.82e-05 ***
## educ         0.01239    0.01157   1.071  0.28427    
## blhisp      -0.18287    0.06534  -2.798  0.00516 ** 
## totchr       0.50063    0.04855  10.311  < 2e-16 ***
## ins         -0.04651    0.05297  -0.878  0.38002    
## Multiple R-Squared:0.1926,   Adjusted R-Squared:0.1906
##    Error terms:
##               Estimate Std. Error t value Pr(>|t|)
## invMillsRatio  -0.4637     0.2826  -1.641    0.101
## sigma           1.2914         NA      NA       NA
## rho            -0.3591         NA      NA       NA
## --------------------------------------------

Método de Máxima Verossimilhança

# ML estimation
fit4<-selection( dambexp ~ age+female+educ+blhisp+totchr+ins+income,
                    lnambx ~ age+female+educ+blhisp+totchr+ins, MEPS2001)
summary(fit4)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 4 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: -5836.219 
## 3328 observations (526 censored and 2802 observed)
## 17 free parameters (df = 3311)
## Probit selection equation:
##              Estimate Std. error t value  Pr(> t)    
## (Intercept) -0.676054   0.194029  -3.484 0.000493 ***
## age          0.087936   0.027421   3.207 0.001342 ** 
## female       0.662665   0.060938  10.874  < 2e-16 ***
## educ         0.061948   0.012029   5.150 2.61e-07 ***
## blhisp      -0.363938   0.061873  -5.882 4.05e-09 ***
## totchr       0.796951   0.071131  11.204  < 2e-16 ***
## ins          0.170137   0.062871   2.706 0.006807 ** 
## income       0.002708   0.001317   2.056 0.039746 *  
## Outcome equation:
##             Estimate Std. error t value  Pr(> t)    
## (Intercept)  5.04406    0.22813  22.111  < 2e-16 ***
## age          0.21197    0.02301   9.213  < 2e-16 ***
## female       0.34814    0.06011   5.791 6.98e-09 ***
## educ         0.01872    0.01055   1.774 0.075986 .  
## blhisp      -0.21857    0.05967  -3.663 0.000249 ***
## totchr       0.53992    0.03933  13.727  < 2e-16 ***
## ins         -0.02999    0.05109  -0.587 0.557221    
##    Error terms:
##       Estimate Std. error t value Pr(> t)    
## sigma  1.27102    0.01838  69.157  <2e-16 ***
## rho   -0.13060    0.14708  -0.888   0.375    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Neste caso a estimativa do parâmetro da razão inversa de Mills é:

IMR<-coef(summary(fit4))["rho",1]*coef(summary(fit4))["sigma",1]
IMR

## [1] -0.1659965

#res<-fit4$estimate
#IMR<-res[16]*res[17]
#IMR

Ajuste com o pacote SSMROB

#Equação de seleção que será ajustada via modelo probit
selectEq <- dambexp ~ age+female+educ+blhisp+totchr+ins+income
#Equação de regressão que será ajustada via modelo linear simples
outcomeEq <- lnambx ~ age+female+educ+blhisp+totchr+ins
fit5<-heckitrob(outcomeEq,selectEq,control=heckitrob.control(tcc=3.2,weights.x1="robCov"))
summary(fit5)

## ------------------------------------------------------------- 
## Robust 2-step Heckman / heckit M-estimation 
## Probit selection equation: 
##                   Estimate   Std.Error t-value  p-value    
## XS(Intercept) -0.700434110 0.196402980  -3.566 3.62e-04 ***
## XSage          0.094588603 0.028149130   3.360 7.79e-04 ***
## XSfemale       0.703608275 0.062981278  11.170 5.61e-29 ***
## XSeduc         0.062307763 0.012118799   5.141 2.73e-07 ***
## XSblhisp      -0.388618039 0.062800450  -6.188 6.09e-10 ***
## XStotchr       0.834052967 0.080225712  10.400 2.58e-25 ***
## XSins          0.172551353 0.064031784   2.695 7.04e-03  **
## XSincome       0.002534592 0.001343944   1.886 5.93e-02   .
## Outcome equation: 
##                  Estimate  Std.Error t-value  p-value    
## XO(Intercept)  5.40932908 0.27290745  19.820 1.96e-87 ***
## XOage          0.20028861 0.02447160   8.185 2.73e-16 ***
## XOfemale       0.25213915 0.06994945   3.605 3.13e-04 ***
## XOeduc         0.01318542 0.01157841   1.139 2.55e-01    
## XOblhisp      -0.15341996 0.06514440  -2.355 1.85e-02   *
## XOtotchr       0.47956288 0.03805277  12.600 2.04e-36 ***
## XOins         -0.06825757 0.05173658  -1.319 1.87e-01    
## imrData$IMR1  -0.68995326 0.25543589  -2.701 6.91e-03  **
## ---
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
## -------------------------------------------------------------

Womenwk - Conjunto de dados do Stata

library(haven)
dados <- read_dta("~/GitHub/Modelos de Heckman/Modelo-de-Heckman/Conjunto_dados/womenwk.dta")
attach(dados)
head(dados)

## # A tibble: 6 × 12
##           c1          c2          u           v county   age education
##        <dbl>       <dbl>      <dbl>       <dbl>  <dbl> <int>     <int>
## 1 -0.4362051 -0.09691816 -0.2181026 -0.37572718      1    22        10
## 2  0.3521407  0.30047631  0.1760704  0.46123438      2    36        10
## 3  1.0772466 -1.59596262  0.5386233 -0.37624397      3    28        10
## 4  1.0212833 -1.71049791  0.5106417 -0.49699912      4    37        10
## 5 -0.4429599  0.30833973 -0.2214800 -0.09251083      5    39        10
## 6 -0.4403342  0.61322883 -0.2201671  0.12598438      6    33        10
## # ... with 5 more variables: married <int>, children <int>, select <dbl>,
## #   wagefull <dbl>, wage <dbl>

Considere o ajuste de um modelo de regressão linear simples considerando somente as pessoas empregadas para analisar o salário em função da educação e da idade.

#Regressão linear simples usando Mínimos Quadrados
fit1<-lm(wage~education+age)
summary(fit1)

## 
## Call:
## lm(formula = wage ~ education + age)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8219  -3.6000  -0.1582   3.7215  19.2203 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.08488    0.88962   6.840 1.20e-11 ***
## education    0.89658    0.04981  18.001  < 2e-16 ***
## age          0.14657    0.01871   7.833 9.66e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.452 on 1340 degrees of freedom
##   (657 observations deleted due to missingness)
## Multiple R-squared:  0.2535, Adjusted R-squared:  0.2524 
## F-statistic: 227.5 on 2 and 1340 DF,  p-value: < 2.2e-16

wage[is.na(wage)]<-0
wage2<-wage
indicadora<-ifelse(wage>0,1,0)
dados$wage2 <- wage2#Criar uma nova coluna com a covariável wage substituindo NA por zero 
dados$indicadora<-indicadora
str(dados)

## Classes 'tbl_df', 'tbl' and 'data.frame':    2000 obs. of  14 variables:
##  $ c1        : num  -0.436 0.352 1.077 1.021 -0.443 ...
##  $ c2        : num  -0.0969 0.3005 -1.596 -1.7105 0.3083 ...
##  $ u         : num  -0.218 0.176 0.539 0.511 -0.221 ...
##  $ v         : num  -0.3757 0.4612 -0.3762 -0.497 -0.0925 ...
##  $ county    : num  1 2 3 4 5 6 7 8 9 0 ...
##  $ age       : int  22 36 28 37 39 33 57 45 39 25 ...
##  $ education : int  10 10 10 10 10 10 10 16 12 10 ...
##  $ married   : int  1 1 1 1 1 1 1 1 1 0 ...
##  $ children  : int  0 0 0 0 1 2 1 0 0 3 ...
##  $ select    : num  16.8 32.4 19.2 21.3 32 ...
##  $ wagefull  : num  12.8 20.3 23.1 24.5 16.1 ...
##  $ wage      : num  NA 20.3 NA NA 16.1 ...
##  $ wage2     : num  0 20.3 0 0 16.1 ...
##  $ indicadora: num  0 1 0 0 1 1 1 1 0 1 ...

Agora o ajuste considerando todas as pessoas, inclusive as que não trabalham:

#Regressão linear simples usando Mínimos Quadrados
fit2<-lm(wage2~education+age)
summary(fit2)

## 
## Call:
## lm(formula = wage2 ~ education + age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.443  -9.809   3.021   8.596  29.882 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -12.16843    1.39815  -8.703   <2e-16 ***
## education     1.06457    0.08442  12.610   <2e-16 ***
## age           0.39077    0.03103  12.593   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.17 on 1997 degrees of freedom
## Multiple R-squared:  0.1726, Adjusted R-squared:  0.1718 
## F-statistic: 208.3 on 2 and 1997 DF,  p-value: < 2.2e-16

Vamos supor agora que o salário por hora é uma função da escolaridade e idade, ao passo que a probabilidade de trabalhar (a probabilidade de o salário ser observado) é uma função do estado civil, o número de crianças em casa, e (implicitamente) o salário ( que é estimado através da inclusão de idade e educação, que pensamos ser o que determinam o salário).

Heckman assume que o salário é a variável dependente e que a primeira lista de variáveis (educ e idade) são os determinantes do salário. As variáveis especificadas na opção de seleção (casado, crianças, educ, e idade) são assumidos para determinar se a variável dependente (a equação de regressão) é observada. Assim, ajustamos o modelo:

$\begin{eqnarray*} wage_{i} &=& \beta_{0} + \beta_{1}educ + \beta_{2}age + u_{1} \end{eqnarray*}$

e assumimos que o salário $(wage)$ é observado se:

$\begin{eqnarray*} \beta_{0}+\beta_{1}married+\beta_{2}children_{i}+\beta_{3}educ_{i}+\beta_{4}age+u_{2}>0, \end{eqnarray*}$

onde $u_{1}$ e $u_{2}$ são correlacionados.

Ajuste com as funções glm e lm do R

# Ajuste com as funções glm e lm do R
#Modelo probit
fit3<-glm(indicadora ~ married+children+education+age, family = binomial(link = "probit"),data=dados)
summary(fit3)

## 
## Call:
## glm(formula = indicadora ~ married + children + education + age, 
##     family = binomial(link = "probit"), data = dados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7594  -0.9414   0.4552   0.8459   2.0427  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.467365   0.192291 -12.831  < 2e-16 ***
## married      0.430857   0.074310   5.798 6.71e-09 ***
## children     0.447325   0.028642  15.618  < 2e-16 ***
## education    0.058365   0.011018   5.297 1.18e-07 ***
## age          0.034721   0.004232   8.204 2.33e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2532.4  on 1999  degrees of freedom
## Residual deviance: 2054.1  on 1995  degrees of freedom
## AIC: 2064.1
## 
## Number of Fisher Scoring iterations: 5

library(sampleSelection)
dados$IMR <- invMillsRatio(fit3)$IMR1#Criar uma nova coluna com a covariável Razão Inversa de Mills

#Modelo Linear Simples com razão de Mill's
fit4 <- lm(wage2 ~ education+age + IMR,
           data = dados[dados$wage2>0, ] )
summary(fit4)

## 
## Call:
## lm(formula = wage2 ~ education + age + IMR, data = dados[dados$wage2 > 
##     0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.4167  -3.5143  -0.1737   3.5506  20.3762 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.73404    1.16621   0.629    0.529    
## education    0.98253    0.05050  19.457  < 2e-16 ***
## age          0.21187    0.02066  10.253  < 2e-16 ***
## IMR          4.00162    0.57710   6.934 6.35e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.359 on 1339 degrees of freedom
## Multiple R-squared:  0.2793, Adjusted R-squared:  0.2777 
## F-statistic:   173 on 3 and 1339 DF,  p-value: < 2.2e-16

Ajuste com o pacote sampleSelection

Processo de duas etapas

#Equação de seleção que será ajustada via modelo probit
selectEq <- indicadora~married+children+education+age
#Equação de regressão que será ajustada via modelo linear simples
outcomeEq <- wage2 ~ education+age
fit5<-heckit(selectEq,outcomeEq,dados)
summary(fit5)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 2000 observations (657 censored and 1343 observed)
## 11 free parameters (df = 1990)
## Probit selection equation:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.467365   0.192563 -12.813  < 2e-16 ***
## married      0.430857   0.074208   5.806 7.43e-09 ***
## children     0.447325   0.028742  15.564  < 2e-16 ***
## education    0.058365   0.010974   5.318 1.16e-07 ***
## age          0.034721   0.004229   8.210 3.94e-16 ***
## Outcome equation:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.73404    1.24833   0.588    0.557    
## education    0.98253    0.05388  18.235   <2e-16 ***
## age          0.21187    0.02205   9.608   <2e-16 ***
## Multiple R-Squared:0.2793,   Adjusted R-Squared:0.2777
##    Error terms:
##               Estimate Std. Error t value Pr(>|t|)    
## invMillsRatio   4.0016     0.6065   6.597 5.35e-11 ***
## sigma           5.9474         NA      NA       NA    
## rho             0.6728         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Método de Máxima Verossimilhança

#Estimação utilizando método de Máxima Verossimilhança do pacote sample selection
fit6<-selection(indicadora ~ married+children+education+age,
                wage2 ~ education+age)
summary(fit6)

## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 3 iterations
## Return code 2: successive function values within tolerance limit
## Log-Likelihood: -5178.304 
## 2000 observations (657 censored and 1343 observed)
## 10 free parameters (df = 1990)
## Probit selection equation:
##              Estimate Std. error t value  Pr(> t)    
## (Intercept) -2.491015   0.189340 -13.156  < 2e-16 ***
## married      0.445171   0.067395   6.605 3.97e-11 ***
## children     0.438707   0.027783  15.791  < 2e-16 ***
## education    0.055732   0.010735   5.192 2.08e-07 ***
## age          0.036510   0.004153   8.790  < 2e-16 ***
## Outcome equation:
##             Estimate Std. error t value Pr(> t)    
## (Intercept)  0.48579    1.07704   0.451   0.652    
## education    0.98995    0.05326  18.588  <2e-16 ***
## age          0.21313    0.02060  10.345  <2e-16 ***
##    Error terms:
##       Estimate Std. error t value Pr(> t)    
## sigma  6.00479    0.16572   36.23  <2e-16 ***
## rho    0.70350    0.05123   13.73  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

A estimativa do parâmetro que acompanha a covariável razão inversa de Mills é:

IMR<-coef(summary(fit6))["rho",1]*coef(summary(fit6))["sigma",1]
IMR

## [1] 4.224401

Ajuste com o pacote SSMROB

library(ssmrob)
#Equação de seleção que será ajustada via modelo probit
selectEq <- indicadora ~ married+children+education+age
#Equação de regressão que será ajustada via modelo linear simples
outcomeEq <- wage2 ~ education+age
fit7<-heckitrob(outcomeEq,selectEq,control=heckitrob.control(tcc=3.2,weights.x1="robCov"))
summary(fit7)

## ------------------------------------------------------------- 
## Robust 2-step Heckman / heckit M-estimation 
## Probit selection equation: 
##                  Estimate   Std.Error t-value  p-value    
## XS(Intercept) -2.46541991 0.193678111 -12.730 4.06e-37 ***
## XSmarried      0.42133992 0.074544389   5.652 1.58e-08 ***
## XSchildren     0.44264144 0.029190486  15.160 6.13e-52 ***
## XSeducation    0.05884727 0.011046541   5.327 9.97e-08 ***
## XSage          0.03479865 0.004246128   8.195 2.50e-16 ***
## Outcome equation: 
##                Estimate  Std.Error t-value  p-value    
## XO(Intercept) 0.7458634 1.24334254  0.5999 5.49e-01    
## XOeducation   0.9798627 0.05471929 17.9100 1.04e-71 ***
## XOage         0.2111380 0.02197695  9.6070 7.45e-22 ***
## imrData$IMR1  4.0074119 0.60700579  6.6020 4.06e-11 ***
## ---
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
## -------------------------------------------------------------

Maiores informações podem ser obtidas em:

Heckman (1976), Heckman (1977), Mroz (1987), Winship and Mare (1992), Puhani (2000), Greene (2003), Cameron and Trivedi (2009) e Marchenko and Genton (2012).

Referências

Cameron, A Colin, and Pravin K Trivedi. 2009. “Microeconomics Using Stata.” College Station, TX: Stata Press Publications.

Greene, William H. 2003. Econometric Analysis. Pearson Education India.

Heckman, James J. 1976. “The Common Structure of Statistical Models of Truncation, Sample Selection and Limited Dependent Variables and a Simple Estimator for Such Models.” In Annals of Economic and Social Measurement, Volume 5, Number 4, 475–92. NBER.

———. 1977. “Sample Selection Bias as a Specification Error (with an Application to the Estimation of Labor Supply Functions).” National Bureau of Economic Research Cambridge, Mass., USA.

Hill, R Carter, William E Griffiths, and Guay C Lim. 2008. Principles of Econometrics. Vol. 5. Wiley Hoboken, NJ.

Marchenko, Yulia V, and Marc G Genton. 2012. “A Heckman Selection-T Model.” Journal of the American Statistical Association 107 (497). Taylor & Francis Group: 304–17.

Mroz, Thomas A. 1987. “The Sensitivity of an Empirical Model of Married Women’s Hours of Work to Economic and Statistical Assumptions.” Econometrica: Journal of the Econometric Society. JSTOR, 765–99.

Puhani, Patrick. 2000. “The Heckman Correction for Sample Selection and Its Critique.” Journal of Economic Surveys 14 (1). Wiley Online Library: 53–68.

Winship, Christopher, and Robert D Mare. 1992. “Models for Sample Selection Bias.” Annual Review of Sociology. JSTOR, 327–50.