Aula de Regressão Logística Simples no R
Aula de Regressão Logística Simples no R
Arthur Lerner
2023-06-22
Este primeiro exemplo tratará da regressão logística simples, portanto, utilizando somente uma variável independente, neste caso numérica.
Trata-se de uma amostra com 100 pessoas. A variável dependente é a ocorrência ou não (1 ou 0) de doença coronária cardíaca (CHD), associando-se com a idade (AGE) dos indivíduos, criando assim um modelo de regressão logística.
Limpar o R
rm(list=ls())Instalar os pacotes caso ainda não tenha instalado tirar o # da frente de install.packages para funcionar
#Instalando pacotes
#install.packages(c("ggplot2","readr","mfx","caret","pRoc", "ResourceSelection","modEvA","foreign","stargazer", "DataExplorer", "yardstick"))#Carregando pacotes exigidos:
library(readr)
library(dplyr)
library(ggplot2)
library(mfx)
library(caret)
library(pROC)
library(ResourceSelection)
library(modEvA)
library(foreign)
library(stargazer)
library(DataExplorer)
library(yardstick)
library(kableExtra)Ler os dados disponíveis na internet
chd <- read_delim("https://github.com/Smolski/livroavancado/raw/master/cdh.csv",
";", escape_double = FALSE, col_types = cols(CHD = col_factor(levels = c())),
trim_ws = TRUE)Visualizar os dados
kable(head(chd, 10)) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"))| AGE | AGRP | CHD |
|---|---|---|
| 20 | 1 | 0 |
| 23 | 1 | 0 |
| 24 | 1 | 0 |
| 25 | 1 | 0 |
| 25 | 1 | 1 |
| 26 | 1 | 0 |
| 26 | 1 | 0 |
| 28 | 1 | 0 |
| 28 | 1 | 0 |
| 29 | 1 | 0 |
summary(chd)## AGE AGRP CHD
## Min. :20.00 Min. :1.00 0:57
## 1st Qu.:34.75 1st Qu.:2.75 1:43
## Median :44.00 Median :4.00
## Mean :44.38 Mean :4.48
## 3rd Qu.:55.00 3rd Qu.:7.00
## Max. :69.00 Max. :8.00str(chd)## spc_tbl_ [100 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ AGE : num [1:100] 20 23 24 25 25 26 26 28 28 29 ...
## $ AGRP: num [1:100] 1 1 1 1 1 1 1 1 1 1 ...
## $ CHD : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 1 1 ...
## - attr(*, "spec")=
## .. cols(
## .. AGE = col_double(),
## .. AGRP = col_double(),
## .. CHD = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE)
## .. )
## - attr(*, "problems")=<externalptr>glimpse(chd)## Rows: 100
## Columns: 3
## $ AGE <dbl> 20, 23, 24, 25, 25, 26, 26, 28, 28, 29, 30, 30, 30, 30, 30, 30, 3…
## $ AGRP <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
## $ CHD <fct> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…Mostrar o gráfico de dispersão das observações Eixo x é idade e eixo y é possui ou não a doença CHD (1 possui e 0 não possui)
ggplot(chd, aes(x=AGE, y=CHD)) +
geom_point() +
stat_smooth(method="glm", method.args=list(family="binomial"), se=FALSE)
m1=glm(CHD~AGE, family = binomial(link="logit"), data = chd)
summary(m1)##
## Call:
## glm(formula = CHD ~ AGE, family = binomial(link = "logit"), data = chd)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9718 -0.8456 -0.4576 0.8253 2.2859
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.30945 1.13365 -4.683 2.82e-06 ***
## AGE 0.11092 0.02406 4.610 4.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 136.66 on 99 degrees of freedom
## Residual deviance: 107.35 on 98 degrees of freedom
## AIC: 111.35
##
## Number of Fisher Scoring iterations: 4Se observa o intercepto com o valor de -5,309, sendo que para a análise aqui proposta da relação entre CHD e AGE não obtém-se um significado prático para este resultado.No entanto, a variável de interesse é idade, que no modelo de regressão obteve o coeficiente de 0,1109. Pelo fato de ser positivo informa que quando a idade (AGE) se eleva, elevam-se as chances de ocorrência de CHD.De igual forma, nota-se que há significância estatística a p = 0,001 na utilização da variável AGE para o modelo, mostrando que possui importância ao modelo de regressão proposto.
Por fim, o modelo é utilizado para construção da predição de todos os valores das idades de todos os indivíduos desta amostra. Para isto, será criada um novo objeto contendo somente a variável dependente do modelo (AGE) e em seguida, é criada nova coluna constando os valores preditos. Assim, pode ser plotado um gráfico completo com todas as probabilidades desta base de dados:
# Filtrando a idade dos indivíduos
IDADE<-chd[,1]# Criando campo de predição para cada idade dos indivíduos
chd$PRED=predict(m1, newdata=IDADE, type="response")Plotando a probabilidade predita pelo modelo
ggplot(chd, aes(x=AGE, y=PRED)) +
geom_point()
Estimando a Razão de Chances
logitor(CHD~AGE,data = chd)## Call:
## logitor(formula = CHD ~ AGE, data = chd)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## AGE 1.117307 0.026882 4.6102 4.022e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O modelo de regressão logística, porém, traz os resultados dos estimadores na forma logarítma, ou seja, o log das chances da variável idade no modelo é 0,1109. No entanto, para uma interpretação mais enriquecida da relação da idade com o CHD é necessária a transformação deste coeficiente, ou seja, que seja efetuada a exponenciação da(s) variavel(eis) da regressão. Assim, obtém-se a razão das chances (OR - Odds Ratio em inglês) para as variáveis independentes.
Uma maneira prática de se obter a razão de chances no RStudio é utilizando o pacote mfx. Novamente o intercepto não nos interessa nesta análise mas sim a variável AGE. Como demonstrado abaixo, o resultado da razão de chances da variável AGE foi de 1,1173, o que pode assim ser interpretado: para cada variação unitária na idade (AGE), as chances de ocorrência de CHD aumentam 1,1173 vezes. Dito de outra forma, para cada variação unitária em AGE, aumentam-se 11,73% ((1,1173-1)*100) as chances da ocorrência de CHD.
Determinando o Intervalo de Confiança
exp(cbind(OR=coef(m1), confint(m1)))## OR 2.5 % 97.5 %
## (Intercept) 0.004944629 0.0004412621 0.0389236
## AGE 1.117306795 1.0692223156 1.1758681A determinação do intervalo de confiança do modelo proposto é relevante para que seja analizada a estimativa do intervalo de predição do coeficiente da variável independente, a um nível de confiança de 95%. Desta forma, em 95% dos casos, o parâmetro dos coeficientes estará dentro deste intervalo.
De forma prática é possível determinar os intervalos de confiança com o comando confint() commo observado abaixo, sendo que o coeficiente AGE toma o valor de 1,1173, podendo variar de 1,0692 a 1,1758.
Predição de Probabilidades
media = data.frame(AGE=mean(chd$AGE))
media## AGE
## 1 44.38A partir dos coeficientes do modelo de regressão logística é possível, portanto, efetuar a predição da variável categórica CHD, ou seja, saber a chance de ocorrer CHD com relação à uma determinada idade (AGE). No exemplo abaixo, primeiramente utilizamos a idade média das observações (44,38 anos), criando assim um novo data.frame chamado media.
media$pred.prob = predict(m1, newdata=media, type="response")
media## AGE pred.prob
## 1 44.38 0.4044944Para utilizar o valor da idade média na função de regressão obtida (m1), utiliza-se a função predict(), de acordo com valor da média encontrada (data.frame media).
O resultado mostra que para a idade média da amostra, 44,38 anos, há uma probabilidade de 40,44% na ocorrência da doença CHD. Esta ferramenta permite também a comparação pelo pesquisador das diferentes probabilidades entre as diversas idades (variável AGE).
Matriz de Confusão
chd$pdata <- as.factor(
ifelse(
predict(m1,
newdata = chd,
type = "response")
>0.5,"1","0"))
caret::confusionMatrix(chd$pdata, chd$CHD, positive="1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 45 14
## 1 12 29
##
## Accuracy : 0.74
## 95% CI : (0.6427, 0.8226)
## No Information Rate : 0.57
## P-Value [Acc > NIR] : 0.0003187
##
## Kappa : 0.4666
##
## Mcnemar's Test P-Value : 0.8445193
##
## Sensitivity : 0.6744
## Specificity : 0.7895
## Pos Pred Value : 0.7073
## Neg Pred Value : 0.7627
## Prevalence : 0.4300
## Detection Rate : 0.2900
## Detection Prevalence : 0.4100
## Balanced Accuracy : 0.7319
##
## 'Positive' Class : 1
## A matriz de confusão retoma uma excelente acurácia total do modelo em 74%, sendo que o modelo consegue acertos de 70,7% na predição de valores positivos ou dos “eventos” (29/41) e 76,3% na predição de valores negativos ou os “não eventos” (45/59).
Curva ROC
roc1=plot.roc(chd$CHD,fitted(m1))
plot(roc1,
print.auc=TRUE,
auc.polygon=TRUE,
grud=c(0.1,0.2),
grid.col=c("green","red"),
max.auc.polygon=TRUE,
auc.polygon.col="lightgreen",
print.thres=TRUE)
O teste Hosmer e Lemeshow
hl=hoslem.test(chd$CHD,fitted(m1),g=10)
hl##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: chd$CHD, fitted(m1)
## X-squared = 100, df = 8, p-value < 2.2e-16O teste de Hosmer e Lemeshow é utilizado para demonstrar a qualidade do ajuste do modelo, ou seja, se o modelo pode explicar os dados observados. Para este teste, os dados são divididos de acordo com as probabilidades previstas em 10 grupos iguais, sendo que os números previstos e os reais são comparados com a estatística do qui-quadrado. Hair et al. (2009) sugerem um tamanho de amostra de pelo menos 50 casos para a realização deste teste.
A hipótese nula H0 do qui-quadrado (p=0,05) deste teste é a de que as proporções observadas e esperadas são as mesmas ao longo da amostra. O modelo apresenta dificuldade de ajuste em função de que rejeita a hipótese nula a p=0,05.
Um grande valor de qui-quadrado (com pequeno valor de p < 0,05) indica um ajuste ruim e pequenos valores de qui-quadrado (com maior valor de p próximo a 1) indicam um bom ajuste do modelo de regressão logística.
Pseudo R2
RsqGLM(m1)
## $CoxSnell
## [1] 0.2540516
##
## $Nagelkerke
## [1] 0.3409928
##
## $McFadden
## [1] 0.2144684
##
## $Tjur
## [1] 0.2705749
##
## $sqPearson
## [1] 0.2725518Semelhante ao coeficiente de determinação R2 da regressão múltipla, a medida de pseudo R2 representam o ajuste geral do modelo proposto. Sua interpretação, portanto, é semelhante à regressão múltipla.
O modelo CoxSnell, por exemplo,explica 25% da variação na variável dependente. No entanto, é chamado de pseudo porque não é exatamente o R2 da regressão linear.
Comentários
Postar um comentário