Análise de Dados com o Software R:
Métodos Estatísticos, Computacionais e Econométricos

Prof. Adriano Azevedo Filho (azevedofilho@usp.br)

Análise de regressão II: introdução à regressão com resposta qualitativa

sumário geral | anterior | próximo

Conteúdo do Módulo

1 - Introdução a modelos de regressão com resposta qualitativa

2 - Estimativas nos modelos Logito e Probito: método da máxima verossimilhança

3 - Análise do conjunto de dados de inadimplência (dinad)

4 - Referências sugeridas

1 - Introdução a modelos de regressão com resposta qualitativa

Considere o conjunto de dados no arquivo dadosinad.csv descrito a seguir, assumindo que contem dados de operações financeiras de empréstimo realizada por um banco. Essas operações caracterizam o status do/a tomador/a do empréstimo como adimplentes (tomador/a do/a empréstimo honrou os pagamentos) ou inadimplentes (tomador do empréstimo não honrou os pagamentos). Cada linha do arquivo contém uma operação, com os dados da pessoa e da própria operação, assim como a região onde a operação foi realizada.

####
## Arquivo dadosinad.csv - dados de operações de empréstimo com
##                         status do/a tomador/a (adimplente/inadimplente)
####
## regiao   - região onde a operação foi realizada (1 ou 2)
## sexo     - sexo do/a tomador/a (M - masculino F - feminino) 
## idade    - idade em anos do/a tomador/a
## renda    - indicador de renda do/a tomador/a
## valorop  - valor da operação
## inad     - status do/a tomador/a (1 - inadimplente 0 - adimplente)
####

Leitura dos dados no data frame dinad

### Leitura do arquivo
dinad<-read.csv2("http://ihbs.com.br/html/dadosinad.csv")
head(dinad) ## listando as primeiras linhas do data frame
##   regiao sexo idade renda valorop inad
## 1      1    M  32.9  2190    4300    0
## 2      1    F  39.2  1482    2571    0
## 3      1    F  18.0   260     500    0
## 4      1    M  39.5  2656    3720    1
## 5      1    M  42.5  1607     500    0
## 6      1    M  48.9  3500    1254    1

Problema 1: como esse banco poderia utilizar esse conjunto de dados para entendimento das variáveis que podem “explicar” o fenômeno de inadimplência, visando o planejamento de ações que possam melhorar a gestão do processo de aprovação das operações?

Problema 2: desenvolva um modelo estatístico e estime os parâmetros desse modelo com os dados do data frame dinad, visando estimar/prever a probabilidade de inadimplência para uma nova operação potencial, em que são conhecidos os valores de sexo, regiao, valorop e renda.

1.1 Modelo com variável resposta qualitativa não pode ser uma regressão múltipla

No data frame dinad, a variável que queremos explicar, inad, que indica inadimplência (inad=1) e adimplência (inad=0) na operação de empréstimo, ainda que codificada numéricamente, é uma variável qualitativa. Se assumirmos que essa variável inad representa nosso \(y_i\), a variável resposta de nossa regressão, a qual gostariamos de “explicar” com o conhecimento do último módulo (regressão simples e múltipla), podemos ser tentados a considerar um modelo como

  • \(y_i=b_0+b_1 x_{1i}+b_2 x_{2i}+b_3 x_{3i}+b_4 x_{4i} +\epsilon_i\)
  • considerando as premissas usuais \(\epsilon_i\sim \text{N}(0,\sigma)\), \(\text{cov}(\epsilon_i,\epsilon_j)\) para \(i\neq j\) e outras premissas usuais quanto às variáveis explicativas.

Nesse modelo, \(y_i\) representaria inad e \(x_{1i}\), \(x_{2i}\), \(x_{3i}\), \(x_{4i}\), representariam, respectivamente, as variáveis sexo, valorop, renda e regiao, com as variáveis qualitativas (sexo e regiao) já recodificadas como variáveis “dummy” ou binárias (que correspondem a variáveis Bernoulli).

