De forma a clarear brevemente sua importância, os modelos Logit e Probit (abreviação de regressão logística e probabilística) nos auxiliam na inferência de probabilidade de ocorrência de eventos onde nossa variável dependente é binária (Y ocorre ou não ocorre), e nosso objetivo é compreender como outras variáveis influenciam a ocorrência ou não desses eventos (como nosso X prevê a ocorrência de Y).
Você pode estar se perguntando agora, mas por que não utilizar o modelo linear? Pois bem, em uma regressão linear, P(Y=1|x) é dado por uma especificação linear dos regressores, o que pode resultar em valores menores que 0 ou maiores que 1, que não fazem sentido com a interpretação probabilística dos parâmetros. Os modelos não lineares entretanto permitem isso, de forma que a média condicional de Y dado X seja expressa pela probabilidade de Y acontecer dado X: \(E(Y|X) = P(Y=1|X)\).
Nosso objetivo aqui permanece sendo encontrar os parâmetros desconhecidos (\(\theta\)), entretanto, uma vez que não trabalhamos mais com o modelo linear, passaremos a utilizar um novo estimador – maxima verossimilhança – que busca aproximar \(\hat{\theta}\) de \(\theta\) via maximização da sua distribuição “conhecida”.
Considere inlf (“no mercado de trabalho”) como uma variável binária que indica a participação no mercado de trabalho por uma mulher casada durante 1975: \(inlf = 1\) se a mulher relata ter trabalhado por um salário fora de casa em algum momento durante o ano, e zero caso contrário. Assumimos que a participação no mercado de trabalho depende de outras fontes de renda, incluindo os ganhos do marido (nwifeinc, medidos em milhares de dólares), anos de educação (educ), anos passados de experiência no mercado de trabalho (exper), idade, número de filhos menores de seis anos (kidslt6), e número de filhos entre 6 e 18 anos de idade (kidsge6). Usando os dados em MROZ de Mroz (1987), estimamos o seguinte modelo de probabilidade linear:
\[inlf= \beta_0 - \beta_1*nwifeinc + beta_2*educ + \beta_3*exper - \beta_4*exper^2 - \beta_5*age − \beta_6*kidslt6 + \beta_7*kidsge6\]
No R, temos:
options(scipen = 999) #desliga a notação científica
#packages-----
library(tidyverse) #analise de dados## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractrequire(mfx) #odds ratio## Carregando pacotes exigidos: mfx
## Carregando pacotes exigidos: sandwich
## Carregando pacotes exigidos: lmtest
## Carregando pacotes exigidos: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Carregando pacotes exigidos: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## Carregando pacotes exigidos: betareglibrary(wooldridge) #base de dados##
## Attaching package: 'wooldridge'
##
## The following object is masked from 'package:MASS':
##
## cement#LOGIT
mlogit<- glm(inlf ~ nwifeinc + educ + exper + expersq + age +
kidslt6 + kidsge6,
data = mroz,
family = binomial(link = "logit"))
summary(mlogit)##
## Call:
## glm(formula = inlf ~ nwifeinc + educ + exper + expersq + age +
## kidslt6 + kidsge6, family = binomial(link = "logit"), data = mroz)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.425452 0.860365 0.495 0.62095
## nwifeinc -0.021345 0.008421 -2.535 0.01126 *
## educ 0.221170 0.043439 5.091 0.00000035527344 ***
## exper 0.205870 0.032057 6.422 0.00000000013446 ***
## expersq -0.003154 0.001016 -3.104 0.00191 **
## age -0.088024 0.014573 -6.040 0.00000000153845 ***
## kidslt6 -1.443354 0.203583 -7.090 0.00000000000134 ***
## kidsge6 0.060112 0.074789 0.804 0.42154
## ---
## 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: 803.53 on 745 degrees of freedom
## AIC: 819.53
##
## Number of Fisher Scoring iterations: 4#PROBIT
mprobit <- glm(inlf ~ nwifeinc + educ + exper + expersq + age +
kidslt6 + kidsge6,
data = mroz,
family = binomial(link = "probit"))
summary(mprobit)##
## Call:
## glm(formula = inlf ~ nwifeinc + educ + exper + expersq + age +
## kidslt6 + kidsge6, family = binomial(link = "probit"), data = mroz)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2700736 0.5080782 0.532 0.59503
## nwifeinc -0.0120236 0.0049392 -2.434 0.01492 *
## educ 0.1309040 0.0253987 5.154 0.000000255045646 ***
## exper 0.1233472 0.0187587 6.575 0.000000000048500 ***
## expersq -0.0018871 0.0005999 -3.145 0.00166 **
## age -0.0528524 0.0084624 -6.246 0.000000000422204 ***
## kidslt6 -0.8683247 0.1183773 -7.335 0.000000000000221 ***
## kidsge6 0.0360056 0.0440303 0.818 0.41350
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1029.7 on 752 degrees of freedom
## Residual deviance: 802.6 on 745 degrees of freedom
## AIC: 818.6
##
## Number of Fisher Scoring iterations: 4Diferentemente do modelo linear, não podemos reportar diretamente os resultados da regressão como o valor dos estimadores, já que a relação entre eles é dada pela função de distribuição. Dessa forma, para obter interpretação dos parâmetros, precisamos obter os efeitos marginais, assim:
Probit: \(\frac{\delta E(Y|X)}{\delta X} = \Phi(X'\beta)*\beta\) , onde \(\Phi(z) = \frac{1}{\sqrt2\pi} exp(- \frac{1}{2}z^2)\),\(Z \sim N(0,1)\) em que \(\Phi(z)\) é a distribuição acumulada.
Logit: Seja \(\Delta(X'\beta)= \frac{exp(X'\beta)}{1 + exp(X'\beta)}\), temos \(\frac{\delta\Delta(X'\beta)}{\delta (X'\beta)} = \frac{d\Delta(X'\beta)}{d (X'\beta)} * \frac{d (X'\beta)}{dX}\)
No R:
library(mfx)
# Modelo Logit
logit.mfx <- logitmfx(inlf ~ nwifeinc + educ + exper + expersq + age +
kidslt6 + kidsge6,
data = mroz)
summary(logit.mfx)## Length Class Mode
## mfxest 28 -none- numeric
## fit 31 glm list
## dcvar 0 -none- character
## call 3 -none- call# Modelo Probit
probit.mfx <- probitmfx(inlf ~ nwifeinc + educ + exper + expersq + age +
kidslt6 + kidsge6,
data = mroz)
summary(probit.mfx)## Length Class Mode
## mfxest 28 -none- numeric
## fit 31 glm list
## dcvar 0 -none- character
## call 3 -none- callOs resultados dos efeitos marginais fornecem a mudança na probabilidade do evento de interesse (neste caso, participação na força de trabalho) para uma unidade de mudança nas variáveis explicativas.
A qualidade da previsão é dada por:
logit.fitted <- as.numeric(mlogit$fitted.values >= .5)
corr.pred.logit <- as.numeric(logit.fitted==mroz$inlf)
mean(corr.pred.logit)## [1] 0.7357238probit.fitted <- as.numeric(mprobit$fitted.values >= .5)
corr.pred.probit <- as.numeric(probit.fitted==mroz$inlf)
mean(corr.pred.probit)## [1] 0.7343958Uma vez que o \(R^2\) é uma medida de ajuste que leva em consideração uma relação linear entre a variável dependente e seus regressores, precisamos adotar uma nova medida que permita uma análise de ajustamento. O pseudo-R2 (McFadden) portanto calcula a razão entre a log-verossimilhança do modelo sem preditores, i.e, \(\beta_0\) e a log-verossimilhança do modelo completo. Assim:
\[pseudo-R2 = 1 - \frac{ln(g^{max})}{ln(g^{max}_0)}\] Em que \(g^{max}\) representa a probabilidade maximizada pelo modelo. Valores próximos a 0 indicam que o modelo escolhido não melhora a previsão em comparação com o modelo nulo, ou seja, os regressores adicionados não são bons preditores de Y.
Obtendo o pseudo-R2 no R:
logLik(mlogit)## 'log Lik.' -401.7652 (df=8)logLik(mprobit)## 'log Lik.' -401.3022 (df=8)require(mfx)
# Calculando a razão de chances (odds ratios)
logitor(inlf ~ nwifeinc + educ + exper + expersq + age +
kidslt6 + kidsge6,
data = mroz)## Call:
## logitor(formula = inlf ~ nwifeinc + educ + exper + expersq +
## age + kidslt6 + kidsge6, data = mroz)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## nwifeinc 0.9788810 0.0082435 -2.5346 0.011256 *
## educ 1.2475360 0.0541921 5.0915 0.000000355273436 ***
## exper 1.2285929 0.0393847 6.4220 0.000000000134459 ***
## expersq 0.9968509 0.0010129 -3.1041 0.001909 **
## age 0.9157386 0.0133450 -6.0403 0.000000001538446 ***
## kidslt6 0.2361344 0.0480729 -7.0898 0.000000000001343 ***
## kidsge6 1.0619557 0.0794229 0.8038 0.421539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretação:
OR = 1: Não há diferença nas chances de ocorrência de Y entre os grupos comparados.
OR > 1: As chances de ocorrência do evento são maiores no grupo de interesse em comparação ao grupo de referência. Por exemplo, se OR = 2, as chances de ocorrência de Y no grupo de interesse (ex. homens e mulheres) são duas vezes maiores do que no grupo de referência.
OR < 1: As chances de ocorrência do evento são menores no grupo de interesse em comparação ao grupo de referência. Por exemplo, se OR = 0.5, as chances de ocorrência de Y no grupo de interesse são metade das chances no grupo de referência.