Começando com um preditor
Para entendermos a construção da regressão logística, vamos começar com um banco de dados sobre o desfecho dibatetes. Vamos avaliar se os preditores glicose (glucose) e idade (age) são bons preditores para esse desfecho. Abaixo segue a linha de código da recuperação do banco de dados.
O banco de dados contem 768 observações. Figura 1 mostra a relação entre glicose e diagnóstico de diabetes.
Código
PimaIndiansDiabetes %>%
ggplot(aes(x = glucose, y = as.numeric(diabetes == "pos") * 0.5, color = diabetes)) +
geom_jitter(height = 0.2, width = 0.2, alpha = 0.3, size = 4) +
scale_color_manual(values = c("neg" = "tomato3", "pos" = "red")) +
labs(
title = "Nivel de glicose e diagnostico",
x = "Nivel de glicose"
) +
scale_y_continuous(breaks = c(0, 0.5), labels = c("No Diabetes", "Diabetes")) +
theme_minimal() +
theme(
legend.position = "none",
axis.line.x = element_line(color = "black"),
plot.title = element_text(hjust = 0.5),
axis.title.y = element_blank()
)
Portanto, vamos criar um modelo que nos permite predizer a probabilidade de ser diagnosticado com diabetes baseado na glicemia e idade.
Código
library(plotly)
model <- glm(diabetes ~ glucose, data = PimaIndiansDiabetes, family = binomial)
# Generate data for the curve
glucose_range <- data.frame(glucose = seq(min(PimaIndiansDiabetes$glucose, na.rm = TRUE),
max(PimaIndiansDiabetes$glucose, na.rm = TRUE),
length.out = 500))
glucose_range$prob <- predict(model, newdata = glucose_range, type = "response")
# Create interactive plot
plot_ly(data = glucose_range,
x = ~glucose,
y = ~prob,
type = 'scatter',
mode = 'lines',
name = "Probability Curve") %>%
layout(
title = "Probability of Diabetes by Glucose Level",
xaxis = list(title = "Glucose Level"),
yaxis = list(title = "Probability of Diabetes")
)Observe na Figura 2 a curva sigmóide (em forma de S) clássica da regressão logística.
Criando um modelo simples
Primeiro, vamos criar um modelo simples.
Código
Call:
glm(formula = diabetes ~ glucose, family = binomial, data = PimaIndiansDiabetes)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.350080 0.420827 -12.71 <2e-16 ***
glucose 0.037873 0.003252 11.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 993.48 on 767 degrees of freedom
Residual deviance: 808.72 on 766 degrees of freedom
AIC: 812.72
Number of Fisher Scoring iterations: 4Observando os resultados da regressão logística, vamos nos concentrar em apenas alguns números (por enquanto). Primeiramente, temos os coeficientes (no resultado chamado estimate). O que são eles? Bem, temos um modelo matemático que possui “espaços reservados” onde valores específicos podem ser usados para que o modelo descreva melhor a relação do preditor, no caso glicose, e o desfecho (diagnóstico de diabetes). Você notará que, neste caso, temos dois coeficientes:
Intercepto (ou Coeficiente de Interceptação): Ele nos diz como o modelo se comporta quando o nível de glicose é zero (que é onde o modelo intercepta o eixo y). Neste exemplo, esse é um número sem significado prático, porque representaria uma pessoa com o nível de glicose zero. E nesse caso, tal valor não é possível. Portanto, podemos ignorar esse número.
Coeficiente para glicose: é o coeficiente da nossa variável preditora, a glicose. Ele nos dirá (a partir dos dados que temos) o que acontece quando o valor da glicose muda. Especificamente, como isso afeta as chances de uma pessoa ser diagnosticada com diabetes. A mudança pode ser positiva (como é o caso neste exemplo), negativa ou zero. Um número positivo significa que, à medida que os níveis de glicose aumentam, as chances de ser diagnosticado com diabetes também aumentam. Se fosse zero, não haveria mudança (a glicose não teria impacto). Se fosse negativo, significaria que, à medida que a glicose aumenta, as chances de ser diagnosticado com diabetes diminuiriam.
Você notará que usei um termo genérico e um tanto vago, “chances” (em inglês, odds), em vez de “probabilidade”. Isso ocorre porque o número real que o modelo produz é o que chamamos de “log-odds” (ou o logaritmo da razão de chances). Não se preocupe com isso por enquanto. Tudo o que você precisa focar é no fato de que o coeficiente é positivo. Você também notará que os resultados incluem um valor-p (p value), que neste caso é muito pequeno. Isso significa que este resultado é estatisticamente significativo, o que significa que glicose é um preditor importante para o fato de ser ou não diagnosticado com diabetes.
Analisaremos outros valores fornecidos nos resultados à medida que construirmos o modelo e adicionarmos variáveis adicionais.
Bem, agora vamos ver como é a relação do diagnóstico de diabetes com idade no gráfico abaixo.
Código
PimaIndiansDiabetes %>%
ggplot(aes(x = age, y = as.numeric(diabetes == "pos") * 0.5, color = diabetes)) +
geom_jitter(height = 0.2, width = 0.2, alpha = 0.3, size = 4) +
scale_color_manual(values = c("neg" = "darkblue", "pos" = "red")) +
labs(
title = "Idade e diagnostico",
x = "Idade (anos)"
) +
scale_y_continuous(breaks = c(0, 0.5), labels = c("No Diabetes", "Diabetes")) +
theme_minimal() +
theme(
legend.position = "none",
axis.line.x = element_line(color = "black"),
plot.title = element_text(hjust = 0.5),
axis.title.y = element_blank()
)
A Figura 3 não demonstra uma relação clara entre o diagnóstico e idade. Entretanto, isso não descarta a possibilidade e importância. Portanto, vamos incluir idade no modelo junto com glicose para analisar.
Código
Call:
glm(formula = diabetes ~ glucose + age, family = binomial, data = PimaIndiansDiabetes)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.912449 0.462620 -12.78 < 2e-16 ***
glucose 0.035644 0.003290 10.83 < 2e-16 ***
age 0.024778 0.007374 3.36 0.000778 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 993.48 on 767 degrees of freedom
Residual deviance: 797.36 on 765 degrees of freedom
AIC: 803.36
Number of Fisher Scoring iterations: 4Podemos observar que agora temos três coeficientes (estimate), um para o intercepto que significa a probabilidade se os preditores são zero. No caso, glicose e idade não serão considerados pela falta de sentido nos valores nulos. Temos glicose e idade, cada um com seu coeficiente. O coeficiente de glicose foi de 0.04 . Observem que esse número é positivo, indicando que as chances de ser diagnósticado com diabetes aumentam a medida que a glicose aumenta. O coeficiente de idade foi de 0.02`. O coeficiente para idade também foi positivo indicando chance maior de ser diagnosticado com diabetes se a idade aumenta. Vejam que o valor de p é menor do que 0.05, o que significa que os preditores são importantes para o desfecho diabetes.
Avaliando os pressupostos do modelo
Os pressupostos são condições em que a regressão logística binária espere dos dados observados. Ao longo da demonstração, vou inserir uma nota demonstrando, quando possível, uma forma mais rápida de demonstrar o pressuposto.
Temos os seguintes pressupostos:
Tamanho amostral adequado
Não é um pressuposto formal, mas é importante. Regressão Logística, por ser baseada na função de Máxima Verossimilhança, exige um número suficiente de eventos no grupo menos frequente. Como regra prática temos:
É recomendado ter pelo menos 10 eventos por cada variável preditora incluída no modelo. Exemplo: se no seu estudo sobre o desfecho de diagnóstico de diabete, se você tivesse 50 participantes no estudo, você poderia adicionar no máximo 5 variáveis preditoras (50/10).
Por que é importante: amostras pequenas (ou número insuficiente de eventos) levam a erros-padrão grandes e coeficientes instáveis.
Relação linear na escala log
A regressão logística binária espera que exista uma relação linear entre variáveis preditoras contínuas e o logarítmo da razão das chances.
Para provar essa relação linear, uma das formas (aqui apresento uma forma visual) é plotar os resíduos do desvio (residuais deviance) em relação à variável contínua (no nosso caso age). Para isso, abaixo, vou criar um vetor chamado residuos_desvio e plotar contra nossa variável contínua (age).
Código
Portanto, observe na Figura 3 que os resíduos estão uniformemente distribuídos ao longo de uma linha imaginária horizontal no valor zero do eixo y. Se houvesse uma forma de U ou S, teríamos que corrigir a linearidade, fazendo uma tranformação da variável age em log(age) ou raiz quadrada(age) ou 1/(age).
Independência das observações
As observações (linhas de dados) devem ser independentes umas das outras. O resultado de uma observação não deve influenciar o resultado de outra. Isto é importante porque a função da verossimilhança utilizada para estimar os coeficientes assume que os erros (resíduos) de cada observação são indendentes.
Se os dados são agrupados ou repetidos no tempo (dados longitudinais), este pressuposto é violado, exigindo o uso de modelos mais complexos, como Modelos Mistos Generalizados (GLMM).
Para verificar a independência dos resíduos, podemos utilizar o teste Durbin Watson. Se o valor desse teste estiver entre 1.5 e 2.5, isso sugere que os resíduos são independentes.
lag Autocorrelation D-W Statistic p-value
1 0.01382267 1.970949 0.686
Alternative hypothesis: rho != 0Observe no teste que o valor D-W Statistic foi de 1.97, estando dentro da faixa de 1.5 - 2.5, indicando ausência de multicolinearidade.
Se o teste de Durbin-Watson estiver fora da faixa considerada e o valor de p for significativo, você tem um problema de dependência. Nesse caso pode-se:
- Incluir efeito temporal: se a dependência for devido ao tempo, inclua a variável tempo/ordem no modelo como um preditor.
2. Usar modelos avançados: se a dependência for estrutural (dados hierárquicos, clusters ou medições repetidas no mesmo indivíduo), você deve recorrer a Modelos de Equações de Estimação Generalizada (GEE) ou Modelos Mistos Generalizados (GLMM), que são projetados especificamente para modelar a estrutura de dependência.
É possível analisar a independência dos resíduos realizando uma inspeção visual. Para isso, é necessário calcularmos os resíduos dos desvios, já calculado para gerar a Figura 4. E precisamos calcular os valores preditos. Abaixo vamos fazer isso juntos e plotar o gráfico.
Código
Veja que na Figura 5 não observamos um padrão na forma de funil, S ou U, indicando que podemos considerar a independência dos resíduos.
O R possui um pacote chamado performance e dentro dele temos a função check_autocorrelation. Essa função já faz o cálculo Durbin Watson e te entrega uma mensagem como essa. Para isso, é só executar:
Ausência de Multicolinearidade
A multicolinearidade só é avaliada em modelos multilinears (com mais de uma variável preditora). As variáveis preditoras (as variáveis independentes, ou X’s) não devem ser altamente correlacionadas entre si.
Isto é importante porque alta multicolinearidade torna os erros-padrão dos coeficientes instáveis e inflacionados, dificultando a interpretação e superestimando o valor de p, podendo levar a um erro do tipo II.
Para verificar a multicolinearidade, podemos calcular o Fator de Inflação da Variância (VIF). Geralmente, valores de VIF acima de 5 ou 10 indicam problemas.
Abaixo, vamos verificar o VIF do nosso modelo1, utilizando a função vif() do pacote car.
Podemos observar que o valor de VIF para glicose 1.02 e de idade 1.02 são menores do que 10, indicando ausência de multicolinearidade.
Se existir multicolinearidade, não há muito o que fazer. O ideal é retirar uma das variáveis altamente correlacionadas e avaliar o modelo. Não existe uma forma de escolher qual variável remover. É necessário você remover e avaliar qual delas que, ao permanecer no modelo, o deixa melhor.
Ausência de valores influentes e outliers
Valores influentes são aqueles que, se removidos do modelo, influenciam os valores dos coeficientes dos preditores. Outliers são valores que destoam do restante.
Para verificar valores influentes, vamos calcular:
Distância de Cook: o quanto uma observação influenciaria o meu modelo se eu a removesse. Uma distância de cook alta significa que aquele valor é influente. Distâncias acima de 1 são influentes. Distâncias acima de 4/N é sinal de alerta. N é o número de observações no estudo.
Leaverage: valor ideal é 1+ K / N. K é o número de preditores e N o número de observações no estudo. Se 2 ou 3 vezes maior do que isso, atenção. O leaverage significa o quão extremo é uma observação em relação aos valores das variáveis preditores. Ele mede o tanto uma observação pode puxar a curva S da regressao logística.
-
DFBETA: ele mede o quanto o coeficiente de um determinado preditor muda quando uma determinada observação é retirada. Uma observação é geralmente considerada influente se exceder \(\pm 2 / \sqrt{N}\). No entanto, grande parte dos pesquisadores e algumas literaturas consideram que o DFBETA deve ser menor do que 1.
Para verificar os valores extremos (outliers):
- Resíduos padronizados: são considerados valores extremos aqueles que estão fora da faixa \(-2/2\). Se você quiser ser mais conservador, pode adotar a faixai entre \(-3/3\).
Calculando distância de Cook
A forma mais simples é plotar a distância de cook e observar valores que sejam maiores do que 1
Observe na Figura 6 que não existe nenhum valor acima de 1. Você também pode verificar se tem algum valor de alterta (aqueles maiores do que \(4/N\)). No caso, o valor seria 0.01`
Você pode avaliar se, no seu modelo, existem valores considerados influentes pela distância de cook. Para isso, você pode utilizar a função check_outliers()do pacote performance. Vamos fazer:
Calculando o Leaverage
Leverage (ou alavancagem) é uma medida que indica quão distante uma observação está das outras observações em termos das variáveis preditoras. Em outras palavras, o leverage mede o potencial de uma observação influenciar os resultados do modelo devido à sua posição no espaço das variáveis independentes. Portanto, observações com alto leverage estão longe do centro das variáveis preditoras. Isso significa que elas têm o potencial de exercer uma influência desproporcional sobre o ajuste do modelo.
Podemos também calcular esse parâmetro de uma forma semelhante à distância de Cook e apresentar em um plot.
A Figura 7 demonstra o valor de Leaverage de cada unidade de análise do estudo. Pode-se avaliar pontos influentes buscando valores acima de (1 + K / N), ou seja, valores acima de 1.
Calculando o DFBETA
O DFBETA pode ser calculado pela função abaixo.
Veja que a função acima gerou uma tabela em que eu tenho o DFBETA do interepto e dos dois preditores (glicose e idade). Agora, verificar se tem valores maiores do que \(\pm 2 / \sqrt{N}\) para cada preditor.
Código
E agora, para idade
Encontrando outliers
Podemos procurar outliers calculando os resíduos padronizados e vendo quais estão fora da faixa -2,2 ou -3,3.
Na Figura 8 podemos observar que temos dois valores acima de 2, talvez no valor de 3. Esses valores não preocupam porque vimos que não temos valores influentes no modelo.
Podemos calcular os outliers dos resíduos padronizados da mesma forma que calculamos a distância de Cook, mas modificando um pouco os argumentos da função. Vejam abaixo:
Código
10 outliers detected: cases 76, 124, 183, 343, 350, 454, 460, 503, 667,
685.
- Based on the following method and threshold: zscore (3).
- For variable: (Whole model).
-----------------------------------------------------------------------------
Outliers per variable (zscore):
$glucose
Row Distance_Zscore
76 76 3.78119
183 183 3.78119
343 343 3.78119
350 350 3.78119
503 503 3.78119
$age
Row Distance_Zscore
124 124 3.040681
454 454 3.295778
460 460 4.061069
667 667 3.125714
685 685 3.040681Descrição dos resultados
Para descrição dos resultados, podemos utilizar uma função denominada report(). Além disso, vou, antes, calcular os pseudoR chamado Nagelkerke para que você entendam o seu significado. Você pode pegar o resultado do report e adaptar ao seu texto. Vejam abaixo.
We fitted a logistic model (estimated using ML) to predict diabetes with
glucose and age (formula: diabetes ~ glucose + age). The model's explanatory
power is moderate (Tjur's R2 = 0.25). The model's intercept, corresponding to
glucose = 0 and age = 0, is at -5.91 (95% CI [-6.85, -5.03], p < .001). Within
this model:
- The effect of glucose is statistically significant and positive (beta = 0.04,
95% CI [0.03, 0.04], p < .001; Std. beta = 1.14, 95% CI [0.94, 1.35])
- The effect of age is statistically significant and positive (beta = 0.02, 95%
CI [0.01, 0.04], p < .001; Std. beta = 0.29, 95% CI [0.12, 0.46])
Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald z-distribution approximation.Ajustamos um modelo logístico (estimado por máxima verossimilhança) para prever diabetes usando glicose e idade (fórmula: diabetes ~ glicose + idade). O modelo apresenta capacidade moderada de dsicrimanção entre grupos (R² de Tjur = 0,25), indicando um diferença média de 22% entre os desfechos tem ou não diabetes. Além disso, o modelo explica uma proporção da variabilidade dos dados em relação ao modelo nulo (R² de Nagelkerke = 0,31), sugerindo que glicose e idade contribuem significativamente para reduzir a incerteza na previsão. Dentro deste modelo:
- O efeito da glicose é estatisticamente significativo e positivo (beta = 0,04, IC 95% [0,03; 0,04], p < 0,001; beta padronizado = 1,14, IC 95% [0,94; 1,35]).
- O efeito da idade é estatisticamente significativo e positivo (beta = 0,02, IC 95% [0,01; 0,04], p < 0,001; beta padronizado = 0,29, IC 95% [0,12; 0,46]).
Avaliando o efeito dos preditores
Para avaliar os efeitos dos preditores (glucose e age) precisamos tranformá-los. Lembra que os coeficientes estão no formato de log de Odds ou log da chance, o que dificulta a intepretação. Para isso, temos que fazer o exponencial do valor de cada um dos coeficientes (de glucose e age). Para isso, podemos utilizar uma função chamada tab_model() do pacote sjPlot em que ele calcula isso e já entrega uma tabela formatada.
| diabetes | |||
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.00 | 0.00 – 0.01 | <0.001 |
| glucose | 1.04 | 1.03 – 1.04 | <0.001 |
| age | 1.03 | 1.01 – 1.04 | 0.001 |
| Observations | 768 | ||
| R2 Tjur | 0.247 | ||
Intepretação dos Odds Ratios
O efeito da glicose é estatisticamente significativo e positivo, com um odds ratio de 1,04 (IC 95%: [1,03; 1,04], p < 0,001). Isso significa que, para cada aumento de 1 unidade na glicose, a odds de ter diabetes aumenta em 4% (mantendo a idade constante). Em outras palavras, indivíduos com níveis mais altos de glicose têm uma chance 4% maior de desenvolver diabetes em comparação àqueles com níveis 1 unidade menores, ajustado para idade.
O efeito da idade também é estatisticamente significativo e positivo, com um odds ratio de 1,03 (IC 95%: [1,01; 1,04], p = 0,001). Isso indica que, para cada aumento de 1 ano na idade, a odds de ter diabetes aumenta em 3% (mantendo a glicose constante). Ou seja, indivíduos mais velhos têm uma chance 3% maior de desenvolver diabetes em comparação àqueles 1 ano mais novos, ajustado para glicose.