Atividade Regressão Logistica

Autor

Enoque Filho

Data de Publicação

14 de julho de 2024

O objetivo desse trabalho foi aplicar o método de regressão logistica em um conjunto de dados contendo uma variável binária (Sobreivente | Não sobrevivente) e outras variáveis independentes (Sexo, Classe e Idade).

SETUP

Código

# Pacotes ----------------------------------
  library(tidyverse) # manipulação dos dados e ggplot2
  library(ggpubr) # gráficos 
  library(flextable) # tabelas
  library(gtsummary) # tabelas
  library(patchwork) # gráficos
  library(broom) # 

# Configurações ----------------------------
  theme_set(theme_test())
  gtsummary::theme_gtsummary_language("pt")

IMPORTAÇÃO DOS DADOS

Após a importação dos dados as variáveis qualitativas foram recategorizadas, adicionando seus respectivies níveis.

Código

dados <-  
  read.csv("dados.csv", sep = '\t') %>% 
  mutate(
    sobrevivente = factor(ifelse(sobrevivente == 0, "Não", "Sim"),
                          levels = c("Sim", "Não")),  # ordem dos níveis
    sexo = factor(ifelse(sexo == 0, "Feminino",  "Masculino"),
                  levels = c("Feminino","Masculino")), # ordem dos níveis
    classe = factor(
      case_when(
        classe == 1 ~ "Alta",
        classe == 2 ~ "Média",
        classe == 3 ~ "Baixa",),
      levels = c("Alta", "Média", "Baixa")) # ordem dos níveis
    )

A variavel dependente será Sobrevivente tendo como nível de refererência o Sim. Portanto o modelo será ajustado para estimar razão de chances de Sobrevivência Não / Sobrevivencia Sim.

ANÁLISE DESCRITIVA

Código

plot_sobrevivente <-  dados %>% 
  count(sobrevivente) %>% 
  mutate(lab = paste0(n,' (', round(prop.table(n)*100, 2), '%)')) %>%
  ggplot(aes(x = sobrevivente, y = n, label = lab)) +
  geom_col(fill = 'darkgreen', alpha = 0.6) +
  geom_text(vjust = -0.5, size = 3) +
  labs(y = 'Frequência')

plot_sexo <-  dados %>% 
  count(sexo) %>% 
  mutate(lab = paste0(n,' (', round(prop.table(n)*100, 2), '%)')) %>%
  ggplot(aes(x = sexo, y = n, label = lab)) +
  geom_col(fill = 'darkgreen', alpha = 0.6) +
  geom_text(vjust = -0.5, size = 3) +
  labs(y = 'Frequência')

plot_classe <-  dados %>% 
  count(classe) %>% 
  mutate(lab = paste0(n,' (', round(prop.table(n)*100, 2), '%)')) %>%
  ggplot(aes(x = classe, y = n, label = lab)) +
  geom_col(fill = 'darkgreen', alpha = 0.6) +
  geom_text(vjust = -0.5, size = 3) +
  labs(y = 'Frequência')

plot_sobrevivente + plot_sexo + plot_classe

Código

  p1 <- ggplot(dados, aes(y = '', x =  idade)) +
    geom_boxplot(fill = "gray70", alpha = 0.5) +
    theme(axis.ticks.x = element_blank(),
          axis.text.y  = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text    = element_blank(),
          axis.title.x = element_blank()
    ) +
    labs(y = '')

  p2 <- ggplot(dados, aes(x = idade, after_stat(density))) +
    geom_histogram(color = "white", fill = "gray70", bins = 30, alpha = 0.5) +
    geom_density(color = "black", linewidth = 0.5) +
    geom_vline(aes(xintercept = mean(   idade, na.rm = TRUE)),
               linetype = 'dashed', color = 'red', linewidth = 0.8) +
    labs(y = 'Densidade')


  (p1 / p2) + 
    plot_layout(heights = c(1,5)) &
    scale_x_continuous(
      limits = c(min(pull(dados, idade), na.rm = TRUE) - 3,
                 max(pull(dados, idade), na.rm = TRUE) + 3)
      )

Código

dados %>% 
  gtsummary::tbl_cross(
    sexo, sobrevivente, 
    statistic = "{n} ({p}%)",
    percent = "column"
    ) %>% 
  bold_labels() %>% 
  as_flex_table()

	sobrevivente
	Sim	Não	Total
sexo
Feminino	292 (68%)	96 (16%)	388 (37%)
Masculino	135 (32%)	523 (84%)	658 (63%)
Total	427 (100%)	619 (100%)	1,046 (100%)

Código

count(dados, sexo, sobrevivente) %>% 
  ggbarplot(y = "n", x = "sexo", fill = "sobrevivente", 
            position = position_fill(), label = TRUE, lab.pos = "in", 
            orientation = "horizontal", lab.hjust = 1.6)

Código

