Monitoria 4 de Econometria II - Variável Dependente Limitada

Regressão com Variável Dependente Binária

Neste capítulo, discutimos uma classe especial de modelos de regressão que visam explicar uma variável dependente limitada. Em particular, consideramos modelos onde a variável dependente é binária. Veremos que em tais modelos, a função de regressão pode ser interpretada como uma função de probabilidade condicional da variável dependente binária.

Revisamos os seguintes conceitos:

o modelo de probabilidade linear
modelo Probit
modelo Logit
estimativa de máxima verossimilhança de modelos de regressão não linear

É claro que também veremos como estimar os modelos acima usando R e discutiremos uma aplicação onde examinaremos a questão de saber se existe discriminação racial no mercado hipotecário dos EUA.

Modelo de Probabilidade Linear

1.1 Modelo de Probabilidade Linear

O modelo de regressão linear

$Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \dots + \beta_k X_{ki} + u_i$

com uma variável dependente binária $Y_i$ é chamado de modelo de probabilidade linear. No modelo de probabilidade linear temos

$E(Y\vert X_1,X_2,\dots,X_k) = P(Y=1\vert X_1, X_2,\dots, X_3)$

onde,

$PP(Y = 1 \vert X_1, X_2, \dots, X_k) = \beta_0 + \beta_1 + X_{1i} + \beta_2 X_{2i} + \dots + \beta_k X_{ki}.$

Assim, $\beta_j$ pode ser interpretado como a mudança na probabilidade tal que $Y_i =1$, mantendo constantes os outros $k-1$ regressores constantes. Assim como na regressão múltipla comum, $\beta_j$pode ser estimado usando OLS e as fórmulas robustas de erro padrão podem ser usadas para testes de hipóteses e cálculo de intervalos de confiança.

Na maioria dos modelos de probabilidade linear, $R^2$ não tem interpretação significativa, uma veque a linha de regressão nunca pode ajustar os dados perfeitamente se a variável dependente for binária e os regressores forem contínuos. Isso pode ser visto no aplicativo abaixo.

É essencial usar erros padrão robustos, uma vez que em um modelo de probabilidade linear $u_i$ são sempre heterocedásticos.

Modelos de probabilidade linear são facilmente estimados em R usando a função $lm()$ .

Aplicação: Dados Hipotecários

Seguindo o livro, começamos por carregar o conjunto de dados HDMA que fornece dados relativos a pedidos de hipoteca apresentados em Boston no ano de 1990.

# load `AER` package and attach the HMDA data
library(AER)
library(tidyverse)
library(stargazer)
data(HMDA)

summary(HMDA)

##   deny          pirat            hirat            lvrat        chist   
##  no :2095   Min.   :0.0000   Min.   :0.0000   Min.   :0.0200   1:1353  
##  yes: 285   1st Qu.:0.2800   1st Qu.:0.2140   1st Qu.:0.6527   2: 441  
##             Median :0.3300   Median :0.2600   Median :0.7795   3: 126  
##             Mean   :0.3308   Mean   :0.2553   Mean   :0.7378   4:  77  
##             3rd Qu.:0.3700   3rd Qu.:0.2988   3rd Qu.:0.8685   5: 182  
##             Max.   :3.0000   Max.   :3.0000   Max.   :1.9500   6: 201  
##  mhist    phist          unemp        selfemp    insurance  condomin  
##  1: 747   no :2205   Min.   : 1.800   no :2103   no :2332   no :1694  
##  2:1571   yes: 175   1st Qu.: 3.100   yes: 277   yes:  48   yes: 686  
##  3:  41              Median : 3.200                                   
##  4:  21              Mean   : 3.774                                   
##                      3rd Qu.: 3.900                                   
##                      Max.   :10.600                                   
##   afam      single     hschool   
##  no :2041   no :1444   no :  39  
##  yes: 339   yes: 936   yes:2341  
##                                  
##                                  
##                                  
##

A variável que estamos interessados em modelar é deny, um indicador que indica se o pedido de hipoteca de um requerente foi aceito (deny = no) ou negado (deny = yes). Um regressor que deveria ter poder para explicar se um pedido de hipoteca foi negado é o pirar, o tamanho do total antecipado de pagamentos mensais do empréstimo em relação à renda do requerente. É simples traduzir isso no modelo de regressão simples

\[\begin{align} deny = \beta_0 + \beta_1 \times P/I \ ratio + u. \tag{1.1} \end{align}\]

Estimamos este modelo como qualquer outro modelo de regressão linear usando lm(). Antes de fazermos isso, a variável deny deve ser convertida em uma variável numérica usando as.numeric(), dado que lm() não aceita variável dependente da classe factor.

HMDA = HMDA %>% 
  mutate(deny = ifelse(deny == "yes", 1, 0))

table(HMDA$deny)

## 
##    0    1 
## 2095  285

# estimate a simple linear probabilty model
denymod1 <- lm(deny ~ pirat, data = HMDA)
denymod1

## 
## Call:
## lm(formula = deny ~ pirat, data = HMDA)
## 
## Coefficients:
## (Intercept)        pirat  
##    -0.07991      0.60353

A seguir, plotamos os dados e a linha de regressão

# plot the data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Scatterplot Mortgage Application Denial and 
                                    the Payment-to-Income Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.8)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add the estimated regression line
abline(denymod1, 
       lwd = 1.8, 
       col = "steelblue")