Esse modelo, na notação geral vetorial caracterizada no módulo anterior (veja tópico 9 do módulo anterior), poderia ser definido (considere \(m=4\)) por

  • \(\def\X{{\mathbf X}} \def\x{{\mathbf x}} \def\b{{\mathbf b}} \def\y{{\mathbf y}} \def\I{{\mathbf I}} \def\bfepsilon{{\boldsymbol{\epsilon}}} y_i=\b^T\x_i+\epsilon_i\) onde \(\b=\left(\begin{array}{c} b_0\\ b_1\\ \vdots\\ b_m\\ \end{array} \right)\; \text{e}\; \x_i=\left(\begin{array}{c} 1\\ x_{i1}\\ \vdots\\ x_{im}\\ \end{array} \right)\) para \(\;\;i=1,\ldots,n\)

Uma consequência importante desse modelo, com as premissas usualmente utilizadas, é

  • \(y_i|\x_i\sim \text{N}(\b^T\x_i,\sigma)\) com \(E(y_i|\x_i)=\b^t\x_i\) e \(V[y_i|\x_i]=\sigma^2\), com \(y_i|\x_i\), \(i=1,\ldots,n\), independentes.

Mas isso não é possível! nesse caso. \(y_i\) é uma variável qualitativa, que já foi condificada como variável Bernoulli (binária), só podendo assumir valores 1 ou 0. Isso exclui a possibilidade de ter uma distribuição Normal como premissa. O uso desse modelo, para essa situação, não é apropriado do ponto de vista técnico.

1.2 Solução através de modelos fundamentados na distribuição Bernoulli: Logito e Probito

No último tópico vimos que não é possível utilizar as premissas usuais da regressão linear dentro de uma situação em que a variável resposta é qualitativa, assumindo valores 1 e 0. A solução é considerar uma distribuição apropriada para essa situação: distribuição Bernoulli. Nessa formulação

  • \(y_i|\x_i\sim \text{Bernoulli}[p_i(\x_i)]\), com \(y_i|\x_i\) independentes, onde \(p_i(\x_i)\) é \(P(y_i|\x_i)=1\) essa notação deixa clara a dependência de \(p_i\) em \(\x_i\) (como se \(p_i\) fosse uma função de \(\x_i\)). E para toda a v.a. Bernoulli, nesse caso, temos que

  • \(E[y_i|\x_i]=p_i(\x_i)\)

Há várias possibilidades utilizadas para especificação de \(p_i(\x_i)\). As mais comuns são:

  • Modelo Logito: \(p_i(\x_i)=F_L(\b^T\x_i)=\displaystyle \frac{1}{1+e^{\displaystyle -\b^T\x_i}}\)

  • \(F_L(\cdot)\) é a chamada função logística avaliada no ponto \(\b^T\x_i\). Essa função também pode ser interpretada como a função cumulativa da distribuição de probabilidade logistica padronizada, com parâmetros \(\mu=0\) e \(s=1\) (a variância será \(\pi^2/3\)). Essa distribuição é parecida com a Normal padronizada, com um coeficiente de curtose um pouco maior.

  • Modelo Probito: \(p_i(\x_i)=\Phi(\b^T\x_i)\)

  • \(\Phi(\cdot)\) é a função cumulativa da Normal padronizada, avaliada no ponto \(\b^T\x_i\).

Os dois modelos tem similaridades. Eles operacionalizam um tipo de curva com formato de S que é chamada genericamente de sigmóide e é usado para caracterizar fenômenos em diversas disciplinas. Para comparação, a figura a seguir mostra as funções \(F_L(x)\) e \(\Phi(x)\) no domínio \(x\in [-4,\;4]\):

Como \(p_i(\x_i)\) é uma probabilidade, seu valor deve estar contido no intervalo \([0,\,1]\) e é possível demostrar que essas duas funções dão essa garantia.

