Modelo de selección de Heckman

Fri, Sep 4, 2020 9-minute read

Recordemos que en el modelo Tobit existen dos categorías de observaciones: 1) aquellas cuyo valor de la variable dependiente es positivo y 2) aquellas cuyo valor de la variable dependiente está censurado y es, de hecho, igual a cero. Es decir:

y={y si y>00 si y0y=\begin{cases} y^* & \ \text{si} \ y^*>0 \\ 0 & \ \text{si} \ y^*\leq 0\\ \end{cases}

En el modelo, el conjunto de variables (x\mathbf{x}) que intentan explicar el valor de la variable dependiente cuando esta es positiva y las razones por las cuales esta cae dentro de una categoría u otra (y=0y=0 o y>0y>0) son las mismas. Esto constituye una de las principales limitaciones del modelo Tobit. De hecho, esta es una de las principales razones por las cuales el modelo Tobit se utiliza raramente en la práctica.

Para generalizar el modelo Tobit, en el modelo de selección de Heckman consideramos dos variables latentes en lugar de una. Sean estas ww^* y yy^*. La variable ww^* define el proceso de censura. Esto es:

wi={1 si wi>00 si wi0.w_i=\begin{cases} 1 & \ \text{si} \ w_i^*>0 \\ 0 & \ \text{si} \ w_i^*\leq 0\\ \end{cases}.

La variable yy^* es equivalente a la variable latente en el modelo Tobit. Esto es:

yi={yi si wi>00 si wi0.y_i=\begin{cases} y_i^* & \ \text{si} \ w_i^*>0 \\ 0 & \ \text{si} \ w_i^*\leq 0\\ \end{cases}.

Finalmente:

wi=zγ+εyi=xβ+u,\begin{aligned} w_i^*&= \mathbf{z}\boldsymbol{\gamma}+\boldsymbol{\varepsilon}\\ y_i^*&= \mathbf{x}\boldsymbol{\beta}+\mathbf{u}\\ \end{aligned},

donde zx\mathbf{z}\supseteq\mathbf{x}, es decir, z\mathbf{z} puede incluir algunas variables no incluidas en x\mathbf{x} y γ\boldsymbol{\gamma} y β\boldsymbol{\beta} son vectores de coeficientes a estimar. Si la distribución conjunta de ε\boldsymbol{\varepsilon} y u\mathbf{u} es normal bivariada con parámetro de correlación ρ\rho, tenemos:

[εu]N(0,Σ); Σ=[1ρσρσ1]\begin{bmatrix}\boldsymbol{\varepsilon}\\\mathbf{u}\end{bmatrix}\sim\mathcal{N}(0,\boldsymbol{\Sigma}); \ \boldsymbol{\Sigma}=\begin{bmatrix} 1 & \rho\sigma \\ \rho\sigma & 1 \end{bmatrix} Aquí: E(yix)=xβ\mathbb{E}(y_i^*|\mathbf{x})=\mathbf{x}\boldsymbol{\beta} y, para la variable observada yiy_i:

E(yiyi>0,x)=xβ+ρσλ(zγ)E(yix)=Φ(zγ)×(xβ+ρσλ(zγ)),\begin{aligned} \mathbb{E}(y_i|y_i>0, \mathbf{x}) &= \mathbf{x}\boldsymbol{\beta} + \rho\sigma\lambda(\mathbf{z}\boldsymbol{\gamma})\\ \mathbb{E}(y_i|\mathbf{x}) &= \Phi(\mathbf{z}\boldsymbol{\gamma}) \times \left(\mathbf{x}\boldsymbol{\beta}+\rho\sigma\lambda(\mathbf{z}\boldsymbol{\gamma})\right), \end{aligned}

donde λ(zγ)=ϕ(zγ)Φ(zγ)\lambda(\mathbf{z}\boldsymbol{\gamma})=\frac{\phi(\mathbf{z}\boldsymbol{\gamma})}{\Phi(\mathbf{z}\boldsymbol{\gamma})} es la inversa del cociente de Mills, también conocido como razón de riesgo.

Estimación del modelo

Existen dos formas de estimar los coeficientes del modelo de selección. La primera es mediante máxima verosimilitud con información completa. La función de verosimilitud en este caso tiene tres elementos. El primero es la probabilidad de observar a la variable dependiente (es decir, la probabilidad de observar un yi>0y_i>0); el segundo es la probabilidad de no observar a la variable dependiente (es decir, la probabilidad de observar un yi=0y_i=0); finalmente, el tercer elemento es la distribución que sigue yiy_i cuando yi>0y_i>0. Suponiendo que el término de error sigue una distribución normal, tenemos que la función de verosimilitud es:

L=i=1N{Φ(zγ)}1wi{Φ(zγ+[yixβ]ρσ1ρ2)}wi{1σ2πe12(yixβσ)2}wiL = \prod_{i=1}^{N}\left\{\Phi(-\mathbf{z}\boldsymbol{\gamma})\right\}^{1-w_i}\left\{\Phi\left(\frac{\mathbf{z}\boldsymbol{\gamma}+\left[y_i-\mathbf{x}\boldsymbol{\beta}\right]\frac{\rho}{\sigma}}{\sqrt{1-\rho^2}}\right)\right\}^{w_i}\left\{\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{y_i-\mathbf{x}\boldsymbol{\beta}}{\sigma}\right)^2}\right\}^{w_i} La segunda forma de estimar es computacionalmente más simple. La idea es tratar el problema de selección como si se tratase de un problema de variables omitidas (Heckman, 1979). En este caso maximizamos una función de verosimilitud con información limitada . El procedimiento se realiza en dos etapas. Primero se utiliza un modelo Probit para estimar la probabilidad de yi=0y_i=0 y de yi>0y_i>0. Con estas, se calcula la inversa del cociente de Mills, λ^=ϕ(zγ^)Φ(zγ^)\widehat{\lambda}=\frac{\phi(\mathbf{z}\widehat{\boldsymbol{\gamma}})}{\Phi(\mathbf{z}\widehat{\boldsymbol{\gamma}})}. En la segunda etapa se estima el modelo original por el método de los MCO, pero se agrega λ^\widehat{\lambda} como regresor (la variable “omitida”):

E(yiyi>0,x)=xβ+ρσλ^\mathbb{E}(y_i|y_i>0, \mathbf{x}) = \mathbf{x}\boldsymbol{\beta}+\rho\sigma\widehat{\boldsymbol{\lambda}} Naturalmente, cuando ρ=0\rho=0, la inversa del cociente de Mills desaparece y, por lo tanto, el modelo de selección se hace innecesario.

Ejemplo

Para estimar los coeficientes en un modelo de selección de Heckman utilizando R, necesitamos los siguientes paquetes:

# Instalar paquetes si es necesario
# install.packages("sampleSelection")

# Leer paquetes
library(sampleSelection)
## Warning: package 'sampleSelection' was built under R version 4.0.2
## Loading required package: maxLik
## Warning: package 'maxLik' was built under R version 4.0.2
## Loading required package: miscTools
## Warning: package 'miscTools' was built under R version 4.0.2
## 
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
## 
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/

La base de datos Mroz87 está incluida en el paquete `sampleSelection e incluye los datos laborales de una muestra de mujeres casadas. El año del levantamiento de la información es 1975.

data(Mroz87)

Estimamos una ecuación para el logaritmo del salario. Este es función de la escolaridad (medida en años), la experiencia (medida en años), el cuadrado de la experiencia y una variable dummy que vale 1 si la observación i-ésima vivía en una ciudad considerada grande.

El problema la estimación de los coeficientes en el modelo de Mroz es que solo observamos el salario de las mujeres que, en ese momento, estaban trabajando en el mercado formal. Una gran proporción de mujeres casadas en 1975 no formaban parte de la fuerza laboral (o, si lo estaban, no eran trabajadores formales). Entonces, en la muestra tenemos una gran cantidad de ceros. En la base de datos, la variable lfp es una variable dummy que vale 1 si la observación i-ésima participaba en el mercado formal de trabajo.

Estimamos primero el modelo por el método de los MCO. Nota que restringimos la muestra a las mujeres participantes en el mercado laboral:

ols <- lm(log(wage) ~ educ + exper + I( exper^2 ) + city, data=subset(Mroz87, lfp==1))
summary(ols)
## 
## Call:
## lm(formula = log(wage) ~ educ + exper + I(exper^2) + city, data = subset(Mroz87, 
##     lfp == 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.10084 -0.32453  0.05292  0.36261  2.34806 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.5308476  0.1990253  -2.667  0.00794 ** 
## educ         0.1057097  0.0143280   7.378 8.58e-13 ***
## exper        0.0410584  0.0131963   3.111  0.00199 ** 
## I(exper^2)  -0.0007973  0.0003938  -2.025  0.04352 *  
## city         0.0542225  0.0680903   0.796  0.42629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6667 on 423 degrees of freedom
## Multiple R-squared:  0.1581, Adjusted R-squared:  0.1501 
## F-statistic: 19.86 on 4 and 423 DF,  p-value: 5.389e-15

Ahora, estimamos el modelo utilizando el procedimiento de dos etapas de Heckman. Para esto utilizamos la función heckit() del paquete sampleSelection. Nota que la sintaxis de esta función es muy similar a la de la estimación por MCO, pero que ahora tenemos dos ecuaciones en lugar de una. En la primer ecuación establecemos los determinantes de lfm, que es la variable que indica si una mujer está o no en el mercado laboral. Esta es una variable binaria. (En un modelo de selección de Heckman siempre tiene que serlo.) La segunda ecuación es la del (logaritmo) del salario.

Mroz87$kids = ( Mroz87$kids5 + Mroz87$kids618 > 0 )
heckman = heckit( lfp ~ age + I( age^2 ) + kids + huswage + educ,
                 log(wage) ~ educ + exper + I( exper^2 ) + city, data=Mroz87 )
summary(heckman)
## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 753 observations (325 censored and 428 observed)
## 14 free parameters (df = 740)
## Probit selection equation:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.6067983  1.4072534  -3.274 0.001111 ** 
## age          0.2043835  0.0663822   3.079 0.002155 ** 
## I(age^2)    -0.0026174  0.0007791  -3.359 0.000821 ***
## kidsTRUE    -0.4270345  0.1309974  -3.260 0.001166 ** 
## huswage     -0.0419745  0.0121387  -3.458 0.000576 ***
## educ         0.1313128  0.0227687   5.767 1.18e-08 ***
## Outcome equation:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.7524012  0.3926416  -1.916  0.05572 .  
## educ         0.1152050  0.0204064   5.646 2.35e-08 ***
## exper        0.0427268  0.0133502   3.200  0.00143 ** 
## I(exper^2)  -0.0008469  0.0003984  -2.126  0.03383 *  
## city         0.0446548  0.0692400   0.645  0.51917    
## Multiple R-Squared:0.1589,   Adjusted R-Squared:0.149
##    Error terms:
##               Estimate Std. Error t value Pr(>|t|)
## invMillsRatio   0.1501     0.2293   0.655    0.513
## sigma           0.6720         NA      NA       NA
## rho             0.2234         NA      NA       NA
## --------------------------------------------

Nota que la primera ecuación en el procedimiento de dos etapas de Heckman hemos incluido a todas las variables que podrían afectar el proceso de selección. En nuestro ejemplo tenemos a las variables que podrían explicar las razones por las cuales una mujer se encuentra o no en el mercado laboral. Es recomendable que en esta ecuación tengamos, al menos, una variable que actúe como instrumento, es decir que afecte la decisión de estar o no en el mercado laboral pero que no explique directamente al salario. Esta es la función de la variable kids, la cual es una dummy que vale uno si la observación i-ésima tiene hijos menores a 16 años. Una mujer con hijos podría decidir no trabajar, pero, una vez que trabaja, el tener hijos no tiene por qué afectar la paga (que, en este caso, se mide por horas).

En cuanto a la segunda ecuación, la interpretación de los coeficientes es idéntica a la del modelo estimado por MCO. Por ejemplo, cada año de escolaridad adicional aumenta el salario en 11.5% (recuerda que la variable dependiente está medida en logaritmos).

El paquete sampleSelection también nos permite estimar el modelo por máxima verosimilitud con información completa. Esto se hace utilizando la función selection. Esta estima ambas ecuaciones del modelo de selección simultáneamente (y no secuencialmente como en el caso del procedimiento en dos etapas de Heckman):

mselection <- selection(lfp ~ age + I( age^2 ) + kids + huswage + educ, 
                        log(wage) ~ educ + exper + I( exper^2 ) + city, data=Mroz87)
summary(mselection)
## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 4 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: -916.2869 
## 753 observations (325 censored and 428 observed)
## 13 free parameters (df = 740)
## Probit selection equation:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.5123815  1.4162713  -3.186 0.001503 ** 
## age          0.1992498  0.0670579   2.971 0.003061 ** 
## I(age^2)    -0.0025569  0.0007875  -3.247 0.001220 ** 
## kidsTRUE    -0.4238143  0.1312568  -3.229 0.001298 ** 
## huswage     -0.0428919  0.0121978  -3.516 0.000464 ***
## educ         0.1325479  0.0229000   5.788 1.05e-08 ***
## Outcome equation:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.6585505  0.2962345  -2.223  0.02651 *  
## educ         0.1111762  0.0170956   6.503 1.45e-10 ***
## exper        0.0420022  0.0132165   3.178  0.00154 ** 
## I(exper^2)  -0.0008252  0.0003943  -2.093  0.03669 *  
## city         0.0486690  0.0683836   0.712  0.47687    
##    Error terms:
##       Estimate Std. Error t value Pr(>|t|)    
## sigma  0.66590    0.02533  26.288   <2e-16 ***
## rho    0.13033    0.22289   0.585    0.559    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------
# Nota. La función `selection()` también nos permite estimar el modelo en dos etapas usando el procedimiento de Heckman. Basta con agregar a la función el argumento `2step`:

# mselectionh <- selection(lfp ~ age + I( age^2 ) + kids + huswage + educ, 
#                         log(wage) ~ educ + exper + I( exper^2 ) + city, data=Mroz87, method="2step")
# summary(mselectionh)

Como puede constatarse, la estimación del modelo de selección no tiene grandes efectos sobre los parámetros de interés.Los coeficientes estimados en todas las ecuaciones no son entre sí muy diferentes. Esto podría implicar que la estimación por MCO podría resultar, después de todo, en estimadores no tan sesgados. De hecho, esto resulta ser el caso, puesto que el valor de ρ\rho en la estimación del modelo de selección no es estadísticamente diferente de cero (valor p de 0.559).

Última actualización: 07-10-2020.