Homens têm mais casos extraconjugais?

Probits

Probits são, essencialmente, modelos lineares generalizados (GLM) em que a variável de resposta assume valores binários. Os parâmetros estimados no dizem em que medida uma variação marginal em uma variável explicativa altera o \(z\)-score da variável dependente, não a sua probabilidade condicional de assumir \(1\).

Seja \(\mathbb{P}\) o operador para probabilidade, \(\Phi\) a função de distribuição cumultiva de uma normal padrão, \(X\) um vetor de variáveis aleatórias e \(\beta\) o vetor de parâmetros a ser estimado. Um modelo probit assume a seguinte forma:

\[\mathbb{P}(Y=1 | X) = \Phi (X^T \beta)\]

É possível também exprimir o modelo como \(Y^{*} = X^T \beta + \epsilon\) e mostrar que as duas formas são equivalentes é um exercício interessante.

Não irei entrar nos pormenores de como se pode estimar os parâmetros desse modelo. A função glm disponível no R usa estimação por Máxima Verossimilhança, até onde sei.

A amostra

Vamos explorar um pouco os dados com a ajuda do ggplot2 antes de sair por aí estimando parâmetros:

library(wooldridge)

data("affairs")

head(affairs) ## Confirmando que funcionou

##   id male age yrsmarr kids relig educ occup ratemarr naffairs affair
## 1  4    1  37    10.0    0     3   18     7        4        0      0
## 2  5    0  27     4.0    0     4   14     6        4        0      0
## 3  6    1  27     1.5    0     3   18     4        4        3      1
## 4 11    0  32    15.0    1     1   12     1        4        0      0
## 5 12    0  27     4.0    1     3   17     1        5        3      1
## 6 16    1  57    15.0    1     5   18     6        5        0      0
##   vryhap hapavg avgmarr unhap vryrel smerel slghtrel notrel
## 1      0      1       0     0      0      0        1      0
## 2      0      1       0     0      0      1        0      0
## 3      0      1       0     0      0      0        1      0
## 4      0      1       0     0      0      0        0      0
## 5      1      0       0     0      0      0        1      0
## 6      1      0       0     0      1      0        0      0

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(scales)

affairs %>%
  ggplot(aes(x=yrsmarr))+
  geom_histogram(binwidth = 1, fill = "white", color="blue")+
  xlab("Anos de Casamento")+
  ylab("Frequência")

affairs %>%
  ggplot(aes(x=ratemarr))+
  geom_histogram(binwidth = 1, fill = "white", color="blue")+
  xlab("Avaliação do Casamento (1 = infeliz, 5 = muito feliz)")+
  ylab("Frequência")

affairs %>%
  ggplot(aes(x=affair))+
  geom_histogram(aes(y=..density..), bins = 2,
                  fill = "white", color="blue")+
  scale_y_continuous(labels = percent)+
  xlab("Teve um caso no ano anterior (1 = sim, 0 = não)")+
  ylab("")

Estimando um modelo probit

O R “vanilla” já conta com ferramentas para estimar modelos desse tipo:

affairs$rel <- ifelse(affairs$relig > 2, 1, 0) #dummy de religiosidade
affairs$feliz <- ifelse(affairs$ratemarr > 2, 1, 0) #dummy de casamento feliz

probit <- glm(affair ~ male + kids + feliz + rel + educ + yrsmarr + age,
                family = binomial(link = "probit"),
                  data = affairs)
summary(probit)

## 
## Call:
## glm(formula = affair ~ male + kids + feliz + rel + educ + yrsmarr + 
##     age, family = binomial(link = "probit"), data = affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5816  -0.7550  -0.6043  -0.3185   2.3605  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.01521    0.47890   0.032  0.97466    
## male         0.20968    0.12964   1.617  0.10579    
## kids         0.23355    0.16199   1.442  0.14937    
## feliz       -0.74577    0.15740  -4.738 2.16e-06 ***
## rel         -0.25022    0.12236  -2.045  0.04086 *  
## educ         0.01169    0.02620   0.446  0.65554    
## yrsmarr      0.05518    0.01857   2.971  0.00297 ** 
## age         -0.02583    0.01036  -2.494  0.01264 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 625.27  on 593  degrees of freedom
## AIC: 641.27
## 
## Number of Fisher Scoring iterations: 4

O problema desse modelo no atual estado é que os parâmetros representam mudanças no \(z\)-score da variável explicada e não são facilmente interpretados. Seria mais interessante ter variações diretas na probabilidade condicional \(\mathbb{P}(Y=1 | X)\). Temos ferramentas para isso, no pacote mfx.

library(mfx)

probitmfx(affair ~ male + kids + feliz + rel + educ + yrsmarr + age,
              data = affairs)

## Call:
## probitmfx(formula = affair ~ male + kids + feliz + rel + educ + 
##     yrsmarr + age, data = affairs)
## 
## Marginal Effects:
##              dF/dx  Std. Err.       z     P>|z|    
## male     0.0645982  0.0399656  1.6163  0.106020    
## kids     0.0690378  0.0458475  1.5058  0.132115    
## feliz   -0.2637864  0.0600435 -4.3933 1.117e-05 ***
## rel     -0.0787878  0.0393153 -2.0040  0.045070 *  
## educ     0.0035905  0.0080486  0.4461  0.655527    
## yrsmarr  0.0169520  0.0056885  2.9800  0.002882 ** 
## age     -0.0079351  0.0031755 -2.4988  0.012461 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "male"  "kids"  "feliz" "rel"

Aqui sim, os parâmetros podem ser lidos como variações na probabilidade condicional de \(Y=1\). Observe, leitor que a variável explicativa male não é estatisticamente significante. Talvez precisamos de uma amostra maior.