O modelo Logito é mais popular que o Probito, talvez por ser mais tratável do ponto de vista algébrico. O modelo Probito depende de aproximações, visto que a função cumulativa da distribuição Normal não tem uma caracterização algébrica exata. Mas os programas estatísticos especializados, como o R e outros, já tem procedimentos específicos para lidar com essas dificuldades e qualquer um desses 2 modelos podem ser estimados sem grandes dificuldades, sendo potenciais candidatos para situações envolvendo variáveis qualitativas com 2 categorias.

2 - Estimativa dos modelos Logito e Probito: método da máxima verossimilhança

Considere o modelo Logito, por exemplo, onde temos

  • \(y_i|\x_i\sim \text{Bernoulli}[p_i(\x_i)]\) onde \(p_i(\x_i)=\displaystyle \frac{1}{1+e^{\displaystyle -\b^T\x_i}}\;\;\;\;\) para \(i=1,\ldots,n.\)

  • problema: como obter \(\hat \b\) ou seja, os estimadores dos parâmetros \(\hat b_0\), \(\hat b_1\), \(\ldots\), \(\hat b_m\) desse modelo (ou alternativamente do modelo Probito), a partir das observações \((y_i,\x_i)\), \(i=1,\ldots,n\)?

Ainda que seja possível a utilização do método dos quadrados mínimos, na sua versão não linear, para obtenção de estimativas de parâmetros de modelos desse tipo (veja aqui um exemplo disso), essa não é a prática usual.

2.1 Método da máxima verossimilhança

A estimativa dos parâmetros dos modelos Probito e Logito (e muitos outros), que são não-lineares, usa o chamado

  • método da máxima verossimilhança: nesse método, de forma simplificada, os estimativas serão obtidas pelos valores dos parâmetros que maximizam a probabilidade (ou densidade de probabilidade, quando for o caso) de observarmos \(y_i\) e \(\x_i\), \(i=1,\ldots,n\), (os valores observados no conjunto de dados) dentro do modelo utilizado.

O processo de otimização utilizado pelo método da máxima verossimilhança, em muitos casos, não leva, necessariamente, à “fórmulas” que estabelecem os valores dos estimadores (como na aplicação do método dos quadrados mínimos na regressão linear). Esse processo se vale de técnicas de otimização fundamentadas em métodos numéricos/computacionais iterativos, que tentam diretamente encontrar os valores dos parâmetros que maximizam a chamada função de verossimilhança, definida no contexto do método da máxima verossimilhança.

  • No caso dos modelos Logito e Probito, técnicas de otimização baseados no método de Newton-Raphson e suas variantes, são possíveis alternativas para obtenção dos parâmetros. Métodos como esses (e outros) estão implementados em funções de programas estatísticos com recursos para estimar modelos como o Logito e o Probito.

As inferências relacionadas ao método da maxima verossimilhança são aproximadas e dependem de resultados avançados, envolvendo normalidade, relacionados ao Teorema do Limite Central.

No R a principal função que estima parâmetros de modelos Logito e Probito (além de outros) chama-se:

  • glm: função do R utilizada para estimar parâmetros dos modelos Logito e Probito. Essa função é muito parecida com a função lm, em termos das definições das fórmulas e opções disponíveis.

Um detalhamento técnico mais profundo sobre método da máxima verossimilhança não está no escopo deste módulo. Faremos apenas uma breve exposição do método para efeito de uma familiarização introdutória sobre o assunto.

2.2 Aplicação do método da máxima verossimilhança a um problema mais simples