dados %>% 
  gtsummary::tbl_cross(
    classe, sobrevivente, 
    statistic = "{n} ({p}%)",
    percent = "column"
    ) %>% 
  bold_labels() %>% 
  as_flex_table()

	sobrevivente
	Sim	Não	Total
classe
Alta	181 (42%)	103 (17%)	284 (27%)
Média	115 (27%)	146 (24%)	261 (25%)
Baixa	131 (31%)	370 (60%)	501 (48%)
Total	427 (100%)	619 (100%)	1,046 (100%)

Código

count(dados, classe, sobrevivente) %>% 
  ggbarplot(y = "n", x = "classe", 
            fill = "sobrevivente", position = position_fill(), 
            label = TRUE, lab.pos = "in", lab.hjust = 1.6,
            orientation = "horizontal")

Código

ggviolin(
  dados,
  x = "sobrevivente", 
  y = "idade", orientation = "horizontal", 
  fill = "sobrevivente", 
  add = "boxplot"
  )

MODELAGEM

Ajustes de diferentes modelos

Modelo completo e modelo escolhido pelo método stepwise

Código

# Modelo com todas as variáveis
  modelo_completo <- glm(
    data = dados, 
    formula = sobrevivente ~ ., 
    family = binomial()
    )

# Modelo stepwise
  modelo_step <- modelo_completo %>% step()

Start:  AIC=992.45
sobrevivente ~ sexo + idade + classe

         Df Deviance     AIC
<none>        982.45  992.45
- idade   1  1013.80 1021.80
- classe  2  1101.34 1107.34
- sexo    1  1255.69 1263.69

Critério para selecionar o melhor modelo.

O método stepwise apontou para o modelo completo.

Código

list(modelo_completo, modelo_step) %>% 
  map(.f = function(x){
    broom::glance(x) %>% 
      select(AIC, BIC, Deviance = deviance)
    }) %>% 
  bind_rows() %>% 
  mutate(.before = 1, 
         Modelo = c("modelo_completo","modelo_step")) %>% 
  arrange(AIC) %>% 
  flextable() %>% 
  bold(i = 1, part = "header") %>% 
  autofit()

Modelo	AIC	BIC	Deviance
modelo_completo	992.4531	1,017.217	982.4531
modelo_step	992.4531	1,017.217	982.4531

DIAGNÓSTICO DO MODELO

Código

tabela_modelo = 
tbl_regression(modelo_completo, exponentiate = TRUE) %>%
    add_glance_table(
          include = c(AIC, BIC, deviance))  %>%
    modify_column_unhide(column = std.error) %>% 
  bold_labels() %>% 
  as_flex_table()

tabela_modelo

Características	OR	SE	95% IC	Valor-p
sexo
Feminino	—	—	—
Masculino	12.2	0.166	8.83, 16.9	<0.001
idade	1.03	0.006	1.02, 1.05	<0.001
classe
Alta	—	—	—
Média	3.60	0.226	2.32, 5.63	<0.001
Baixa	9.87	0.226	6.39, 15.5	<0.001
AIC	992
BIC	1,017
Deviance	982
Abreviações: IC = Intervalo de Confiança, OR = Razão de chances, SE = Erro padrão

Análise de Deviance

Código

anova(modelo_completo, test = "Chisq") %>% 
  tidy() %>%
  mutate(p.value = scales::pvalue(p.value)) %>% 
  flextable() %>% 
  bold(part = "header") %>% 
  colformat_double(j = 3:5, digits = 3)

term	df	deviance	df.residual	residual.deviance	p.value
NULL			1,045	1,414.620
sexo	1	312.612	1,044	1,102.008	<0.001
idade	1	0.669	1,043	1,101.339	0.413
classe	2	118.886	1,041	982.453	<0.001

Análise de Residuos

Pacote performance

Autocorrelação

Código

modelo_completo %>% performance::check_autocorrelation()

OK: Residuals appear to be independent and not autocorrelated (p = 0.062).

Outliers

Código

modelo_completo %>% performance::check_outliers()

OK: No outliers detected.
- Based on the following method and threshold: cook (0.5).
- For variable: (Whole model)

Multicolinearidade

Código

modelo_completo %>% performance::check_collinearity()

# Check for Multicollinearity

Low Correlation

   Term  VIF   VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
   sexo 1.05 [1.01, 1.19]     1.03      0.95     [0.84, 0.99]
  idade 1.35 [1.26, 1.47]     1.16      0.74     [0.68, 0.79]
 classe 1.41 [1.32, 1.54]     1.19      0.71     [0.65, 0.76]

Teste de resíduos

Código

# Resíduos de Pearson
tibble(
  Pearson = residuals(modelo_completo, type = "pearson"),
  Deviance = residuals(modelo_completo, type = "deviance"),
) %>% 
  pivot_longer(everything(), names_to = "Residuo") %>% 
  summarise(SQR = sum(value^2), .by = Residuo) %>% 
  mutate(Pchisq = pchisq(SQR, df = modelo_completo$df.residual, lower.tail = F)) %>% 
  flextable() %>% 
  autofit() %>% 
  bold(part = "header")