De acordo com o modelo estimado, a razão entre pagamentos e rendimentos de 1 está associado a uma probabilidade esperada de negação do pedido de hipoteca de aproximadamente $50%$ . O modelo indica que existe uma relação positiva entre o razão pagamento/rendimento e a probabilidade de um pedido de hipoteca ser negado, pelo que os indivíduos com um elevado razão entre pagamentos de empréstimos e rendimento têm maior probabilidade de serem rejeitados.

Podemos usar coeftest() para obter erros padrão robustos para ambas as estimativas de coeficientes.

# print robust coefficient summary
coeftest(denymod1, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.079910   0.031967 -2.4998   0.01249 *  
## pirat        0.603535   0.098483  6.1283 1.036e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A linha de regressão estimada é

\[ \begin{align} \widehat{deny} = -\underset{(0,032)}{0,080} + \underset{(0,098)}{0,604} P/I \ ratio. \tag{11.2} \end{align} \]

O verdadeiro coeficiente na razão $P/I$ é estatisticamente diferente de 0 no nível $1\%$. Sua estimativa pode ser interpretada da seguinte forma: um aumento de 1 ponto percentual no índice $P/I$ leva a um aumento na probabilidade de recusa de um empréstimo em $0,604 \cdot 0,01 = 0,00604 \approx 0,6\%$.

Se aumentarmos o modelo simples (1.1) com um regressor de cor preta adicional que é igual a $1$ se o requerente for um afro-americano e igual a $0$ caso contrário. Tal especificação é a base para investigar se existe discriminação racial no mercado hipotecário: se ser negro tem uma influência significativa (positiva) na probabilidade de recusa de um empréstimo quando controlamos factores que permitem uma avaliação objectiva dea qualidade de crédito de um requerente, este é um indicador de discriminação.

HMDA = HMDA %>% 
  rename(black = afam)

# estimate the model
denymod2 <- lm(deny ~ pirat + black, data = HMDA)
coeftest(denymod2, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.090514   0.033430 -2.7076  0.006826 ** 
## pirat        0.559195   0.103671  5.3939 7.575e-08 ***
## blackyes     0.177428   0.025055  7.0815 1.871e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A função de regressão estimada é

\[\begin{align} \widehat{deny} =& \, -\underset{(0,029)}{0,091} + \underset{(0,089)}{0,559} P/I \ ratio + \underset{(0,025)}{ 0,177} black \tag{1.3} \end{align}\]

O coeficiente em $black$ é positivo e significativamente diferente de zero no nível de $0,01\%$ . A interpretação é que, mantendo constante a razão P/I, ser preto aumenta a probabilidade de um pedido de hipoteca ser negado em cerca de $17,7\%$ . Esta constatação é compatível com a discriminação racial. No entanto, pode ser distorcido pela omissão de preconceitos variáveis, pelo que a discriminação pode ser uma conclusão prematura.

Regressão Probit e Logit

O modelo de probabilidade linear tem uma grande falha: ele assume que a função de probabilidade condicional é linear. Isso não restringe $P(Y=1\vert X_1,\dots,X_k)$ a ficar entre $0$ e $1$. Podemos ver isso facilmente em nossa reprodução feita anteriormente: para $P/I \ ratio \geq 1.75$, (1.2) prevê que a probabilidade de negação de um pedido de hipoteca seja maior que $1$. Para aplicações com Razão $P/I$ próxima de $0$, a probabilidade prevista de negação é ainda negativa, de modo que o modelo não tem interpretação significativa aqui.

Esta circunstância exige uma abordagem que utilize uma função não linear para modelar a função de probabilidade condicional de uma variável dependente binária. Os métodos comumente usados são regressão Probit e Logit.

Probit

Na regressão Probit, a função de distribuição normal padrão cumulativa $\Phi(\cdot)$ é usada para modelar a função de regressão quando a variável dependente é binária, ou seja, assumimos

\[ \begin{align} E(Y\vert X) = P(Y=1\vert X) = \Phi(\beta_0 + \beta_1 X). \tag{11.4} \end{align} \]

$\beta_0 + \beta_1 X$ desempenha o papel de um quantil $z$.

Lembre-se que

\[ \Phi(z) = P(Z \leq z) \ , \ Z \sim \mathcal{N}(0,1) \]

tal que o coeficiente Probit $\beta_1$ em (11.4) é a mudança em $z$ associada a uma mudança de uma unidade em $X$ . Embora o efeito sobre $z$ de uma mudança em $X$ seja linear, a ligação entre $z$ e a variável dependente $Y$ é não linear, uma vez que $\Phi$ é uma função não linear de $X$.

Como a variável dependente é uma função não linear dos regressores, o coeficiente em $X$ não tem interpretação simples. De acordo com o Conceito Chave 8.1, a mudança esperada na probabilidade de $Y = 1$ devido a uma mudança na razão P/I pode ser calculada da seguinte forma:

Calcule a probabilidade prevista de que $Y=1$ para o valor original de $X$ .
Calcule a probabilidade prevista de que $Y=1$ para $X + \Delta X$x’.
Calcule a diferença entre as duas probabilidades previstas.

É claro que podemos generalizar (11.4) para a regressão Probit com regressores múltiplos para mitigar o risco de enfrentar vieses de variáveisomitidas. Os fundamentos da regressão Probit estão resumidos no Conceito 2.

Conceito 2: Modelo Probit, Probabilidades Previstas e Efeitos Estimados: Suponha que Y seja uma variável binária. O modelo

$Y= \beta_0 + \beta_1 + X_1 + \beta_2 X_2 + \dots + \beta_k X_k + u$

com

$P(Y = 1 \vert X_1, X_2, \dots ,X_k) = \Phi(\beta_0 + \beta_1 + X_1 + \beta_2 X_2 + \dots + \beta_k X_k)$

é o modelo Probit populacional com múltiplos regressores $X_1, X_2, \dots, X_k e \Phi(\cdot)$ é a função de distribuição normal padrão cumulativa.

A probabilidade prevista de que $Y=1$ dado $X_1, X_2, \dots, X_k$ pode ser calculada em duas etapas:

Calcule $z = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k$
Procure $\Phi(z)$ chamando pnorm() .

$\beta_j$ é o efeito sobre $z$ de uma mudança de uma unidade no regressor $X_j$ , mantendo constantes todos os outros $k-1$ regressores.

O efeito na probabilidade prevista de uma mudança num regressor pode ser calculado como no Conceito Chave 8.1.

Em R , os modelos Probit podem ser estimados usando a função glm() do pacote stats . Usando a família de argumentos , especificamos que queremos usar uma função de ligação Probit.

# estimate the simple probit model
denyprobit <- glm(deny ~ pirat, 
                  family = binomial(link = "probit"), 
                  data = HMDA)

coeftest(denyprobit, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -2.19415    0.18901 -11.6087 < 2.2e-16 ***
## pirat        2.96787    0.53698   5.5269 3.259e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O modelo estimado é

\[ \begin{align} \widehat{P(deny\vert P/I \ ratio}) = \Phi(-\underset{(0,19)}{2,19} + \underset{(0,54)}{2,97} P/I \ ratio). \tag{2} \end{align} \]

Tal como no modelo de probabilidade linear, descobrimos que a relação entre a probabilidade de recusa e o razão pagamentos/rendimento é positiva e que o coeficiente correspondente é altamente significativo.

O trecho de código a seguir reproduz o modelo probit

# plot data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Probit Model of the Probability of Denial, Given P/I Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.85)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add estimated regression line
x <- seq(0, 3, 0.01)
y <- predict(denyprobit, list(pirat = x), type = "response")

lines(x, y, lwd = 1.5, col = "steelblue")

A função de regressão estimada tem um formato de “S” esticado, típico para o CDF de uma variável aleatória contínua com PDF simétrica como a de uma variável aleatória normal. A função é claramente não linear e se achata para valores grandes e pequenos de $P/I \ ratio$. A forma funcional também garante que as probabilidades condicionais previstas de uma negação estejam entre $0$ e $1$ .

Usamos a função predict() para calcular a mudança prevista na probabilidade de negação quando a $P/I ratio$ é aumentada de $0.3$ para $0.4$.

# 1. compute predictions for P/I ratio = 0.3, 0.4
predictions <- predict(denyprobit, 
                       newdata = data.frame("pirat" = c(0.3, 0.4)),
                       type = "response")

# 2. Compute difference in probabilities
diff(predictions)

##          2 
## 0.06081433

Descobrimos que se prevê que um aumento na razão pagamento/rendimento de $0.3$ para $0.4$ aumente a probabilidade de negação em aproximadamente $6.1%$ .

Continuamos usando um modelo Probit aumentado para estimar o efeito da raça na probabilidade de negação de um pedido de hipoteca.

Expandindo o Modelo

#As a remainder, we have renamed the variable "afam" to "black" in Chapter 11.1
#colnames(HMDA)[colnames(HMDA) == "afam"] <- "black"
#We continue by using an augmented
denyprobit2 <- glm(deny ~ pirat + black, 
                   family = binomial(link = "probit"), 
                   data = HMDA)

coeftest(denyprobit2, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##              Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -2.258787   0.176608 -12.7898 < 2.2e-16 ***
## pirat        2.741779   0.497673   5.5092 3.605e-08 ***
## blackyes     0.708155   0.083091   8.5227 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A equação do modelo estimada é

\[ \begin{align} \widehat{P(deny\vert P/I \ ratio, black)} = \Phi (-\underset{(0,18)}{2,26} + \underset{(0,50)}{2,74} P/ I \ ratio + \underset{(0,08)}{0,71} black). \tag{11.6} \end{align} \]

Embora todos os coeficientes sejam altamente significativos, tanto os coeficientes estimados sobre a razão pagamentos/rendimento como o indicador para a ascendência afro-americana são positivos. Mais uma vez, os coeficientes são difíceis de interpretar, mas indicam que, em primeiro lugar, os afro-americanos têm uma maior probabilidade de recusa do que os candidatos brancos, mantendo constante o razão pagamentos/rendimento e, em segundo lugar, os requerentes com uma elevada razão pagamentos/rendimento enfrentam maior risco de ser rejeitado.

Qual é o tamanho da diferença estimada nas probabilidades de recusa entre dois requerentes hipotéticos com a mesma relação entre pagamentos e rendimentos? Como antes, podemos usar predict() para calcular essa diferença.

# 1. compute predictions for P/I ratio = 0.3
predictions <- predict(denyprobit2, 
newdata = data.frame("black" = c("no", "yes"), 
"pirat" = c(0.3, 0.3)),
type = "response")

# 2. compute difference in probabilities
diff(predictions)

##         2 
## 0.1578117

Neste caso, a diferença estimada nas probabilidades de negação é de cerca de $15.8\%$.

Logit

O Conceito 3 resume a função de regressão Logit.

Definição 3: Regressão Logística

A função de regressão Logit populacional é

\[\begin{align*} P(Y=1\vert X_1, X_2, \dots, X_k) =& \, F(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k) \\ =& \, \frac{1}{1+e^{-(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k)}}. \end{align*}\]

A ideia é semelhante à regressão Probit, exceto que um CDF diferente é usado:

\[F(x) = \frac{1}{1+e^{-x}}\]

é o CDF de uma variável aleatória padrão distribuída logisticamente.

Quanto à regressão Probit, não existe uma interpretação simples dos coeficientes do modelo e é melhor considerar as probabilidades previstas ou diferenças nas probabilidades previstas. Aqui, novamente, estatísticas t e intervalos de confiança baseados em aproximações normais de grandes amostras podem ser calculados como de costume.

É bastante fácil estimar um modelo de regressão Logit usando R .

denylogit <- glm(deny ~ pirat,
                 family = binomial(link = "logit"),
                 data = HMDA)

coeftest(denylogit, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -4.02843    0.35898 -11.2218 < 2.2e-16 ***
## pirat        5.88450    1.00015   5.8836 4.014e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O bloco de código subsequente reproduz a comparação do modelo logit com o modelo probit

# plot data
plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Probit and Logit Models of the Probability of Denial, Given P/I Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.9)

# add horizontal dashed lines and text
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Mortgage denied")
text(2.5, -0.1, cex= 0.8, "Mortgage approved")

# add estimated regression line of Probit and Logit models
x <- seq(0, 3, 0.01)
y_probit <- predict(denyprobit, list(pirat = x), type = "response")
y_logit <- predict(denylogit, list(pirat = x), type = "response")

lines(x, y_probit, lwd = 1.5, col = "steelblue")
lines(x, y_logit, lwd = 1.5, col = "black", lty = 2)

# add a legend
legend("topleft",
       horiz = TRUE,
       legend = c("Probit", "Logit"),
       col = c("steelblue", "black"), 
       lty = c(1, 2))

Ambos os modelos produzem estimativas muito semelhantes da probabilidade de um pedido de hipoteca ser negado, dependendo da relação entre pagamento e rendimento do requerente.

####Expandindo o modelo

Seguindo o livro, estendemos o modelo Logit simples de negação de hipoteca com o regressor adicional $black$

# estimate a Logit regression with multiple regressors
denylogit2 <- glm(deny ~ pirat + black, 
                  family = binomial(link = "logit"), 
                  data = HMDA)

coeftest(denylogit2, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -4.12556    0.34597 -11.9245 < 2.2e-16 ***
## pirat        5.37036    0.96376   5.5723 2.514e-08 ***
## blackyes     1.27278    0.14616   8.7081 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nós obtemos

\[ \begin{align} \widehat{P(deny=1 \vert ratio \ P/I , black)} = F(-\underset{(0,35)}{4,13} + \underset{(0,96)}{5,37} P/ I \ ratio + \underset{(0,15)}{1,27} black). \tag{3} \end{align} \]

Quanto ao modelo Probit (11.6) todos os coeficientes do modelo são altamente significativos e obtemos estimativas positivas para os coeficientes de $ratio \ P/I$ e $black$. Para comparação, calculamos a probabilidade prevista de negação para dois candidatos hipotéticos que diferem em raça e possui $ratio \ P/I = 0.3$.

# 1. compute predictions for P/I ratio = 0.3
predictions <- predict(denylogit2, 
                       newdata = data.frame("black" = c("no", "yes"), 
                                            "pirat" = c(0.3, 0.3)),
                       type = "response")
predictions

##          1          2 
## 0.07485143 0.22414592

# 2. Compute difference in probabilities
diff(predictions)

##         2 
## 0.1492945

Descobrimos que o candidato branco enfrenta uma probabilidade de negação de apenas $7.5\%$ , enquanto o afro-americano é rejeitado com uma probabilidade de $22.4\%$ , uma diferença de $14.9$ pontos percentuais.

Comparação dos Modelos

O modelo Probit e o modelo Logit fornecem apenas aproximações para a função de regressão populacional desconhecida $E(Y\vert X)$. Não é óbvio como decidir qual modelo usar na prática. O modelo de probabilidade linear tem a clara desvantagem de não ser capaz de capturar a natureza não linear da função de regressão populacional e pode prever probabilidades fora do intervalo $[0,1]$ . Os modelos Probit e Logit são mais difíceis de interpretar, mas capturam melhor as não linearidades do que a abordagem linear: ambos os modelos produzem previsões de probabilidades que estão dentro do intervalo $[0,1]$. As previsões dos três modelos são frequentemente próximas umas das outras. O livro sugere usar o método mais fácil de usar no software estatístico de sua escolha. Como vimos, é igualmente fácil estimar os modelos Probit e Logit usando R. Não podemos, portanto, dar nenhuma recomendação geral sobre qual método usar.

Estimação e Inferência nos Modelos Logit e Probit

Estimativa de Máxima Verossimilhança

No EMV procuramos estimar os parâmetros desconhecidos escolhendo-os de forma que a probabilidade de obtenção da amostra observada seja maximizada. Essa probabilidade é medida por meio da função de verossimilhança, a distribuição de probabilidade conjunta dos dados tratados em função dos parâmetros desconhecidos. Dito de outra forma, as estimativas de máxima verossimilhança dos parâmetros desconhecidos são os valores que resultam num modelo com maior probabilidade de produzir os dados observados.

Como as estimativas de máxima verossimilhança são normalmente distribuídas em grandes amostras, a inferência estatística para coeficientes em modelos não lineares como regressão Logit e Probit pode ser feita usando as mesmas ferramentas usadas para modelos de regressão linear: podemos calcular os intervalos de confiança das estatísticas-$t$.

Muitos pacotes de software usam um algoritmo EMV para estimativa de modelos não lineares. A função glm() usa um algoritmo denominado mínimos quadrados reponderados iterativamente .

Medidas de ajuste

É importante estar ciente de que os usuais $\bar{R}^2$ e $R^2$ são inválidos para modelos de regressão não linear. A razão para isto é simples: ambas as medidas assumem que a relação entre a(s) variável(es) dependente(s) e explicativa(s) é linear. Isto obviamente não se aplica aos modelos Probit e Logit. Assim não precisa estar entre $0$ e $1$ e não há interpretação significativa. No entanto, o software estatístico por vezes reporta estas medidas de qualquer maneira.

Existem muitas medidas de ajuste para modelos de regressão não linear e não há consenso sobre qual delas deve ser relatada. A situação é ainda mais complicada porque não existe uma medida de ajuste que seja geralmente significativa. Para modelos com uma variável de resposta binária como deny pode-se usar a seguinte regra: Se $Y_i = 1$ e $\widehat{P(Y_i|X_{i1}, \dots, X_{ik})} > 0.5$ ou se $Y_i = 0$ e $\widehat{P(Y_i|X_{i1}, \dots, X_{ik})} < 0.5$ , considere o $Y_i$ como previsto corretamente. Caso contrário, diz-se que $Y_i$ foi previsto incorretamente. A medida de ajuste é a parcela de observações previstas corretamente. A desvantagem de tal abordagem é que ela não reflete a qualidade da previsão: se $\widehat{P(Y_i = 1|X_{i1}, \dots, X_{ik}) = 0.51}$ ou $\widehat{P( Y_i =1|X_{i1}, \dots, X_{ik}) = 0.99}$ não é refletido, apenas prevemos $Y_i = 1$.

Uma alternativa a estas últimas são as chamadas medidas pseudo-$R^2$. Para medir a qualidade do ajuste, essas medidas comparam o valor da (log-)verossimilhança maximizada do modelo com todos os regressores (o modelo completo) com a probabilidade de um modelo sem regressores (modelo nulo, regressão em uma constante).

Por exemplo, considere uma regressão Probit. O $\text{pseudo-}R^2$ é dado por $\text{pseudo-}R^2 = 1 - \frac{\ln(f^{max}_{completo})}{\ln(f^{max}_{null})}$ onde $f^{max}_j \in [0,1]$ denota a probabilidade maximizada para o modelo $j$.

O raciocínio por trás disso é que a probabilidade maximizada aumenta à medida que regressores adicionais são adicionados ao modelo, de forma semelhante à diminuição do SSR quando regressores são adicionados em um modelo de regressão linear. Se o modelo completo tiver uma probabilidade maximizada semelhante à do modelo nulo, o modelo completo não melhora realmente um modelo que usa apenas as informações na variável dependente, então $\text{pseudo-}R^2 \approx 0$. Se o modelo completo se ajustar muito bem aos dados, a probabilidade maximizada deve ser próxima de 1 , de modo que $\ln(f^{max}_{full}) \approx 0$ e $\text{pseudo-}R^2 \approx 1$.

summary() não reporta $\text{pseudo-}R^2$ para modelos estimados por glm() mas podemos usar as entradas desvio residual (deviance) e desvio nulo (null.deviance). Estes são calculados como

\[\text{deviance} = -2 \times \left[\ln(f^{max}_{saturated}) - \ln(f^{max}_{full}) \right]\]

\[\text{null deviance} = -2 \times \left[\ln(f^{max}_{saturated}) - \ln(f^{max}_{null}) \right]\]

onde $f^{max}_{saturated}$ é a probabilidade maximizada para um modelo que assume que cada observação tem seu próprio parâmetro (há $n+1$ parâmetros a serem estimados, o que leva a um ajuste perfeito). Para modelos com uma variável dependente binária, afirma-se que

\[\text{pseudo-}R^2 = 1 - \frac{\text{deviance}}{\text{null deviance}} = 1- \frac{\ln(f^{max}_{full})}{ \ln(f^{max}_{null})}\]

Agora calculamos o $\text{pseudo-}R^2$ para o modelo Probit aumentado de negação de hipotecas.

# compute pseudo-R2 for the probit model of mortgage denial
pseudoR2 <- 1 - (denyprobit2$deviance) / (denyprobit2$null.deviance)
pseudoR2

## [1] 0.08594259

Outra forma de obter o $\text{pseudo-}R^2$ é estimar o modelo nulo usando glm() e extrair as probabilidades logarítmicas maximizadas para o modelo nulo e completo usando a função logLik().

# compute the null model
denyprobit_null <- glm(formula = deny ~ 1, 
                       family = binomial(link = "probit"), 
                       data = HMDA)

# compute the pseudo-R2 using 'logLik'
1 - logLik(denyprobit2)[1]/logLik(denyprobit_null)[1]

## [1] 0.08594259

Outros modelos com variável dependente limitada

library(sampleSelection)
library(stargazer)
library(censReg)

Modelo de regressão censurado (Tobit)

O modelo tobit, também chamado de modelo de regressão censurado, é projetado para estimar relações lineares entre variáveis quando há censura à esquerda ou à direita na variável dependente (também conhecida como censura de baixo e de cima, respectivamente). A censura superior ocorre quando os casos com um valor igual ou superior a algum limite, todos assumem o valor desse limite, de modo que o valor verdadeiro pode ser igual ao limite, mas também pode ser maior. No caso de censura por baixo, os valores que ficam em ou abaixo de algum limite são censurados.

O conjunto de dados Mroz87 são dados de 1975 sobre o salário das mulheres casadas e a participação na força de trabalho, extraídos de um conhecido artigo de 1987 de Thomas Mroz. Mroz retirou os seus dados do Panel Study of Income Dynamics (PSID), um importante inquérito nacional. O conjunto de dados está disponível como parte do pacote sampleSelection, que você deve ter instalado e carregado. Para ver quais variáveis estão nos dados, você utilizar o comando help(Mroz87) .

data(Mroz87)

Mroz87_lf = Mroz87 %>% 
  filter(lfp==1)

Podemos observar o histograma da distribuição de horas trabalhadas pelas mulheres

Mroz87 %>% 
  ggplot(aes(hours)) +
  geom_histogram()

MQO

# Estimate OLS regression
ols <- lm(hours ~ nwifeinc + educ + exper + I(exper^2) + age + kids5 + kids618,
          data = Mroz87)

Tobit

# Estimate Tobit model using censReg package
library(censReg)
TobitRes <- censReg(hours ~ nwifeinc + educ + exper + I(exper^2) + age + kids5 + kids618,
                    data = Mroz87)

# Notes and interpretation: 325 observations are censored from the left.
# Interpretation of the coefficients:
# The women who have a higher family income (nwifeinc), are older (age), and 
#have a kid less than 6 years old have fewer hours worked in 1975.
# The women who have more education (educ) have more hours worked in 1975.

stargazer(ols, TobitRes,  
          type = "text")

## 
## =======================================================
##                             Dependent variable:        
##                     -----------------------------------
##                                    hours               
##                               OLS            censored  
##                                             regression 
##                               (1)               (2)    
## -------------------------------------------------------
## nwifeinc                    -3.447           -8.814**  
##                             (2.544)           (4.459)  
##                                                        
## educ                       28.761**          80.646*** 
##                            (12.955)          (21.583)  
##                                                        
## exper                      65.673***        131.564*** 
##                             (9.963)          (17.279)  
##                                                        
## I(exper2)                  -0.700**          -1.864*** 
##                             (0.325)           (0.538)  
##                                                        
## age                       -30.512***        -54.405*** 
##                             (4.364)           (7.419)  
##                                                        
## kids5                     -442.090***       -894.022***
##                            (58.847)          (111.878) 
##                                                        
## kids618                     -32.779           -16.218  
##                            (23.176)          (38.641)  
##                                                        
## logSigma                                     7.023***  
##                                               (0.037)  
##                                                        
## Constant                 1,330.482***        965.305** 
##                            (270.785)         (446.436) 
##                                                        
## -------------------------------------------------------
## Observations                  753               753    
## R2                           0.266                     
## Adjusted R2                  0.259                     
## Log Likelihood                              -3,819.095 
## Akaike Inf. Crit.                            7,656.189 
## Bayesian Inf. Crit.                          7,697.806 
## Residual Std. Error   750.179 (df = 745)               
## F Statistic         38.495*** (df = 7; 745)            
## =======================================================
## Note:                       *p<0.1; **p<0.05; ***p<0.01

# First, the Tobit coefficient estimates have the same sign as the 
#corresponding OLS estimates,
# and the statistical significance of the estimates is similar. 
#(Possible exceptions are the coefficients
# on nwifeinc and kidsge6, but the t statistics have similar magnitudes.)
# Second, though it is tempting to compare the magnitudes of the OLS and Tobit 
#estimates, this is not very informative.
# We must be careful not to think that, because the Tobit coefficient on 
#kidslt6 is roughly twice that of the OLS coefficient

# Partial effects
Tobit_margEff <- margEff(TobitRes)
summary(Tobit_margEff)

##            Marg. Eff. Std. Error t value  Pr(>|t|)    
## nwifeinc     -5.32644    2.69073 -1.9796 0.0481217 *  
## educ         48.73409   12.96341  3.7594 0.0001837 ***
## exper        79.50423   10.30497  7.7151 3.886e-14 ***
## I(exper^2)   -1.12651    0.32326 -3.4848 0.0005213 ***
## age         -32.87692    4.45770 -7.3753 4.383e-13 ***
## kids5      -540.25683   66.62393 -8.1091 2.220e-15 ***
## kids618      -9.80053   23.36134 -0.4195 0.6749580    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Interpretation of the marginal effects:
# The a woman's family income (nwifeinc) increases by one-unit, her work 
#hours decreases by 5.32, on average.
# # If a woman has one more year of education, her work hours increases by 48.73 hours.
# If a woman is older by one year (age), her work hours reduces by 32.87 hours.
# If a woman has one extra kid who is less than 6 years old, her work 
#hours decrease by 540 hours, on average (It is a large effect.)

Modelo de Contagem (Regressão de Poisson)

#Example 17.3
data(crime1, package='wooldridge')
attach(crime1)

crime1 %>% 
  ggplot(aes(narr86)) +
  geom_bar()

# The dependent variable, narr86, is the number of times a man is arrested 
#during 1986. This variable is zero for 1,970 of the 2,725 men in the sample, and only eight values of
#narr86 are greater than five. Thus, a Poisson regression model is more appropriate than a linear 
#regression model.

# Estimate linear model
regLPM  <-  lm(narr86 ~ pcnv + avgsen + tottime + ptime86+ qemp86 + inc86 + black + hispan + born60,
               data = crime1)

# Estimate the poisson model
regPoisson <- glm(narr86 ~ pcnv + avgsen + tottime + ptime86+ qemp86 + inc86 + black + hispan + born60,
                  data = crime1,
                  family = poisson)


stargazer(regLPM, regPoisson,
          type = "text",
          keep.stat = "n")

## 
## =========================================
##                  Dependent variable:     
##              ----------------------------
##                         narr86           
##                   OLS          Poisson   
##                   (1)            (2)     
## -----------------------------------------
## pcnv           -0.132***      -0.402***  
##                 (0.040)        (0.085)   
##                                          
## avgsen           -0.011        -0.024    
##                 (0.012)        (0.020)   
##                                          
## tottime          0.012         0.024*    
##                 (0.009)        (0.015)   
##                                          
## ptime86        -0.041***      -0.099***  
##                 (0.009)        (0.021)   
##                                          
## qemp86         -0.051***       -0.038    
##                 (0.014)        (0.029)   
##                                          
## inc86          -0.001***      -0.008***  
##                 (0.0003)       (0.001)   
##                                          
## black           0.327***      0.661***   
##                 (0.045)        (0.074)   
##                                          
## hispan          0.194***      0.500***   
##                 (0.040)        (0.074)   
##                                          
## born60           -0.022        -0.051    
##                 (0.033)        (0.064)   
##                                          
## Constant        0.577***      -0.600***  
##                 (0.038)        (0.067)   
##                                          
## -----------------------------------------
## Observations     2,725          2,725    
## =========================================
## Note:         *p<0.1; **p<0.05; ***p<0.01

# The OLS and Poisson coefficients are not directly comparable, they have very 
#different meanings.
# For example, the OLS coefficient on pcnv implies that, if pcnv changes 0.10, 
#the expected number of arrests
# falls by 0.013 (pcnv is the proportion of prior arrests that led to conviction).

# The Poisson coefficient implies that if pcnv changes 0.10, it reduces expected 
#arrests by about 4% (0.40*0.10 = .040),
# and we multiply this by 100 to get the percentage effect.

# The Poisson coefficient on black implies that other factors being equal, the 
#expected number of arrests for a black man is estimated to be about 
#100*[exp(0.66) - 1] = 93.7% higher than for a white man with the same values #
#for the other explanatory variables.


n <- nobs(regPoisson)
y_hat <- fitted(regPoisson)
u_hat2 <- (crime1$narr86 - y_hat)^2

(sigma2 <- (sum(u_hat2/y_hat))/(n - 10))

## [1] 1.516788

(sigma <- sqrt(sigma2))

## [1] 1.23158

# sigma = 1.232 implies that the standard errors for Poisson regression should 
#be inflated by sigma, i.e., should be 23% larger. 
#For example, s.e.(tottime) = 1.23*0.015 = 0.0185.

Modelos com seleção de amostra (Heckit)

O viés de seleção de amostra em regressões OLS surge quando o processo pelo qual as observações são selecionadas na amostra depende da variável dependente ($Y$) de uma forma que cria uma correlação entre o termo de erro e os regressores. Um exemplo clássico é a estimativa de uma equação de salários em que o salário é uma função de competências observáveis (por exemplo, educação), quando a decisão de trabalhar por remuneração depende do próprio salário esperado.

Para ilustrar, consideremos um cenáriosimplificado, em que as pessoas só optam por trabalhar por remuneração se o seu salário esperado exceder algum valor limite constante. Caso contrário, o trabalho não vale a pena para eles e recorrem a outros recursos (cônjuge, pais, etc.). Em seguida, a amostra de trabalhadores da força de trabalho, para os quais observamos os rendimentos, é selecionada entre aqueles da população que recebem salários acima desse limiar. Entre os trabalhadores altamente qualificados, a maioria ganharia acima do limiar de qualquer maneira, pelo que estão bastante bem representados na amostra. Mas entre os trabalhadores com menos escolaridade, geralmente com baixos salários, a amostra é mais selectiva e não aleatória: apenas aqueles com rendimentos invulgarmente elevados, dados os seus baixos níveis de educação, ultrapassam o limiar e são observados na amostra de trabalhadores remunerados. Neste cenário, a educação estaria negativamente correlacionada com o termo de erro, o que distorce o coeficiente OLS da educação para baixo.

É possível estimar um modelo que corrija o viés de seleção da amostra que estaria presente nos MQO. Essencialmente, o que precisamos é de um modelo de processo de seleção. Então, uma vez estimados os determinantes da seleção, podemos usar isso para corrigir o viés na regressão. Nosso foco neste tutorial está na mecânica de como estimar esse tipo de modelo em R, o que é simples.

Estimaremos uma equação logarítmica dos salários para mulheres casadas, na qual o logaritmo natural dos salários é uma função da educação, experiência, experiência ao quadrado e uma variável dummy para viver numa cidade grande. A questão aqui é que só podemos observar o salário das mulheres que trabalham por um salário. Uma proporção substancial de mulheres casadas em 1975 estava fora da força de trabalho, pelo que uma estimativa OLS da equação salarial seria estimada utilizando uma amostra sujeita a uma selecção considerável da amostra, resultando potencialmente num enviesamento da natureza discutida acima.

Como temos observações no conjunto de dados sobre não participantes (que não trabalham remunerados), podemos implementar modelos que corrigem o processo de seleção. Observe também que criamos uma nova variável que é uma variável dummy para a mulher que tem filhos menores de 18 anos. Na verdade, usaremos isso como um instrumento para a participação na força de trabalho.

Vamos executar três regressões para estimar a equação logarítmica do salário e comparar os resultados. O primeiro é simplesmente MQO, utilizando a amostra de participantes na força de trabalho, para os quais podemos observar o salário. Aqui está o código, que deve parecer bastante familiar. Certifique-se de entender tudo:

data(Mroz87)

MQO

Mroz87_lf = Mroz87 %>% 
  filter(lfp==1)

# OLS: log wage regression on LF participants only
ols1 = lm(log(wage) ~ educ + exper + I( exper^2 ), data=Mroz87_lf)

Heckit

Agora vamos examinar dois métodos alternativos para estimar a equação salarial com seleção. O primeiro é o chamado modelo de correção em duas etapas de Heckman (heckit), nomeado em homenagem ao economista ganhador do Prêmio Nobel, James Heckman. A ideia de Heckman era tratar o problema de seleção como se fosse um problema de variável omitida. Uma equação probit de primeira fase estima o processo de selecção (quem está na força de trabalho?), e os resultados dessa equação são utilizados para construir uma variável que capta o efeito de seleção na equação salarial. Esta variável de correção é chamada de razão de Mills inversa.

# Two-step estimation with LFP selection equation
heck1 = heckit(lfp ~ educ + age + exper + I(exper^2) + nwifeinc + kids5 + kids618,
               log(wage) ~ educ + exper + I(exper^2), data=Mroz87)

A configuração segue o espírito do padrão lm, mas neste caso observe que temos duas equações, separadas por vírgula.

A primeira equação do comando é o processo de seleção, que é estimado como um probit: lfp ~ age + I( age^2 ) + faminc + kids + huswage + educ. A variável dependente é sempre a variável binária a ser selecionada na amostra – neste caso, é lfp , que é uma variável dummy igual a um se a mulher estiver participando da força de trabalho (e, portanto, estiver sendo paga e observada durante o equação salarial). Depois vem uma vírgula, seguida por:
A segunda equação, que é a regressão salarial de interesse: log(wage) ~ exper + I( exper^2 ) + educ + city.

A equação de seleção deve incluir como regressores variáveis que possam afetar o processo de seleção. No nosso exemplo, queremos variáveis que possam afetar de forma plausível se uma mulher casada estaria ou não na força de trabalho. Examine as variáveis incluídas aqui, certifique-se de saber quais são e veja se consegue entender por que elas seriam incluídas na equação lfp .

Em princípio, a equação de seleção e a equação de salários poderiam ter exactamente o mesmo conjunto de regressores. Mas isso geralmente não é uma boa ideia. Para obter uma boa estimativa do modelo de seleção, é muito desejável ter pelo menos uma variável na equação de seleção que atue como um instrumento : isto é, uma variável que se espera que afete o processo de seleção, mas não o processo salarial, exceto por meio de seleção. Neste exemplo, a variável filhos pode desempenhar esse papel, se pudermos assumir que as mulheres com filhos podem ter maior probabilidade de serem mães que ficam em casa, mas para as mães que trabalham, ter filhos não afetaria a sua taxa de pagamento por hora.

Estimação por Máxima Verossimilhança

Um método estatístico alternativo para obter o mesmo resultado do modelo de Heckman é estimar o modelo de seleção usando máxima verossimilhança (ML) para ambas as equações simultaneamente. Discutiremos isso mais adiante em aula. Por enquanto, observe que a sintaxe do comando é virtualmente idêntica à do heckit:

# ML estimation of selection model
ml1 = selection(lfp ~ educ + age + exper + I(exper^2 ) + nwifeinc + kids5 + kids618,
               log(wage) ~ educ + exper + I(exper^2), data=Mroz87)

Comparação

# stargazer table
stargazer(ols1, heck1, ml1,    
          #se=list(cse(ols1),NULL,NULL), 
          title="Married women's wage regressions", type="text", 
          df=FALSE, digits=4)

## 
## Married women's wage regressions
## ==============================================================
##                                Dependent variable:            
##                     ------------------------------------------
##                                     log(wage)                 
##                        OLS         Heckman        selection   
##                                   selection                   
##                        (1)           (2)             (3)      
## --------------------------------------------------------------
## educ                0.1075***     0.1091***       0.1084***   
##                      (0.0141)     (0.0155)        (0.0149)    
##                                                               
## exper               0.0416***     0.0439***       0.0428***   
##                      (0.0132)     (0.0163)        (0.0149)    
##                                                               
## I(exper2)           -0.0008**     -0.0009*        -0.0008**   
##                      (0.0004)     (0.0004)        (0.0004)    
##                                                               
## Constant            -0.5220***    -0.5781*        -0.5527**   
##                      (0.1986)     (0.3050)        (0.2604)    
##                                                               
## --------------------------------------------------------------
## Observations           428           753             753      
## R2                    0.1568       0.1569                     
## Adjusted R2           0.1509       0.1490                     
## Log Likelihood                                    -832.8851   
## rho                                0.0486      0.0266 (0.1471)
## Inverse Mills Ratio            0.0323 (0.1336)                
## Residual Std. Error   0.6664                                  
## F Statistic         26.2862***                                
## ==============================================================
## Note:                              *p<0.1; **p<0.05; ***p<0.01

Como você pode ver, as correções de seleção não tiveram grande efeito em nenhum dos parâmetros de interesse. Os retornos estimados para educação e experiência são um pouco maiores nas regressões de seleção, mas a diferença é muito pequena. Neste caso, parece que a equação OLS na amostra selecionada pode não ter sido muito tendenciosa, afinal.

No painel abaixo dos coeficientes, observe que há uma estimativa para rho nos modelos de seleção: colunas (2) e (3). Esta é uma estimativa da correlação dos erros entre as equações de seleção e de salário. No painel inferior, o coeficiente estimado no índice inverso de Mills é dado para o modelo de Heckman. O facto de não ser estatisticamente diferente de zero é consistente com a ideia de que o viés de seleção não foi um problema sério neste caso.

A tabela final stargazerdo script inclui a opção selection.equation=TRUE, o que significa que a tabela mostrará o modelo probit de primeiro estágio para $lfp$ , se você estiver interessado em quais fatores afetam o próprio processo de seleção.