Neste tópico faremos uma breve introdução ao método da máxima verossimilhança num problema bem simples, a título de ilustração do método. Suponha que \(y_1\), \(y_2\), \(\ldots\), \(y_n\), são \(n\) resultados do lançamento de uma moeda, assumida como sendo uma amostra i.i.d. de uma v.a. Bernoulli com probabilidade \(p\) de resultado 1 (cara) e \(1-p\) de resultado 0 (coroa).

  • Recordando a distribuição Bernoulli: \(f(y)=p^y(1-p)^{1-y}\) para \(y=\{0,\,1\}\).

Para a aplicar o método precisamos construir a chamada função de verossimilhança, que é a caracterização da distribuição conjunta de probabilidade avaliada nos valores observados condicionada pelos parâmetros da distribuição. Nesse caso teríamos

  • \(f(y_1,y_2,\ldots,y_n|p)=f(y_1|p)f(y_2|p)\ldots f(y_n|p)\)

Se substituírmos \(f(y_i|p)\) pela sua definição para esse caso em particular, temos

  • \(f(y_1,y_2,\ldots,y_n|p)=p^{y_1}(1-p)^{1-y_1} p^{y_2}(1-p)^{1-y_2}\ldots p^{y_n}(1-p)^{1-y_n}\)

  • \(f(y_1,y_2,\ldots,y_n|p)=p^{\sum_{i=1}^n y_i}(1-p)^{n-\sum_{i=1}^n y_i}\)

Essa última expressão é exatamente a função de verossimilhança associada a esse problema, que é representada pela notação \(L(y_1,\ldots,y_n|p)\), \(L(p|y_1,\ldots,y_n)\) ou simplesmente \(L(p)\). Assim temos, para essa situação,

  • \(L(p)= \displaystyle p^{\sum_{i=1}^n y_i}(1-p)^{n-\sum_{i=1}^n y_i}\) (função de verossimilhança)

O problema agora é achar \(\hat p\), o valor de \(p\) que maximiza essa expressão. Esse valor será exatamente o estimador de máxima verossimilhança, ou seja,

  • \(\DeclareMathOperator*{\argmax}{argmax} \hat p=\argmax_{p\in [0,1]}\;\; p^{\sum_{i=1}^n y_i}\,(1-p)^{n-\sum_{i=1}^n y_i}\)

Esse problema pode ser solucionado facilmente pelo uso das condições de primeira ordem, derivando-se a expressão a ser maximizada com relação a \(p\) e igualando-se a zero, obtendo-se como solução:

  • \(\hat p= \frac{\sum_{i=1}^n y_i}{n}\) (estimador de máxima verossimilhança - problema simples)

Um problema de otimização como esse é usualmente mais fácil de resolver, do ponto de vista algébrico, pela aplicação do procedimento ao logarítmo da expressão a ser maximizada, em lugar da expressão original, algo que não irá alterar a solução do problema.

2.3 Aplicação do método da máxima verossimilhança a estimativa do modelo Logito

O tópico anterior é uma boa introdução ao método da máxima verossimilhança num problema muito similar ao que temos no Logito, também baseado na distribuição Bernoulli. Para caracterizar o problema para o caso do Logito, a unica diferença fundamental é que em lugar de \(y_i\) com probabilidade \(p\), para \(i=1,\ldots,n\), temos \(y_i|\x_i\) com probabilidade \(p_i (\x_i)\), \(i=1,\ldots,n\) com

  • \(p_i(\x_i)=\displaystyle \frac{1}{1+e^{\displaystyle -\b^T\x_i}}\;\;\;\;\) para \(i=1,\ldots,n.\)

O parâmetro a ser estimado não é mais \(p\) mas os parâmetros do vetor \(\b\), numa função de verossimilhança que será caracterizada por

  • \(L(\b)=\displaystyle \prod_{i=1}^n \left(\frac{1}{1+e^{\displaystyle -\b^T\x_i}}\;\;\;\right)^{y_i}\left(1-\frac{1}{1+e^{\displaystyle -\b^T\x_i}}\;\;\;\right)^{1-y_i}\)

