Exemplos de Aplicação do Modelo de Seleção de Heckman
Fernando
2016-09-22
1 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 regressã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 1ddados[423:433,c(5,6,19,23,14,1,7,21)]## age educ exper kids mtr inlf wage lwage
## 423 32 17 14 1 0.7215 1 4.7826 1.5649840
## 424 36 10 2 1 0.7215 1 2.3118 0.8380265
## 425 40 12 21 1 0.6215 1 5.3061 1.6688570
## 426 43 13 22 1 0.5815 1 5.8675 1.7694290
## 427 33 12 14 1 0.5815 1 3.4091 1.2264480
## 428 30 12 7 1 0.6915 1 4.0816 1.4064890
## 429 49 12 2 1 0.6615 0 NA NA
## 430 30 16 5 1 0.6615 0 NA NA
## 431 30 12 12 1 0.6615 0 NA NA
## 432 41 12 1 1 0.5815 0 NA NA
## 433 45 12 12 1 0.6615 0 NA NAattach(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-15O 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.
1.1 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: 4Como 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)]## age educ exper kids mtr inlf wage lwage IMR
## 423 32 17 14 1 0.7215 1 4.7826 1.5649840 0.4412366
## 424 36 10 2 1 0.7215 1 2.3118 0.8380265 0.8181494
## 425 40 12 21 1 0.6215 1 5.3061 1.6688570 0.6793267
## 426 43 13 22 1 0.5815 1 5.8675 1.7694290 0.6340350
## 427 33 12 14 1 0.5815 1 3.4091 1.2264480 0.5657426
## 428 30 12 7 1 0.6915 1 4.0816 1.4064890 0.6164304
## 429 49 12 2 1 0.6615 0 NA NA 0.8290024
## 430 30 16 5 1 0.6615 0 NA NA 0.4219251
## 431 30 12 12 1 0.6615 0 NA NA 0.5929940
## 432 41 12 1 1 0.5815 0 NA NA 0.6586579
## 433 45 12 12 1 0.6615 0 NA NA 0.7763784Este é 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-16Notemos 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:
1.2 Ajuste com o pacote Sample Selection
1.2.1 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.
1.2.2 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]
#IMR1.3 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
## -------------------------------------------------------------2 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-16Podemos 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-16Mas 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.
2.1 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: 6Com 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.25274431Assim, 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 (???), 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:
2.2 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.\)
2.2.1 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.
2.2.2 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]
#IMRPoderí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.
2.3 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.
3 Meps 2001 Com a Covariável Rendimento
3.1 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: 6MEPS2001$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)3.2 Ajuste com o pacote Sample Selection
3.2.1 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
## --------------------------------------------3.2.2 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]
#IMR3.3 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
## -------------------------------------------------------------4 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)## c1 c2 u v county age education
## 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
## married children select wagefull wage
## 1 1 0 16.79127 12.78277 NA
## 2 1 0 32.43481 20.31285 20.31285
## 3 1 0 19.18507 23.06348 NA
## 4 1 0 21.33601 24.52770 NA
## 5 1 1 31.98987 16.14224 16.14224
## 6 1 2 37.21181 14.95799 14.95799Considere 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-16wage[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-16Vamos 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.
4.1 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: 5library(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-164.2 Ajuste com o pacote sampleSelection
4.2.1 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
## --------------------------------------------4.2.2 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.2244014.3 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
## -------------------------------------------------------------Referências
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.