Regressão Logística (opcional)

Às vezes, queremos prever uma variável dependente binária, ou seja, uma variável que pode assumir apenas dois valores, com base em várias variáveis independentes contínuas ou categóricas. Por exemplo, digamos que estamos interessados em testar se uma listagem é ou não uma jóia depende do preço e do tipo de quarto da listagem.

Não podemos usar ANOVA ou regressão linear aqui porque a variável dependente é uma variável binária e, portanto, normalmente não é distribuída. Outro problema com a ANOVA ou regressão linear é que ela pode prever valores que não são possíveis (por exemplo, nosso valor previsto pode ser 5, mas apenas 0 e 1 fazem sentido para essa variável dependente). Portanto, usaremos regressão logística. A regressão logística primeiro transforma a variável dependente Y com a transformação do logit. A transformação do logit usa o logaritmo natural das chances de que a variável dependente seja igual a 1:

\[ probabilidades=\displaystyle\frac{P(Y=1)}{P(Y=0}=\displaystyle\frac{P(Y=1)}{1-P(Y=1)} \]

então o logit \(P(Y=1))=\ln\displaystyle\frac{P(Y-1)}{1-P(Y=1)}\)

Isso garante que nossa variável dependente seja normalmente distribuída e não seja restrita a ser 0 ou 1.Isso garante que nossa variável dependente seja normalmente distribuída e não seja restrita a ser 0 ou 1.

Vamos realizar a regressão logística:

logistic.model <- glm(gem ~ price * room_type, data=airbnb, family="binomial") # o argumento family = "binomial" diz R para tratar a variavel dependente como uma variavel 0/1
summary(logistic.model) # saida da regressao
## 
## Call:
## glm(formula = gem ~ price * room_type, family = "binomial", data = airbnb)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4909  -0.3819  -0.3756  -0.3660   2.5307  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -2.5270702  0.0570207 -44.318   <2e-16 ***
## price                       -0.0007185  0.0004030  -1.783   0.0746 .  
## room_typePrivate room       -0.0334362  0.1072091  -0.312   0.7551    
## room_typeShared room         1.0986318  1.1278159   0.974   0.3300    
## price:room_typePrivate room -0.0011336  0.0012941  -0.876   0.3810    
## price:room_typeShared room  -0.0562610  0.0381846  -1.473   0.1406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 8663.8  on 17650  degrees of freedom
## Residual deviance: 8648.6  on 17645  degrees of freedom
## AIC: 8660.6
## 
## Number of Fisher Scoring iterations: 7

Vemos que o único preditor marginalmente significativo de uma listagem ser ou não uma jóia é o preço da listagem. Você pode relatar o seguinte: “Controlando o tipo de quarto e a interação entre preço e tipo de quarto, havia uma relação negativa marginalmente significativa entre o preço e a probabilidade de uma listagem ser uma jóia (\(\beta\) = -0.0007185, \(\chi\) (17645) = -1,783, p = 0,075).

A interpretação dos coeficientes de regressão na regressão logística é diferente da do caso da regressão linear:

summary(logistic.model)
## 
## Call:
## glm(formula = gem ~ price * room_type, family = "binomial", data = airbnb)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4909  -0.3819  -0.3756  -0.3660   2.5307  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -2.5270702  0.0570207 -44.318   <2e-16 ***
## price                       -0.0007185  0.0004030  -1.783   0.0746 .  
## room_typePrivate room       -0.0334362  0.1072091  -0.312   0.7551    
## room_typeShared room         1.0986318  1.1278159   0.974   0.3300    
## price:room_typePrivate room -0.0011336  0.0012941  -0.876   0.3810    
## price:room_typeShared room  -0.0562610  0.0381846  -1.473   0.1406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 8663.8  on 17650  degrees of freedom
## Residual deviance: 8648.6  on 17645  degrees of freedom
## AIC: 8660.6
## 
## Number of Fisher Scoring iterations: 7