e da mesma forma caracterizada anteriormente, os estimadores \(\hat \b\) serão os valores de \(\b\) que maximizam essa expressão, ou seja,

  • \(\DeclareMathOperator*{\argmax}{argmax} \hat \b=\argmax\;\; L(\b)\)

Ao contrário do problema anterior, em que chegamos a uma fórmula para o estimador de máxima verossimilhança, nesse caso o problema terá que ser solucionado por algorítmos iterativos de otimização, como o método de Newton-Raphson e outros, que são utilizados internamente pelas funções específicas que estimam parâmetros de modelos estatísticos. O modelo Probito, ainda que não explorado neste tópico, terá que ser solucionado por métodos similares. No R, a função glm já internaliza esses procedimentos de otimização para estimativa dos parâmetros de modelos desse tipo.

3 - Análise do conjunto de dados de inadimplência

Após uma breve introdução aos métodos de estimação fundamentados no método da máxima verossimilhança, no tópico anterior, vamos agora operacionalizar esses métodos no contexto do problema introduzido no tópico 1 deste módulo, envolvendo o problema de inadimplência em operações de empréstimo, visando responder as 2 perguntas formuladas naquele tópico. Essas perguntas visam esclarecer aspectos relacionados ao problema de inadimplência, dentro do contexto do conjunto de dados utilizado.

3.1 Uso da matriz de correlação para análise preliminar

O primeiro procedimento que utilizaremos é a obtenção da matriz de correlação com as variáveis do data frame dinad. Para colhermos alguma informação relativa ao efeito das variáveis qualitativas, criaremos versões dessas variaveis representadas por variáveis Bernoulli (ou dummy/binária). No caso da variável sexo criaremos a variável bsexo, com mulher representando o valor 1 e homem o valor 0. No caso da variável regiao, criaremos a variável bregiao, com valor 1 para região 2 e 0 para a região 1. Criaremos as 2 variáveis dentro do data frame dinad para facilitar.

dinad$bsexo=ifelse(dinad$sexo=="F",1,0)
dinad$bregiao=ifelse(dinad$regiao==2,1,0)
head(dinad)
##   regiao sexo idade renda valorop inad bsexo bregiao
## 1      1    M  32.9  2190    4300    0     0       0
## 2      1    F  39.2  1482    2571    0     1       0
## 3      1    F  18.0   260     500    0     1       0
## 4      1    M  39.5  2656    3720    1     0       0
## 5      1    M  42.5  1607     500    0     0       0
## 6      1    M  48.9  3500    1254    1     0       0

Para encontrar a matriz de correlação use

## a função with é usada para evitar o uso do nome completo das 
## variáveis (ex. dinad$inad, dinad$bsexo, etc.) restringindo a operação
## do data frame definido na opção data. A função cor usa a opção
## use="pairwise.complete.obs" que pode ser necessária se existirem
## observações com valores não definidos NA, NaN, etc. Veja documentação
## para detalhes.
with(data=dinad, cor(data.frame(inad,bsexo,bregiao,valorop,renda,idade),
                     use="pairwise.complete.obs"))
##                inad       bsexo      bregiao    valorop        renda
## inad     1.00000000 -0.34410251  0.079967416  0.3964914  0.382793389
## bsexo   -0.34410251  1.00000000  0.016548821 -0.3177317 -0.433086923
## bregiao  0.07996742  0.01654882  1.000000000  0.2109083 -0.005238992
## valorop  0.39649137 -0.31773168  0.210908326  1.0000000  0.772604494
## renda    0.38279339 -0.43308692 -0.005238992  0.7726045  1.000000000
## idade    0.14542470 -0.09812561  0.078477619  0.4697409  0.544071898
##               idade
## inad     0.14542470
## bsexo   -0.09812561
## bregiao  0.07847762
## valorop  0.46974086
## renda    0.54407190
## idade    1.00000000