Residuo	SQR	Pchisq
Pearson	1,084.9937	0.1671558
Deviance	982.4531	0.9020470

Acurácia

Código

modelo_completo %>% performance::performance_accuracy()

# Accuracy of Model Predictions

Accuracy (95% CI): 83.81% [80.14%, 86.01%]
Method: Area under Curve

Porcentagem de predições corretas

Código

modelo_completo %>% performance::performance_pcp()

# Percentage of Correct Predictions from Logistic Regression Model

  Full model: 69.86% [67.08% - 72.64%]
  Null model: 51.68% [48.66% - 54.71%]

# Likelihood-Ratio-Test

  Chi-squared: 432.167
  df:   4.000
  p-value:   0.000

Curva ROC

Código

Epi::ROC(form=dados$sobrevivente~ dados$sexo + dados$idade +  dados$classe, plot="ROC",MI=F)

Poder preditivo do modelo.

Sensibilidade e Especificidade

Código

modelo_completo %>% 
  augment(type.predict = "response") %>%
  mutate(
    Predito = factor(ifelse(.fitted <= 0.63, "Sim", "Não"), 
                     levels = c("Sim","Não"))
    ) %>% 
  tbl_cross(row = sobrevivente, col = Predito, percent =  "row") %>% 
  bold_labels() %>% 
  as_flex_table()

	Predito
	Sim	Não	Total
sobrevivente
Sim	340 (80%)	87 (20%)	427 (100%)
Não	146 (24%)	473 (76%)	619 (100%)
Total	486 (46%)	560 (54%)	1,046 (100%)

Sensibilidade: 76% (Verdadeiro positivo)
Especificidade: 80% (Verdadeiro negativo)

Poder de previsão

Código

modelo_completo %>% 
  augment(type.predict = "response") %>% 
  mutate(Predito = factor(ifelse(.fitted <= 0.63, "Sim", "Não"), levels = c("Sim", "Não"))) %>%
  tbl_cross(row = sobrevivente, col = Predito, percent =  "column") %>% 
  bold_labels() %>%
  as_flex_table()

	Predito
	Sim	Não	Total
sobrevivente
Sim	340 (70%)	87 (16%)	427 (41%)
Não	146 (30%)	473 (84%)	619 (59%)
Total	486 (100%)	560 (100%)	1,046 (100%)

Poder de previsão de não sobrevivente: 84%
Poder de previsão sobrevivente: 70%

INTERPRETAÇÃO DOS RESULTADOS

Código

tabela_modelo

Características	OR	SE	95% IC	Valor-p
sexo
Feminino	—	—	—
Masculino	12.2	0.166	8.83, 16.9	<0.001
idade	1.03	0.006	1.02, 1.05	<0.001
classe
Alta	—	—	—
Média	3.60	0.226	2.32, 5.63	<0.001
Baixa	9.87	0.226	6.39, 15.5	<0.001
AIC	992
BIC	1,017
Deviance	982
Abreviações: IC = Intervalo de Confiança, OR = Razão de chances, SE = Erro padrão

Individuos do sexo Masculino possuem uma chance 12 vezes maior de pertencerem ao grupo de não sobrevivente quando comparado com os do sexo Feminino.
Idade está pouco relacionada a pertença de um ou outro grupo, com uma leve tendencia a estar relacionada com a pertença ao grupo de não sobrevivente.
Individuos de classe Média ou Baixa estão mais relacionadas a pertença do grupo de não sobreiventes quando comparadas com a classe Alta. a chance é 3,6 vezes maior nos de classe Média e quase 10 vezes maior nos de classe baixa quando comparados com a Classe Alta.

REFERÊNCIAS

Cientística & Podcast Naruhodo: Estatística Psicobio I 2022 #30 https://www.youtube.com/live/yr657G1N7GQ?si=GSnCydg1S17aDx-V
Cientística & Podcast Naruhodo: Estatística Psicobio I 2022 #31 https://www.youtube.com/live/lCLEaVOvo-c?si=JG0NSbI3YoEh4W3R
Cientística & Podcast Naruhodo: Estatística Psicobio I 2022 #32 https://www.youtube.com/live/lCLEaVOvo-c?si=iMqFh39XcGo5R0kt
R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
Sjoberg DD, Whiting K, Curry M, Lavery JA, Larmarange J. Reproducible summary tables with the gtsummary package. The R Journal 2021;13:570–80. https://doi.org/10.32614/RJ-2021-053.
Kassambara A (2023). ggpubr: ‘ggplot2’ Based Publication Ready Plots. R package version 0.6.0, https://CRAN.R-project.org/package=ggpubr.
Lüdecke et al., (2021). performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
Pedersen T (2024). patchwork: The Composer of Plots. R package version 1.2.0, https://CRAN.R-project.org/package=patchwork.
Gohel D, Skintzos P (2023). flextable: Functions for Tabular Reporting. R package version 0.9.4, https://CRAN.R-project.org/package=flextable.