O coeficiente de regressão do preço é -0.0007185. Isso significa que um aumento de uma unidade no preço levará a um aumento de -0.0007185 nas chances de log de joia ser igual a 1 (ou seja, de uma listagem sendo uma joia). Por exemplo:

\[ \textrm{logit}(P(Y = 1 | \textrm{price} = 60 )) = \textrm{logit}(P(Y = 1 | \textrm{price} = 59 )) - 0.0007185 \]

\[ \mbox{logit}(P(Y = 1 | \mbox{price} = 60 )) = \mbox{logit}(P(Y = 1 | \mbox{price} = 59 )) - 0.0007185 \] \[ \Leftrightarrow \textrm{logit}(P(Y = 1 | \textrm{price} = 60 )) - \textrm{logit}(P(Y = 1 | \textrm{price} = 59 )) = -0.0007185 \] \[ \Leftrightarrow \textrm{ln(probs}(P(Y = 1 | \textrm{price} = 60 )) - \textrm{ln(probs}(P(Y = 1 | \textrm{price} = 59 )) = -0.0007185 \] \[ \Leftrightarrow \textrm{ln}(\dfrac{ \textrm{probs}(P(Y = 1 | \textrm{price} = 60 )}{\textrm{probs}(P(Y = 1 | \textrm{price} = 59 )}) = -0.0007185 \]

\[ \Leftrightarrow \dfrac{ \textrm{probs}(P(Y = 1 | \textrm{price} = 60 )}{\textrm{probs}(P(Y = 1 | \textrm{price} = 59 )} = e^{-0.0007185} \] \[ \Leftrightarrow \textrm{probs}(P(Y = 1 | \textrm{price} = 60 ) = e^{-0.0007185} * \textrm{probs}(P(Y = 1 | \textrm{price} = 59 ) \]

Assim, o coeficiente de regressão em uma regressão logística deve ser interpretado como o aumento relativo nas chances da variável dependente ser igual a 1, para cada aumento de unidade no preditor, controlando todos os outros fatores em nosso modelo. Nesse caso, as probabilidades de uma listagem ser uma gema devem ser multiplicadas por $ e ^ {0,0007185} = 0,999 $ ou diminuídas em 0,1%, para cada aumento de preço unitário. Em outras palavras, listagens mais caras têm menos probabilidade de serem jóias. No exemplo específico acima, as chances de ser uma joia de uma listagem com preço de 60 são \(e^{(-0.0007185×5)}= 0,996\) vezes a chance de ser uma jóia de uma listagem com preço de 59.

Medindo o ajuste de uma regressão logística: porcentagem classificada corretamente

Nosso modelo usa o preço e o tipo de quarto da listagem para prever se a listagem é uma joia ou não. Ao fazer previsões, é natural nos perguntarmos se nossas previsões são boas. Em outras palavras, usando preço e tipo de quarto, com que frequência prevemos corretamente se uma listagem é uma jóia ou não? Para ter uma idéia da qualidade de nossas previsões, podemos pegar o preço e o tipo de quarto das listagens em nosso conjunto de dados, prever se as listagens são gemas e comparar nossas previsões com o status real da gema das listagens. Vamos primeiro fazer as previsões:

airbnb <- airbnb %>% 
  mutate(prediction = predict(logistic.model, airbnb))
# Crie uma nova coluna chamada previsao no quadro de dados do airbnb e armazene nela a previsao,
# baseado em logistic.model, para os dados do airbnb
# De uma olhada nessas previsoes:
head(airbnb$prediction)
##         1         2         3         4         5         6 
## -4.790232 -4.448355 -4.049498 -4.619293 -4.106477 -4.847211
# Compare com as observacoes:
head(airbnb$gem)
## [1] no gem no gem no gem no gem no gem no gem
## Levels: no gem gem

Você vê o problema? As previsões são logits, ou seja, logaritmos de chances de que as listagens sejam gemas, mas as observações simplesmente nos dizem se uma listagem é uma gema ou não. Para uma comparação significativa entre previsões e observações, precisamos transformar os logits em uma decisão: gema ou não gema. É fácil transformar logits em probabilidades usando o exponencial do logit. A relação entre probabilidades e probabilidades é a seguinte:

\[ \textrm{probs} = \dfrac{P(Y = 1)}{P(Y = 0)} = \dfrac{P(Y = 1)}{1 - P(Y = 1)} \] \[ \Leftrightarrow \dfrac{probs}{1 + \dfrac{P(Y = 1)}{1 - P(Y = 1)}} = \dfrac{\dfrac{P(Y = 1)}{1 - P(Y = 1)}}{1 + \dfrac{P(Y = 1)}{1 - P(Y = 1)}} \]

\[ \Leftrightarrow \dfrac{probs}{1 + probs} = \dfrac{\dfrac{P(Y = 1)}{1 - P(Y = 1)}}{\dfrac{1 - P(Y = 1) + P(Y = 1)}{1 - P(Y = 1)}} \]

\[ \Leftrightarrow \dfrac{probs}{1 + probs} = \dfrac{P(Y = 1)}{1 - P(Y = 1) + P(Y = 1)} \] \[ \Leftrightarrow \dfrac{probs}{1 + probs} = P(Y = 1) \] Agora vamos calcular, para cada listagem, a probabilidade de a listagem ser uma jóia (gem):

airbnb <- airbnb %>% 
  mutate(prediction.logit = predict(logistic.model, airbnb),
         prediction.odds = exp(prediction.logit),
         prediction.probability = prediction.odds / (1+prediction.odds))

# Inspecionando as probabilidades preditas
head(airbnb$prediction.probability)
##           1           2           3           4           5           6 
## 0.008242034 0.011562542 0.017132489 0.009763496 0.016198950 0.007789092

Os primeiros números são probabilidades muito baixas, mas também existem probabilidades mais altas e todas as previsões estão entre 0 e 1, como deveriam. Agora precisamos decidir qual probabilidade é suficiente para prevermos que uma listagem é uma gem ou não.

Uma escolha óbvia é uma probabilidade de 0,5: uma probabilidade maior que 0,5 significa que prevemos que é mais provável que uma listagem seja uma gem do que não. Vamos converter probabilidades em previsões:

airbnb <- airbnb %>% 
  mutate(prediction = case_when(prediction.probability<=.50 ~ "no gem",
                                prediction.probability> .50 ~ "gem"))

# Inspecinando as predicoes
head(airbnb$prediction)
## [1] "no gem" "no gem" "no gem" "no gem" "no gem" "no gem"

Uma etapa final é comparar previsões com observações:

table(airbnb$prediction, airbnb$gem)
##         
##          no gem   gem
##   no gem  16471  1180

Normalmente, vemos uma tabela 2x2, mas vemos uma tabela com um valor previsto nas linhas e dois valores observados nas colunas. Isso ocorre porque todas as probabilidades previstas estão abaixo de 0,50 e, portanto, sempre previmos que não há gemas. Vamos reduzir o limite para prever que uma listagem é uma gem:

airbnb <- airbnb %>% 
  mutate(prediction = case_when(prediction.probability<=.07 ~ "no gem",
                                prediction.probability> .07 ~ "gem"),
         prediction = factor(prediction, levels = c("no gem","gem")))# verifique se nenhuma joia eh o primeiro nivel do nosso fator

# Observando a tabela
table(airbnb$prediction,airbnb$gem)
##         
##          no gem   gem
##   no gem  11240   854
##   gem      5231   326

Podemos ver que, com esta regra de decisão (prever gems sempre que a probabilidade prevista de gems for superior a 7%), obtivemos 11240 + 326 corretas e 5231 + 854 previsões erradas, que é uma taxa de acerto de (11240 + 326) / (11240 + 326 + 5231 + 854) = 65,5%.