A matriz de correlação mostra um “raio-x” dos dados, sugerindo relações importantes entre as variáveis envolvidas. Algumas regularidades importantes: (a) as variáveis bsexo, valorop e renda são as que apresentam maior correlação com inad, a primeira negativa e as duas outras positivas; (b) as mulheres tendem a ganhar menos que os homens e a utilizar operações menores (por que?); (c ) renda e valorop apresentam alta correlação positiva (0,77).

3.2 Estimativa do modelo completo usando a função glm

Inicialmente estimaremos os parâmetros do modelo Logito, definido nos tópicos anteriores, com \(y_i\) sendo representado pela variável inad e \(\x_i\) pelas variáveis bsexo, bregiao, valorop, renda e idade.

### Modelo Logito completo
mod1Logito<-glm(inad~bsexo+bregiao+valorop+renda+idade,data=dinad,
             family=binomial("logit"))
summary(mod1Logito)
## 
## Call:
## glm(formula = inad ~ bsexo + bregiao + valorop + renda + idade, 
##     family = binomial("logit"), data = dinad)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1505  -0.8388  -0.5368   0.9972   2.1665  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.0619279  0.4403916  -2.411 0.015895 *  
## bsexo       -1.0501541  0.2748068  -3.821 0.000133 ***
## bregiao      0.1883850  0.2595468   0.726 0.467947    
## valorop      0.0003375  0.0001148   2.940 0.003286 ** 
## renda        0.0003365  0.0002330   1.444 0.148703    
## idade       -0.0135086  0.0121407  -1.113 0.265852    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 487.63  on 359  degrees of freedom
## Residual deviance: 403.29  on 354  degrees of freedom
## AIC: 415.29
## 
## Number of Fisher Scoring iterations: 3

Os resultados desse tipo de estimação são relativamente diferentes dos apresentados no contexto da regressão linear múltipla. O R2 e o R2 ajustado como medidas de ajuste do modelo, não se aplicam nesse caso (há o conceito de pseudo-R2, mas não é de fato necessário). Uma possível medida de ajuste é o AIC (Akaike Information Criterion). De uma forma simplificada, quanto menor o valor do AIC melhor é o ajuste. As estimativas dos valores-p sugerem que somente o termo constante, e os parâmetros associados às variáveis bsexo e valorop são estatisticamente significativos.

Vamos estimar um novo modelo, somente considerando essas 2 variáveis explicativas bsexo e **valorop* :

### Modelo Logito completo
mod2Logito<-glm(inad~bsexo+valorop,data=dinad,
             family=binomial("logit"))
summary(mod2Logito)
## 
## Call:
## glm(formula = inad ~ bsexo + valorop, family = binomial("logit"), 
##     data = dinad)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1931  -0.8681  -0.5321   1.0401   2.0682  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.048e+00  2.597e-01  -4.036 5.44e-05 ***
## bsexo       -1.181e+00  2.572e-01  -4.592 4.38e-06 ***
## valorop      4.319e-04  7.671e-05   5.631 1.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 487.63  on 359  degrees of freedom
## Residual deviance: 405.91  on 357  degrees of freedom
## AIC: 411.91
## 
## Number of Fisher Scoring iterations: 3
coef(mod2Logito)
##   (Intercept)         bsexo       valorop 
## -1.0481020614 -1.1812829808  0.0004319309

O AIC é um pouco mais baixo, 411,9 que no caso anterior do modelo mod1Logito (AIC 415,3).

Vamos analisar em mais profundidade esse modelo (que eventualmente ainda poderia ser aprimorado). O que esse resultado esta dizendo é que

  • \(\hat E(y|\x)=\hat y|\x = \hat p(\x_i)=\displaystyle \frac{1}{1+e^{\displaystyle - (\hat b_0+\hat b_1 \text{bsexo}+\hat b_2 \text{valorop})}}\)

