Aula de laboratório - Escolha Qualitativa

Autor

Rafael Bressan

Data de Publicação

7 de agosto de 2023

library(wooldridge)
library(fixest)
library(modelsummary)
library(ggplot2)
data("mroz")

Modelo de Probabilidade Linear

Vamos analizar participação na força de trabalho de mulheres casadas com os dados mroz do pacote wooldridge. Nossa variável dependente será o status de participação na força de trabalho, $inlf\in\{0, 1\}$.

Se uma variável dummy é usada como variável dependente, nós podemos usar MQO para estimar a relação desta com os regressores normalmente. A probabilidade condicional de ocorrência do evento é:

\[P(y=1|x)=E[y|x]=\beta_0 + \beta_1 x\]

A interpretação dos coeficientes é direta, um aumento de uma unidade em $x$ provoca um aumento médio na probabilidade do evento de $\beta_1$. Modelos MPL sofrem de heterocedasticidade e, portanto, devemos utilizar erros padrão robustos.

Cheque a estrutura dos dados com uma tabela de estatísticas descritivas.

	Prompt
	`crie uma tabela de estatísticas descritivas para todas as variáveis do conjunto de dados mroz utilizando a função datasummary do pacote modelsummary. As estatísticas devem ser número de observações, número de valores únicos, média e desvio padrão.`

Código

datasummary(All(mroz) ~ Ncol + NUnique + Mean + SD, data = mroz)

	Ncol	NUnique	Mean	SD
inlf	(N = 753)	2	0.57	0.50
hours	(N = 753)	306	740.58	871.31
kidslt6	(N = 753)	4	0.24	0.52
kidsge6	(N = 753)	9	1.35	1.32
age	(N = 753)	31	42.54	8.07
educ	(N = 753)	13	12.29	2.28
wage	(N = 753)	374	4.18	3.31
repwage	(N = 753)	164	1.85	2.42
hushrs	(N = 753)	422	2267.27	595.57
husage	(N = 753)	31	45.12	8.06
huseduc	(N = 753)	15	12.49	3.02
huswage	(N = 753)	697	7.48	4.23
faminc	(N = 753)	613	23080.59	12190.20
mtr	(N = 753)	22	0.68	0.08
motheduc	(N = 753)	9	9.25	3.37
fatheduc	(N = 753)	9	8.81	3.57
unem	(N = 753)	7	8.62	3.11
city	(N = 753)	2	0.64	0.48
exper	(N = 753)	40	10.63	8.07
nwifeinc	(N = 753)	706	20.13	11.63
lwage	(N = 753)	374	1.19	0.72
expersq	(N = 753)	40	178.04	249.63

Agora faça a regressão da dummy de presença na força de trabalho em relação as seguintes variáveis: nwifeinc, educ, exper, exper^2, age, kidslt6, kidsge6.

	Prompt
	`faça uma regressão linear da variável inlf em relação as variáveis nwifeinc, educ, exper, exper^2, age, kidslt6, kidsge6. Os erros padrão reportados devem ser robustos. Utilize o pacote fixest.`

Código

fml <- inlf~nwifeinc + educ + exper + I(exper^2) + age + kidslt6 + kidsge6
mpl <- feols(fml,
    data = mroz,
    vcov = "HC1"
)
msummary(mpl)

	(1)
(Intercept)	0.586
	(0.152)
nwifeinc	-0.003
	(0.002)
educ	0.038
	(0.007)
exper	0.039
	(0.006)
I(exper^2)	-0.001
	(0.000)
age	-0.016
	(0.002)
kidslt6	-0.262
	(0.032)
kidsge6	0.013
	(0.014)
Num.Obs.	753
R2	0.264
R2 Adj.	0.257
AIC	863.8
BIC	900.8
RMSE	0.42
Std.Errors	Heteroskedasticity-robust

Por construção, um modelo MPL pode render previsão de probabilidade fora do intervalo $\left[0,1\right]$ para uma certa combinação de regressores (que pode ser realista ou não).

Como exemplo, supomos duas mulheres, a primeira com 20 anos, sem experiência, 5 anos de educação, 2 crianças com menos de 6 anos e renda de $ 100 (em milhares). A segunda mulher possui 52 anos, 30 anos de experiência e 17 de educação, sem crianças e renda.

