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;
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 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}).\)
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}\]
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.
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:
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
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
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.
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:
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.
#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.
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
#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
#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
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
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.