Na expressão não explicitamos o índice \(i\) de cada observação por se tratar de um modelo geral (válido para qualquer \(\x\)). No caso \(\hat b_0=-1{,}0481021\), \(\hat b_1=-1{,}1812830\) e \(\hat b_2=0{,}0004319\).

Para uma operação de uma pessoa do sexo feminino (bsexo é 1) e valorop 3000, a estimativa de \(\hat E[\text{inad}|\text{bsexo}=1,\text{valorop}=3000]\), que seria uma estimativa da probabilidade de inadimplência para essa situação, poderia ser obtida pela substituição desses valores na última expressão, obtendo-se

## Estimativa da probabilidade de inadimplência
## bsexo=1 (mulher) valorop=3000 (valor da operação)
1/(1+exp(-(-1.0481021-1.1812830*1+0.0004319* 3000)))
## [1] 0.2821777

A estimativa dessa probabilidade (condicional) seria 0,2822. Se essa operação, com \(\text{valorop}=3000\) fosse de um homem, teríamos \(\text{bsexo}=0\) e a estimativa seria 0,5616.

## Estimativa da probabilidade de inadimplência
## bsexo=0 (homem) valorop=3000 (valor da operação)
1/(1+exp(-(-1.0481021-1.1812830*0+0.0004319* 3000)))
## [1] 0.5615852

Esses últimos 2 resultados poderiam ter sido obtidos diretamente através da função predict, para as 2 situações que acabamos de examinar:

novobsexo=c(1,0) ## mulher e homem
novovalorop=c(3000,3000)
predict(mod2Logito,newdata=data.frame(bsexo=novobsexo,valorop=novovalorop),type="response")
##         1         2 
## 0.2821965 0.5616080

3.3 Apresentação gráfica dos resultados

É interessante a visualização dos resultados em um gráfico, o qual possibilitará uma melhor visualização dos efeitos das variáveis bsexo e valorop sobre a estimativa da probabilidade de inadimplência.

## diagrama de dispersão valorop vs. inad
with(data=dinad,plot(valorop,inad,col="blue",xlim=c(0,7000),ylim=c(0,1)),ylab="prob")
## resultados para as mulheres
title(expression(paste(hat(P),"(inad|sexo,valorop) - mod2Logito")))
novovalorop=c(0,7000)
novobsexo=rep(1,length(novovalorop))
novoinad<-predict(mod2Logito,newdata=data.frame(bsexo=novobsexo,valorop=novovalorop),type="response")
lines(novovalorop,novoinad,type="l",col="deeppink")
## resultados para os homens
novobsexo=rep(0,length(novovalorop))
novoinad<-predict(mod2Logito,newdata=data.frame(bsexo=novobsexo,valorop=novovalorop),type="response")
lines(novovalorop,novoinad,col="blue")

Esse resultado mostra, de um modo geral, que a estimativa da probabilidade de inadimplências tende a ser maior para os homens e cresce com o aumento no valor da operação (valorop). É importante destacar que isso é relativo à condição especifica dos dados, que envolve por razões de sigilo outros aspectos e não pode ser generalizado para outras situações.

3.4 Avaliação dos efeitos marginais

As estimativas do efeitos marginais dos parâmetros são usualmente definidas por:

  • \(\displaystyle \frac{\partial\, \hat E(y|\x)}{\partial\, x_j}\;\;=\;\;\frac{\partial\, \hat y|\x}{\partial\, x_j}\;\;\;\; j=1,\ldots,m\)

onde, no caso do Logito, temos

  • \(\hat E(y|\x)= \hat y|\x = \hat p(\x)=\displaystyle \frac{1}{1+e^{\displaystyle - (\hat b_0+\hat b_1 x_1 +\ldots+ \hat b_m x_m})}.\)

Num modelo linear essa estimativa corresponderá à estimativa \(\hat b_j\) associada à variável \(x_j\). Num modelo não linear como esse teríamos que encontrar a derivada da função, que pode ser uma função relativamente complicada. Usualmente, essa função depende dos valores das variáveis e uma possibilidade é encontrar o efeito marginal na média da variáveis explicativas (outra situações podem ocorrer).