Faça a previsão de probabilidade de estarem no mercado de trabalho para ambas hipotéticas mulheres.

	Prompt
	`com base na regressão anterior, faça a previsão de probabilidade de estarem no mercado de trabalho para duas mulheres hipotéticas. A primeira com 20 anos, sem experiência, 5 anos de educação, 2 crianças com menos de 6 anos e renda de $ 100 (em milhares). A segunda mulher possui 52 anos, 30 anos de experiência e 17 de educação, sem crianças e renda.`

xpred <- data.frame(
    nwifeinc = c(100, 0),
    educ = c(5, 17),
    exper = c(0, 30),
    age = c(20, 52),
    kidslt6 = c(2, 0),
    kidsge6 = c(0, 0)
)
predict(mpl, newdata = xpred)

[1] -0.4104582  1.0428084

Modelos Logit e Probit

Estes são modelos especializados em variáveis dependentes binárias e, portanto, não sofrem do problema de predição de probabilidade fora do intervalo $\left[0,1\right]$. A um indexador linear, estes modelos aplicam uma função (não-linear) que garante $y\in \left[0,1\right]$. Esta é chamada de função link.

no modelo probit a função link é $G(z)=\Phi(z)$, a função de distribuição da Normal, e;
no modelo logit a função link é $G(z)=\Lambda(z)=\frac{\exp(z)}{1+\exp(z)}$, a função de distribuição Logística.

Desta forma, a probabilidade do evento é definida como:

\[P(y=1|x)=G(\beta_0 + \beta_1 x)\]

No R estimamos esses modelos através da função glm, de Generalized Linear Models, especificando o argumento family. Vejamos um exemplo.

logit <- glm(fml,
    data = mroz,
    family = binomial(link = logit)
)

probit <- glm(fml,
    data = mroz,
    family = binomial(link = probit)
)

msummary(list(Logit = logit, Probit = probit))

	Logit	Probit
(Intercept)	0.425	0.270
	(0.860)	(0.508)
nwifeinc	-0.021	-0.012
	(0.008)	(0.005)
educ	0.221	0.131
	(0.043)	(0.025)
exper	0.206	0.123
	(0.032)	(0.019)
I(exper^2)	-0.003	-0.002
	(0.001)	(0.001)
age	-0.088	-0.053
	(0.015)	(0.008)
kidslt6	-1.443	-0.868
	(0.204)	(0.118)
kidsge6	0.060	0.036
	(0.075)	(0.044)
Num.Obs.	753	753
AIC	819.5	818.6
BIC	856.5	855.6
Log.Lik.	-401.765	-401.302
F	21.785	25.501
RMSE	0.42	0.42

Uma forma visual de apresentar os resultados de uma regressão é através de um gráfico de coeficientes (coefficients plot). Este gráfico mostra a estimativa pontual acompanhada de uma barra denotando o intervalo de confiança (e.g. 90%). É muito utilizado em apresentações sobre a pesquisa ou estudo sendo realizada.

modelplot(list(Logit = logit, Probit = probit),
    conf_level = 0.9
)

Podemos fazer a previsão da probabilidade das mesmas mulheres anteriores estarem no mercado de trabalho.

predict(logit, newdata = xpred, type = "response")

          1           2 
0.005218002 0.950049117

predict(probit, newdata = xpred, type = "response")

          1           2 
0.001065043 0.959869044

As previsões de logit e probit são semelhantes, mas especialmente nas caudas elas podem diferirem de forma significativa. O comando predict também pode retornar o valor do índice linear ($\hat\beta_0 + \hat\beta_1 x_i$) quando o argumento for type = “link”

predict(logit, newdata = xpred, type = "link")

        1         2 
-5.250409  2.945474

predict(probit, newdata = xpred, type = "link")

        1         2 
-3.071468  1.749168

Simule valores para o índice linear entre -5 e 5 e plote qual seriam as respostas em probabilidades dos modelos logit e probit.

	Prompt
	`simule a resposta de um modelo logit e outro probit para valores do índice linear entre -5 e 5 e plote as respostas em probabilidades dos modelos utilizando o ggplot.`

Código

index <- seq(-5, 5, by = 0.5)
logit_response <- plogis(index)
probit_response <- pnorm(index)
df <- data.frame(index, logit_response, probit_response)

ggplot(df, aes(x = index)) +
    geom_line(aes(y = logit_response, color = "Logit")) +
    geom_line(aes(y = probit_response, color = "Probit")) +
    guides(color = guide_legend(title = "Modelo")) +
    labs(
        x = "Índice linear",
        y = "Probabilidade",
        title = "Resposta em Probabilidade"
    ) +
    theme_classic()

Efeitos Parciais

Os parâmetros de modelos não-lineares como logit e probit não possuem uma interpretação direta dos efeitos ceteris paribus do regressor sobre a variável dependente. Uma medida utilizada é o efeito parcial (ou efeito marginal) que é dependente do ponto $x_i$ em que se faz a análise. Este efeito marginal é equivalente a inclinação da curva de resposta do modelo em um determinado ponto.

\[ \begin{align} \frac{\partial\hat y}{\partial x}&=\frac{\partial G(x\hat\beta)}{\partial x}\\ &=\hat\beta \cdot g(x\hat\beta) \end{align} \]

Onde é possível verificar que o efeito marginal depende tanto do valor estimado $\hat\beta$ quanto do ponto onde se faz a análise, $x$. O fato de este efeito diferir por valor do regressor faz com que seja difícil apresentar os resultados de uma maneira concisa. Adota-se, em geral, duas medidas de efeito parcial:

Efeito Parcial na Média (PEA): $PEA(X_j)=g(\bar{x} \beta) \beta_j$
Efeito Parcial Médio (APE): $APE(X_j)=\frac{1}{n} \sum_{i=1}^N g(x_i \beta) \beta_j$

Calcule manualmente o PEA do modelo logit

# Ponto médio dos regressores
regressors <- c("nwifeinc", "educ", "exper", "age", "kidslt6", "kidsge6")
x_mean <- mroz[, regressors] |>
    colMeans()
betas <- coef(logit)
# x_mean deve conter 1 na primeira posicao e exper^2 na quinta
x_mean <- c(1, x_mean[1:3], mean(mroz[, "exper"]^2), x_mean[4:6])
# xb
index_mean <- sum(x_mean * betas)
# Agora o PEA
print("PEA de cada regressor é:")

[1] "PEA de cada regressor é:"

dlogis(index_mean) * betas

  (Intercept)      nwifeinc          educ         exper    I(exper^2) 
 0.1034482321 -0.0051900534  0.0537773087  0.0500569282 -0.0007669166 
          age       kidslt6       kidsge6 
-0.0214030205 -0.3509498193  0.0146162143

# usando o pacote mfx
mfx::logitmfx(fml, data = mroz, atmean = TRUE)

Call:
mfx::logitmfx(formula = fml, data = mroz, atmean = TRUE)

Marginal Effects:
                 dF/dx   Std. Err.       z     P>|z|    
nwifeinc   -0.00519005  0.00204820 -2.5340  0.011278 *  
educ        0.05377731  0.01056074  5.0922 3.539e-07 ***
exper       0.05005693  0.00782462  6.3974 1.581e-10 ***
I(exper^2) -0.00076692  0.00024768 -3.0965  0.001959 ** 
age        -0.02140302  0.00353973 -6.0465 1.480e-09 ***
kidslt6    -0.35094982  0.04963897 -7.0700 1.549e-12 ***
kidsge6     0.01461621  0.01818832  0.8036  0.421625    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Calcule manualmente o APE do modelo logit

g_mean <- predict(logit) |>
    dlogis() |>
    mean()
print("APE de cada regressor é:")

[1] "APE de cada regressor é:"

betas * g_mean

  (Intercept)      nwifeinc          educ         exper    I(exper^2) 
 0.0759771297 -0.0038118135  0.0394965238  0.0367641056 -0.0005632587 
          age       kidslt6       kidsge6 
-0.0157193606 -0.2577536551  0.0107348186

# usando o pacote mfx
mfx::logitmfx(fml, data = mroz, atmean = FALSE)

Call:
mfx::logitmfx(formula = fml, data = mroz, atmean = FALSE)

Marginal Effects:
                 dF/dx   Std. Err.       z     P>|z|    
nwifeinc   -0.00381181  0.00153898 -2.4769  0.013255 *  
educ        0.03949652  0.00846811  4.6641 3.099e-06 ***
exper       0.03676411  0.00655577  5.6079 2.048e-08 ***
I(exper^2) -0.00056326  0.00018795 -2.9968  0.002728 ** 
age        -0.01571936  0.00293269 -5.3600 8.320e-08 ***
kidslt6    -0.25775366  0.04263493 -6.0456 1.489e-09 ***
kidsge6     0.01073482  0.01339130  0.8016  0.422769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1