A obtenção desses efeitos pode se fazer de várias maneiras no R. Uma das formas práticas é a através da função maBina do package erer que é especializado na análise de efeitos marginais, com muitas opções para a análise. Para uso desse package, reestimaremos o modelo mod2Logito com as opções x=TRUE e y=TRUE que são requeridas pelo função maBina (essas opções incluem os valores das variáveis dentro do modelo).

### Modelo Logito completo
require(erer)
mod2Logito<-glm(inad~bsexo+valorop,data=dinad,x=TRUE,y=TRUE,
             family=binomial("logit"))
maBina(w=mod2Logito,digits=7) 
##                 effect     error   t.value  p.value
## (Intercept) -0.2506853 0.0599147 -4.184038 3.61e-05
## bsexo       -0.2700663 0.0540925 -4.992673 9.00e-07
## valorop      0.0001033 0.0000186  5.559678 1.00e-07

O resultado mostra que o efeito do sexo na probabilidade é, aproximadamente, -0,27 (redução em passar de homem para mulher). O efeito marginal de valorop é 0,000103 para cada unidade de aumento dessa variável. Para um aumento de 1000 unidades esse valor seria aproximadamente 0,1 (em probabilidade de inadimplência).

3.5 Estimativa utilizando o modelo Probito

Só para completar esse tópico, faremos agora a estimativa dos parâmetros do modelo mod2Logito usando o modelo Probito e estimaremos os efeitos marginais, usando o mesmo procedimento do último tópico.

### Modelo Logito completo
mod2Probito<-glm(inad~bsexo+valorop,data=dinad,x=TRUE,y=TRUE,
             family=binomial("probit"))
summary(mod2Probito)
## 
## Call:
## glm(formula = inad ~ bsexo + valorop, family = binomial("probit"), 
##     data = dinad, x = TRUE, y = TRUE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2242  -0.8723  -0.5166   1.0513   2.1028  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.374e-01  1.569e-01  -4.063 4.84e-05 ***
## bsexo       -7.208e-01  1.528e-01  -4.716 2.40e-06 ***
## valorop      2.591e-04  4.458e-05   5.811 6.22e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 487.63  on 359  degrees of freedom
## Residual deviance: 405.36  on 357  degrees of freedom
## AIC: 411.36
## 
## Number of Fisher Scoring iterations: 4
maBina(w=mod2Probito,digits=7) 
##                 effect     error   t.value  p.value
## (Intercept) -0.2453583 0.0586160 -4.185862 3.58e-05
## bsexo       -0.2675590 0.0532374 -5.025774 8.00e-07
## valorop      0.0000997 0.0000173  5.763361 0.00e+00

O modelo Probito apresentou um desempenho um pouco melhor pelo critério do AIC (411,4 em lugar de 411,9 no modelo mod1Logito). As estimativas dos efeitos marginais associados às variáveis bsexo e valorop são muito similares, assim como a significância medida pelo valor-p.

3.6 Considerações finais

Soluções possíveis com relação aos problemas 1 e 2 colocados no tópico 1 deste módulo podem ser obtidas com facilidade a partir do desenvolvimento das idéias apresentadas neste tópico e nos tópicos anteriores, sugeridos como um ótimo exercício para os interessados.

Os modelos de resposta qualitativa não se restringem ao Logito e Probito. Há muitos outros. Uma área fértil de desenvolvimento chama-se GLM (modelos lineares generalizados), dentro da qual os modelos Logito e Probito são casos particulares de um modelo mais geral. Isso e outros desenvolvimentos mais avançados relativos a seleção de modelos e mesmo diagnóstico de problemas são assuntos mais avançados que não serão abordados neste módulo

4 - Referências sugeridas