8  Módulo 6 | Modelos estatísticos

Índice

• 9 SETUP
• 10 TABELAS
  • 10.1 Criar a tabela
  • 10.2 Adaptar a tabela
• 11 TESTES ESTATÍSTICOS
  • 11.1 Comparar médias
    • 11.1.1 Test t (comparar médias)
  • 11.2 Comparar proporções
    • 11.2.1 Teste Chi-Quadrado
• 12 MODELOS DE REGRESSÃO
  • 12.1 MODELOS LINEARES
    • 12.1.1 Criar intuição
    • 12.1.2 Visualizar a previsão do modelo
    • 12.1.3 Dataset Real
    • 12.1.4 Regressão Univariada
    • 12.1.5 Regressão Multivariada
    • 12.1.6 Diagnostico de modelos
    • 12.1.7 Visualização dos efeitos do modelo
  • 12.2 MODELOS NÃO LINEARES
    • 12.2.1 Logit model
  • 12.3 SURVIVAL ANALYSIS
    • 12.3.1 Kaplan-Meier plots
    • 12.3.2 Survival times
    • 12.3.3 Cox regression
• 13 FIM

9 SETUP

Limpar o ambiente

rm(list = ls(all.names = TRUE)) # limpa todos os objetos, incluindo os ocultos

Definir a semente (seed)

set.seed(123)

Carregar as bibliotecas

required_packages <- c("janitor",
                       "tidyverse",
                       "rio",
                       "data.table",
                       "broom",
                       "ggplot2",
                       "gt",
                       "gtsummary",
                       "MLDataR",
                       "rstatix",
                       "cards",
                       "survival",
                       "ggstatsplot",
                       "ggsurvfit",
                       "insight",
                       "cardx",
                       "broom.helpers",
                       "skimr"
                       )      

for (pkg in required_packages) {
  # install packages if not already present
  if (!pkg %in% rownames(installed.packages())) {
    install.packages(pkg, repos = "https://cloud.r-project.org", dependencies=T)
  }
  # load packages to this current session 
  library(pkg, character.only = TRUE)
}
remove(required_packages)
remove(pkg)

10 TABELAS

Importante

O R(Posit) é uma ferramenta poderosa para produzir tabelas de forma muito customizavel. O package com maior potencialidade é o GT. ver mais aqui https://gt.rstudio.com/articles/gt.html

Outra sugestão de recurso de aprendizagem https://gt.albert-rapp.de/

Dataset de 918 doentes e com 10 variaveis:

Age: age of the patient presenting with heart disease Sex: gender of the patient RestingBP: blood pressure for resting heart beat Cholesterol: Cholesterol reading FastingBS: blood sample of glucose after a patient fasts RestingECG: Resting echocardiography is an indicator of previous myocardial infarction e.g. heart attack MaxHR: Maximum heart rate Angina: chest pain caused by decreased flood flow https://www.nhs.uk/conditions/angina/ HeartPeakReading: reading at the peak of the heart rate HeartDisease: the classification label of whether patient has heart disease or not

Frequentemente os relatórios pedem tabelas descritivas pela variável de exposição o “gtsummary” ajuda a criar essas tabelas e pode ser complementado pelo package GT

# Carregar o conjunto de dados de doença cardíaca do pacote MLDataR
hd <- MLDataR::heartdisease

# Limpar nomes das colunas para consistência e legibilidade
hd <- janitor::clean_names(hd)
# A função clean_names do pacote janitor padroniza o formato dos nomes das colunas,
# tornando-os mais consistentes e fáceis de trabalhar.

hd <- hd %>%
  mutate(heart_disease_yn = case_when(
    heart_disease == 0 ~ "No",
    heart_disease == 1 ~ "Yes",
    TRUE ~ as.character(heart_disease) # for safety, in case there are other values
  ),
  Age=age,
  Sex=sex,
  "Blood Pressure"=resting_bp,
  Cholesterol=cholesterol,
  Angina=angina)

10.1 Criar a tabela

# Criar uma tabela resumo de variáveis selecionadas dividida pelo estado da doença cardíaca
# O pacote gtsummary é utilizado para criar tabelas prontas para publicação
table1 <- tbl_summary(
  hd,
  include = c(Age, Sex, "Blood Pressure", Cholesterol, Angina), # Selecionar variáveis específicas para o resumo
  by = heart_disease_yn, # Dividir a tabela pelo estado da doença cardíaca
  missing = "no", # Excluir dados ausentes do resumo
) |> 
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Heart Failure**") %>%
  add_n() |>  # Adicionar uma coluna para contar observações não ausentes
  add_p() |>  # Realizar testes estatísticos para comparar grupos
  modify_header(label = "**Variable**") |>  # Personalizar o cabeçalho da tabela
  bold_labels() # Tornar as etiquetas em negrito para ênfase

# Visualizar a tabela criada
table1

Variable	N	Heart Failure		p-value2
		No N = 4101	Yes N = 5081
Age	918	51 (43, 57)	57 (51, 62)	<0.001
Sex	918			<0.001
F		143 (35%)	50 (9.8%)
M		267 (65%)	458 (90%)
Blood Pressure	918	130 (120, 140)	132 (120, 145)	<0.001
Cholesterol	918	227 (197, 267)	217 (0, 267)	<0.001
Angina	918			<0.001
N		355 (87%)	192 (38%)
Y		55 (13%)	316 (62%)
1 Median (Q1, Q3); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test

10.2 Adaptar a tabela

Melhorar a tabela com o GT

table1_f <- table1 |> 
  as_gt() |>  
  tab_header(
    title = md("**Risk model for Heart Failure**") ,
    subtitle = "By Age and Sex"
  ) |>
  tab_source_note(
    source_note = md("Source: MlR dataset")
  ) |>
  fmt_number( decimals = 3) |>   
  opt_stylize(style = 1, color = "gray") |> 
  opt_align_table_header(align = "left") 
# Show the gt Table

table1_f

Risk model for Heart Failure
By Age and Sex
Variable	N	Heart Failure		p-value2
		No N = 4101	Yes N = 5081
Age	918	51 (43, 57)	57 (51, 62)	0.000
Sex	918			0.000
F		143 (35%)	50 (9.8%)
M		267 (65%)	458 (90%)
Blood Pressure	918	130 (120, 140)	132 (120, 145)	0.001
Cholesterol	918	227 (197, 267)	217 (0, 267)	0.000
Angina	918			0.000
N		355 (87%)	192 (38%)
Y		55 (13%)	316 (62%)
Source: MlR dataset
1 Median (Q1, Q3); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test

Gravar a tabela

#table1_f |> gtsave("tab_1.png", expand = 10)
table1_f |> gtsave("tab_1.docx")

11 TESTES ESTATÍSTICOS

11.1 Comparar médias

library(rstatix)

11.1.1 Test t (comparar médias)

Realiza um teste t de duas amostras não pareadas (ou teste t independente) no R para comparar as médias de idade (age) entre dois grupos de sexo (sex) no conjunto de dados hd. Este teste é útil para determinar se há uma diferença estatisticamente significativa entre as médias dos dois grupos.

# Two-samples unpaired test
#:::::::::::::::::::::::::::::::::::::::::
hd %>% t_test(age ~ sex)

# A tibble: 1 × 8
  .y.   group1 group2    n1    n2 statistic    df      p
* <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>  <dbl>
1 age   F      M        193   725     -1.68  299. 0.0945

11.2 Comparar proporções

11.2.1 Teste Chi-Quadrado

Realiza um teste qui-quadrado (χ²) para avaliar a associação entre gênero e ocorrência de doença cardíaca no conjunto de dados. Vou comentar o código em português de Portugal para explicar cada parte.

# Example contingency table
hd_gender <- matrix(c(458, 50, 267, 143), nrow = 2,
                         dimnames = list(Gender = c("Male", "Female"),
                                         Smoking = c("hd", "no_hd")))

# Perform chi-square test
chisq.test(hd_gender)

    Pearson's Chi-squared test with Yates' continuity correction

data:  hd_gender
X-squared = 84.145, df = 1, p-value < 2.2e-16

12 MODELOS DE REGRESSÃO

12.1 MODELOS LINEARES

A regressão linear é um método estatístico utilizado para explorar a relação entre duas variáveis contínuas. Em termos simples, é como encontrar a melhor linha reta que pode ser traçada através de um conjunto de pontos de dados em um gráfico. Esta linha representa a relação média entre as duas variáveis.

Por exemplo, vamos considerar a relação entre idade e pressão arterial. Neste caso, a idade é a variável independente (aquela que pensamos que pode estar causar mudanças), e a pressão arterial é a variável dependente (aquela em que estamos a observar mudanças).

A regressão linear ajuda a entender, em média, quanto a pressão arterial muda para cada ano de idade (coeficiente da regressão para a idade). Isso não significa que a pressão arterial de cada indivíduo seguirá exatamente este padrão, mas em média, podemos esperar uma certa tendência.

O erro entre a previsão do modelo e o valor observado é o erro ou residuo do modelo.

12.1.1 Criar intuição

1. Criar um dataset de 100 observações onde somos nos que definimos os componentes.

# Gerar 100 observações random de 25 a 75 anos com replacement 
age <- sample(25:75, 100, replace = TRUE)

# Criar ruido (normal distribution)
noise <- rnorm(100, mean = 0, sd = 5) # Adjust sd (standard deviation) as needed


# calcular a pressão alterial aumentando 0.5 por cada ano de idade 
blood_pressure <- 100 + 0.5 * age + noise

# Combine into a data frame
data <- data.frame(Age = age, BloodPressure = blood_pressure)

# View the first few rows of the data frame
head(data)

  Age BloodPressure
1  66      125.6392
2  30      117.6448
3  66      145.2720
4  66      124.7533
5  62      134.9790
6  72      134.2935

# Criar um gráfico de dispersão com as variáveis Age (Idade) e BloodPressure (Pressão Sanguínea)
# A função ggplot() inicia a construção do gráfico
# aes() define as variáveis a serem utilizadas nos eixos x e y
# geom_point() adiciona pontos ao gráfico, representando as observações individuais

data |> 
  ggplot(aes(x = Age, y = BloodPressure)) + 
  geom_point(alpha = 0.5) # O parâmetro alpha controla a transparência dos pontos, melhorando a visibilidade em áreas de sobreposição

2. Criar um modelo de regressão linear

# Realizar uma regressão linear utilizando a função lm()
# lm(formula, data): onde 'formula' define o modelo e 'data' é o conjunto de dados
# formula = var_dependente ~ var_independente: sintaxe para especificar o modelo

# Criar o modelo de regressão linear onde 'BloodPressure' é a variável dependente 
# e 'Age' é a variável independente, usando o conjunto de dados 'data'
model1 <- lm(BloodPressure ~ Age, data = data)

# Visualizar um resumo do modelo de regressão linear
# summary(model) fornece detalhes como coeficientes, erro padrão, valor-p, R-quadrado, etc.
summary(model1)

Call:
lm(formula = BloodPressure ~ Age, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.5217  -3.8387   0.2271   4.1783  12.5513 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 99.46561    2.05356   48.44   <2e-16 ***
Age          0.50562    0.03745   13.50   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.108 on 98 degrees of freedom
Multiple R-squared:  0.6503,    Adjusted R-squared:  0.6468 
F-statistic: 182.3 on 1 and 98 DF,  p-value: < 2.2e-16

# Carregar o pacote broom para trabalhar com modelos estatísticos
# library(broom) - Ja estava carregado no inicio do script

# Utilizar a função tidy do pacote broom para converter a saída do modelo de regressão linear (model1) em um data frame 'arrumado'
# conf.int = TRUE inclui o intervalo de confiança para os coeficientes do modelo
model1_tidy <- tidy(model1, conf.int = TRUE)

# Visualizar o data frame resultante
model1_tidy

# A tibble: 2 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   99.5      2.05        48.4 2.93e-70   95.4     104.   
2 Age            0.506    0.0375      13.5 4.31e-24    0.431     0.580

Podemos ver que a estimativa para o coefeciente da idade é de 0.51, um valor muito proximo do que establecemos de 0.5. Pode ser interpretado como que por cada ano de vida a mais a pressão arterial aumenta em méda 0,51 mmHg. Observamos também que o verdadeiro valor pode variar entre 0.44 e 0.58, e nos sabemos este intervalo contem o verdadeiro valor de 0.5.

Vamos agora colocar a previsão do modelo no grafico

# Carregar o pacote broom para trabalhar com modelos estatísticos
# library(broom)

# Utilizar a função augment do pacote broom para adicionar as previsões e resíduos do modelo ao dataset original
# model1_augment será um novo dataset que inclui tanto os dados originais como as previsões do modelo e outras estatísticas relevantes
model1_augment <- augment(model1)

# Visualizar o novo dataset
head(model1_augment)

# A tibble: 6 × 8
  BloodPressure   Age .fitted .resid   .hat .sigma .cooksd .std.resid
          <dbl> <int>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>
1          126.    66    133.  -7.20 0.0189   5.08 0.0195      -1.42 
2          118.    30    115.   3.01 0.0387   5.12 0.00728      0.601
3          145.    66    133.  12.4  0.0189   4.97 0.0583       2.46 
4          125.    66    133.  -8.08 0.0189   5.07 0.0246      -1.60 
5          135.    62    131.   4.17 0.0142   5.12 0.00487      0.821
6          134.    72    136.  -1.58 0.0292   5.13 0.00147     -0.313

12.1.2 Visualizar a previsão do modelo

# Criar um gráfico para diagnosticar o modelo de regressão linear
model1_plot <- ggplot() +
  # Adicionar pontos representando os dados originais
  geom_point(
    data = model1_augment,
    mapping = aes(x = Age, y = BloodPressure),
    alpha = 0.5  # Ajustar a transparência para melhor visualização
  ) +
  # Adicionar uma linha representando as previsões do modelo
  geom_line(
    data = model1_augment,
    mapping = aes(x = Age, y = .fitted)
  ) +
  # Adicionar títulos e etiquetas aos eixos
  labs(
    title = "Diagnóstico do Modelo 1",
    x = "Idade",
    y = "Pressão Sanguínea"
  )

# Visualizar o gráfico
model1_plot

12.1.3 Dataset Real

Dataset de 918 doentes e com 10 variaveis:

Age: age of the patient presenting with heart disease Sex: gender of the patient RestingBP: blood pressure for resting heart beat Cholesterol: Cholesterol reading FastingBS: blood sample of glucose after a patient fasts RestingECG: Resting echocardiography is an indicator of previous myocardial infarction e.g. heart attack MaxHR: Maximum heart rate Angina: chest pain caused by decreased flood flow https://www.nhs.uk/conditions/angina/ HeartPeakReading: reading at the peak of the heart rate HeartDisease: the classification label of whether patient has heart disease or not

# Carregar o conjunto de dados 'heartdisease' do pacote MLDataR
hd <- MLDataR::heartdisease

# Utilizar a função clean_names do pacote janitor para padronizar os nomes das colunas do conjunto de dados
# Esta função transforma todos os nomes das colunas para letras minúsculas e substitui espaços e caracteres especiais por underscores
hd <- janitor::clean_names(hd)

Será que a idade está associada à pressão arterial? Podemos explorar esta possível associação

1. Podemos fazer uma primeira exploração visual

hd |> ggplot(aes(age, resting_bp)) + 
  geom_point(alpha = 0.5)

12.1.4 Regressão Univariada

\[ BP = \text{Intercept} + (\text{Slope}_1 \times \text{Age}) + \text{Error}\]

Criar um modelo de regressão linear

# Utilizar a função lm() para realizar uma regressão linear
# lm(formula, data): 'formula' define o modelo e 'data' é o conjunto de dados
# formula = var_dependente ~ var_independente: formato para especificar o modelo

# Criar o modelo de regressão linear onde 'resting_bp' é a variável dependente
# e 'age' é a variável independente, usando o conjunto de dados 'hd'
model2 <- lm(resting_bp ~ age, data = hd)

# Obter um resumo estatístico do modelo de regressão linear
# summary(model) fornece detalhes como coeficientes, erro padrão, valor-p, R-quadrado, etc.
summary(model2)

Call:
lm(formula = resting_bp ~ age, data = hd)

Residuals:
     Min       1Q   Median       3Q      Max 
-133.140  -11.143   -1.151    9.357   67.359 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 105.67692    3.40781  31.010  < 2e-16 ***
age           0.49933    0.06272   7.961 5.01e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.91 on 916 degrees of freedom
Multiple R-squared:  0.06472,   Adjusted R-squared:  0.0637 
F-statistic: 63.38 on 1 and 916 DF,  p-value: 5.015e-15

# Podemos utilizar a função tidy para criar um dataframe com o modelo
model2_tidy <- tidy(model2, conf.int = TRUE)
model2_tidy

# A tibble: 2 × 7
  term        estimate std.error statistic   p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)  106.       3.41       31.0  6.31e-145   99.0     112.   
2 age            0.499    0.0627      7.96 5.01e- 15    0.376     0.622

Vamos agora colocar a previsão do modelo no grafico

# Criar um dataset que incliu as previsões do modelo  
model2_augment <- augment(model2)

model2_plot <- ggplot() + 
    geom_point( data = model2_augment, 
  mapping = aes(
    x=age, 
    y=resting_bp), alpha=0.3
  ) +
  geom_line(
    data = model2_augment, 
  mapping = aes(
    x=age, 
    y=.fitted)
  ) +
    labs(
      title = "Model 2 diagnostics",
      x = "Age", 
      y = "Resting Blood Pressure") 

model2_plot

12.1.5 Regressão Multivariada

Permite a previsão de uma variável dependente (Pressão Arterial) com base em múltiplas variáveis independentes. Em termos simples, é como ajustar um plano (ou hiperplano, em dimensões superiores) mais adequado aos pontos de dados num espaço multidimensional, onde cada dimensão representa uma das variáveis independentes.

Adicionar o Sexo ao Modelo:

Vamos adicionar o sexo como outra variável independente, juntamente com a idade, para prever a pressão arterial. O sexo é uma variável categórica (tipicamente masculino ou feminino), por isso é tratado de forma ligeiramente diferente na análise de regressão em comparação com variáveis contínuas como a idade.

\[ BP = \text{Intercept} + (\text{Slope}_1 \times \text{Age}) + (\text{Slope}_2 \times \text{Sex}) + \text{Error}\]

Podemos fazer uma primeira exploração visual

# Utilizar ggplot para criar um gráfico de dispersão com os dados de 'hd'
hd |> 
  ggplot(aes(x = age, y = resting_bp, color = sex)) +  # Definir 'age' e 'resting_bp' como variáveis nos eixos x e y, respectivamente, e 'sex' para a cor dos pontos
  geom_point(alpha = 0.7)  # Adicionar pontos ao gráfico com transparência para melhor visualização de sobreposições

Regressão

# Utilizar a função lm() para criar um modelo de regressão linear
# lm(formula, data): 'formula' define o modelo e 'data' é o conjunto de dados
# formula = var_dependente ~ var_independente + var_independente: formato para especificar o modelo com múltiplas variáveis independentes

# Criar o modelo de regressão linear com 'resting_bp' como a variável dependente
# e 'age' (idade) e 'sex' (sexo) como variáveis independentes
# 'as.factor(sex)' converte 'sex' em uma variável categórica
model3 <- lm(resting_bp ~ age + as.factor(sex), data = hd)

# Podemos utilizar a função tidy para criar um dataframe com o modelo
model3_tidy <- tidy(model3, conf.int = TRUE)
model3_tidy

# A tibble: 3 × 7
  term            estimate std.error statistic   p.value conf.low conf.high
  <chr>              <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)      106.       3.54      29.9   1.21e-137   99.0     113.   
2 age                0.500    0.0628     7.96  5.04e- 15    0.377     0.624
3 as.factor(sex)M   -0.412    1.45      -0.284 7.77e-  1   -3.27      2.44 

Mantendo a idade constante o sexo masculino tem em média um valor de pressão arterior 0.4mmHg inferir ao das mulheres. No entanto o verdadeiro valor pode variar entre -3.2 e +2.4mmHg.

12.1.6 Diagnostico de modelos

# Carregar o pacote 'performance' para avaliação de modelos estatísticos
library(performance)

# Verificar as suposições do modelo de regressão linear 'model3'
# Esta função verifica várias suposições importantes, como normalidade dos resíduos, homocedasticidade, entre outras.
check_model(model3)

check_model() é uma função útil para diagnosticar se um modelo de regressão linear cumpre as suposições fundamentais. Fornece gráficos para avaliar aspectos como a normalidade dos resíduos e a homocedasticidade (variância constante dos resíduos).

model_performance(model3)

# Indices of model performance

AIC    |   AICc |    BIC |    R2 | R2 (adj.) |   RMSE |  Sigma
--------------------------------------------------------------
7909.1 | 7909.1 | 7928.4 | 0.065 |     0.063 | 17.894 | 17.924

model_performance() retorna várias métricas de desempenho do modelo, como o R-quadrado, que indica a proporção da variância na variável dependente que é previsível a partir das variáveis independentes, e o RMSE (Root Mean Square Error), que mede a qualidade do ajuste do modelo.

compare_performance(model2, model3, verbose = FALSE)

# Comparison of Model Performance Indices

Name   | Model |  AIC (weights) | AICc (weights) |  BIC (weights) |    R2
-------------------------------------------------------------------------
model2 |    lm | 7907.2 (0.723) | 7907.2 (0.725) | 7921.6 (0.967) | 0.065
model3 |    lm | 7909.1 (0.277) | 7909.1 (0.275) | 7928.4 (0.033) | 0.065

Name   | R2 (adj.) |   RMSE |  Sigma
------------------------------------
model2 |     0.064 | 17.895 | 17.915
model3 |     0.063 | 17.894 | 17.924

compare_performance() compara métricas de desempenho entre dois modelos. Isso é útil para avaliar qual modelo apresenta melhor ajuste ou previsão. O argumento verbose = FALSE simplifica a saída, mostrando apenas as informações mais críticas.

12.1.7 Visualização dos efeitos do modelo

# Carregar o pacote ggstatsplot para visualização avançada de estatísticas
library(ggstatsplot)

# Criar um gráfico dos coeficientes do modelo de regressão 'model3'
p <- ggcoefstats(
  model3,
  conf.int = TRUE,                      # Incluir intervalo de confiança
  conf.level = 0.95,                    # Nível de confiança de 95%
  exclude.intercept = TRUE,             # Excluir o intercepto do gráfico
  xlab = "Effect size in mmHg",         # Etiqueta do eixo x
  ylab = "Variables",                   # Etiqueta do eixo y
  title = "Effect of Age and Sex on Blood Pressure",  # Título do gráfico
  subtitle = "Patients with Heart Failure",           # Subtítulo do gráfico
  caption = "MlR dataset",              # Legenda
  point.args = list(size = 3, color = "grey", na.rm = TRUE),  # Argumentos para os pontos
  errorbar.args = list(height = 0, na.rm = TRUE),            # Argumentos para as barras de erro
  vline = TRUE,                         # Incluir linha vertical em zero
  vline.args = list(linewidth = 1, linetype = "dashed"),     # Argumentos para a linha vertical
  package = "RColorBrewer",             # Pacote para a paleta de cores
  palette = "Dark2",                    # Escolha da paleta de cores
  ggtheme = ggstatsplot::theme_ggstatsplot()  # Tema do gráfico
) + xlim(-5, 5) + scale_y_discrete(labels = c("Sex", "Age"))  # Limites do eixo x e etiquetas personalizadas do eixo y

# Visualizar o gráfico
p

12.2 MODELOS NÃO LINEARES

Neste website temos acesso a uma tabela com os tipos de modelos que podemos usar de acordo com a nossa variavel dependente. https://strengejacke.github.io/regressionmodels/

12.2.1 Logit model

Um modelo de regressão logística, frequentemente referido como modelo logit, é usado quando a variável dependente é binária (ou seja, tem apenas dois possíveis resultados). Neste caso estamos interessados em prever as odds (probabilidade) de insuficiencia cardiaca(um resultado binário: sim 1 ou não 0) com base em duas variáveis independentes: idade e sexo.

Interpretação: - Os coeficientes \(\beta\) num modelo de regressão logística representam a alteração das log odds (logit) para cada alteração de uma unidade na variável preditora (idade ou sexo). - O exponencial do logit é um odds ratio (OR).

# Carregar o pacote glm para modelos generalizados
# library(glm)

# Criar um modelo de regressão logística para prever a ocorrência de doença cardíaca ('heart_disease')
# usando 'age' e 'sex' como variáveis preditoras
model4_logit <- glm(
  heart_disease ~ age + as.factor(sex),  # Definir a fórmula do modelo
  family = "binomial",                   # Especificar a função de ligação binomial para resposta binária
  data = hd                              # Usar o conjunto de dados 'hd'
)

# Carregar o pacote broom para manipulação de modelos estatísticos
library(broom)

# Utilizar a função tidy para criar um resumo do modelo
# 'exponentiate = TRUE' transforma os coeficientes em razões de probabilidades
# 'conf.int = TRUE' inclui os intervalos de confiança
model4_tidy <- tidy(model4_logit, exponentiate = TRUE, conf.int = TRUE)

Resultados em tabela

# Manipular o resumo do modelo 'model4_tidy'
model4_tidy <- model4_tidy |> 
  # Selecionar colunas específicas para o resumo
  select(term, estimate, p.value, conf.low, conf.high) |> 
  # Modificar os nomes dos termos para torná-los mais legíveis
  mutate(
    term = case_when(
      term == "(Intercept)" ~ "Intercept",   # Renomear '(Intercept)' para 'Intercept'
      term == "age" ~ "Age",                 # Renomear 'age' para 'Age'
      term == "as.factor(sex)M" ~ "Male"     # Renomear 'as.factor(sex)M' para 'Male'
    )
  ) |>  
  # Renomear colunas para maior clareza
  rename(
    variable = term,        # Renomear 'term' para 'variable'
    Odds_Ratio = estimate   # Renomear 'estimate' para 'Odds_Ratio'
  )

# Carregar o pacote gt para criar tabelas formatadas
# library(gt)

# Criar uma tabela com o pacote gt a partir do resumo do modelo 'model4_tidy'
tbl_m4 <- 
  gt(model4_tidy) |>
  # Adicionar um cabeçalho com título e subtítulo
  tab_header(
    title = md("**Modelo de Risco para Insuficiência Cardíaca**"),
    subtitle = "Por Idade e Sexo"
  ) |>
  # Adicionar uma nota de rodapé com a fonte dos dados
  tab_source_note(
    source_note = md("Fonte: Conjunto de dados MlR")
  ) |>
  # Renomear os rótulos das colunas
  cols_label(
    variable = "Preditor",
    Odds_Ratio = "Razão de Probabilidades",
    p.value = "Valor p",
    conf.low = "IC Inferior",
    conf.high = "IC Superior",
  ) |>
  # Formatar todos os números com três casas decimais
  fmt_number(decimals = 3) |>
  # Aplicar um estilo visual às células dos rótulos das colunas
  tab_style(
    style = cell_fill(color = "grey95"),
    locations = list(cells_column_labels())
  ) |>
  # Alinhar o cabeçalho da tabela à esquerda
  opt_align_table_header(align = "left") 

# Mostrar a tabela gt
tbl_m4

Modelo de Risco para Insuficiência Cardíaca
Por Idade e Sexo
Preditor	Razão de Probabilidades	Valor p	IC Inferior	IC Superior
Intercept	0.010	0.000	0.004	0.024
Age	1.069	0.000	1.052	1.086
Male	5.161	0.000	3.585	7.539
Fonte: Conjunto de dados MlR

Exportar tabela

tbl_m4 |> gtsave("output/documents/tbl_m4.docx")

12.3 SURVIVAL ANALYSIS

A análise de sobrevivência lida com dados de tempo até a ocorrência de um evento. Este tipo de análise é usado para estudar o tempo até a ocorrência de um evento de interesse, frequentemente referido como “falha” ou “morte”, embora o evento possa ser qualquer ponto final, como recorrência de doença, alta hospitalar ou falha de uma máquina. As principais características da análise de sobrevivência incluem:

Evento de Interesse: Este é um evento específico que o estudo foi desenhado para observar, como morte, recaída de uma doença, etc. Em contextos não médicos, isso pode ser a falha de um sistema mecânico, mudança de emprego, etc.

Censura: Um aspecto único da análise de sobrevivência é lidar com dados censurados. A censura ocorre quando temos informações incompletas sobre o tempo de sobrevivência de alguns indivíduos. Por exemplo, se um estudo termina e um participante ainda não experimentou o evento de interesse, seus dados são considerados censurados à direita. A censura também pode ocorrer se um participante abandonar o estudo.

Nota: a função Surv() no pacote {survival} aceita por padrão VERDADEIRO/FALSO, onde VERDADEIRO é evento e FALSO é censurado; 1/0 onde 1 é evento e 0 é censurado; ou 2/1 onde 2 é evento e 1 é censurado.

Os dados de evento devem ser formatados adequadamente.

Importante

Os dados muitas vezes vêm com datas de início e término em vez de tempos de sobrevivência pré-calculados. O primeiro passo é garantir que esses estejam formatados como datas no R. E depois calcular o tempo até ao evento

# Carregar o pacote survival para análise de sobrevivência
library(survival)

# Carregar o conjunto de dados 'lung' do pacote survival
data(lung, package = "survival")

# Modificar o conjunto de dados 'lung'
lung <- lung %>%
  # Recodificar a variável 'status': 1 alterado para 0, 2 alterado para 1
  mutate(
    status = recode(status, `1` = 0, `2` = 1)
  )

# Visualizar as primeiras linhas do conjunto de dados modificado
head(lung)

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      1  74   1       1       90       100     1175      NA
2    3  455      1  68   1       0       90        90     1225      15
3    3 1010      0  56   1       0       90        90       NA      15
4    5  210      1  57   1       1       90        60     1150      11
5    1  883      1  60   1       0      100        90       NA       0
6   12 1022      0  74   1       1       50        80      513       0

12.3.1 Kaplan-Meier plots

O método Kaplan-Meier é a maneira mais comum de estimar tempos de sobrevivência e probabilidades. Resulta em uma função degrau, onde há uma descida a cada vez que um evento ocorre.

Modelo de Sobrevivência: A função survfit2() com a fórmula Surv(time, status) ~ 1 cria um modelo de sobrevivência usando o tempo (time) e o status do evento (status) do conjunto de dados lung. O termo ~ 1 indica que é um modelo de sobrevivência geral, sem estratificação por variáveis adicionais.

Gráfico de Kaplan-Meier: ggsurvfit() é usado para criar um gráfico de Kaplan-Meier, mostrando a probabilidade de sobrevivência ao longo do tempo.

# Carregar o pacote ggsurvfit para visualização gráfica de dados de sobrevivência
library(ggsurvfit)

# Criar e visualizar uma curva de Kaplan-Meier para o conjunto de dados 'lung'
survfit2(Surv(time, status) ~ 1, data = lung) %>% 
  ggsurvfit() +
  # Adicionar rótulos e título ao gráfico
  labs(
    x = "Dias",  # Etiqueta do eixo x
    y = "Probabilidade de sobrevivência global",  # Etiqueta do eixo y
    title = "Probabilidade de sobrevivência para Câncer de Pulmão"  # Título do gráfico
  ) + 
  # Adicionar intervalo de confiança ao gráfico
  add_confidence_interval() +
  # Adicionar uma tabela de risco ao gráfico
  add_risktable()

Podemos calcular as curvas de sobrevivencia por sexo.

survfit2(Surv(time, status) ~ sex, data = lung) %>% 
  ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability", 
    title = "Survival probability for Lung Cancer"
  ) + 
  add_confidence_interval() 

12.3.2 Survival times

Podemos calcular qual é a probabilidade de sobrevivencia ao ano (365.25 dias).

# Criar um modelo de sobrevivência usando a função survfit
# Surv(time, status) ~ 1 indica um modelo sem variáveis preditoras (análise global)
model_surv <- survfit(Surv(time, status) ~ 1, data = lung)

# Obter um resumo do modelo de sobrevivência no tempo especificado (1 ano)
# times = 365.25 especifica que o resumo deve ser para o tempo de 1 ano
summary(model_surv, times = 365.25)

Call: survfit(formula = Surv(time, status) ~ 1, data = lung)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  365     65     121    0.409  0.0358        0.345        0.486

Ou podemos calcular qual é a mediana de sobrevivencia (quanto tempo é que 50% dos doentes tiveram o evento)

# Carregar os pacotes necessários
# library(survival)
# library(gtsummary)

# Criar um modelo de sobrevivência para o conjunto de dados 'lung'
# usando a função survfit()
model_surv <- survfit(Surv(time, status) ~ 1, data = lung)

# Utilizar a função tbl_survfit() para criar uma tabela resumindo a sobrevivência mediana
model_surv %>% 
  tbl_survfit(
    probs = 0.5,  # Definir a probabilidade para a sobrevivência mediana (50%)
    label_header = "**Sobrevivência mediana (IC 95%)**"  # Definir o cabeçalho da tabela
  )

Characteristic	Sobrevivência mediana (IC 95%)
Overall	310 (285, 363)

12.3.3 Cox regression

Podemos querer quantificar o tamanho do efeito de uma única variável, ou incluir mais de uma variável em um modelo de regressão para levar em conta os efeitos de múltiplas variáveis.

O modelo de regressão de Cox é um modelo semi-paramétrico que pode ser usado para ajustar modelos de regressão univariada e multivariada que têm resultados de sobrevivência.

Podemos ajustar modelos de regressão para dados de sobrevivência usando a função coxph() do pacote {survival}, que leva um objeto Surv() no lado esquerdo e tem uma sintaxe padrão para fórmulas de regressão em R no lado direito.

# Carregar os pacotes necessários
# library(survival)
# library(gtsummary)

# Realizar uma análise de regressão de Cox para o conjunto de dados 'lung'
# Investigando a influência do sexo na sobrevivência
cox_model <- coxph(Surv(time, status) ~ sex, data = lung)

# Criar uma tabela formatada com os resultados da regressão de Cox
cox_model|> 
  tbl_regression(exp = TRUE)  # 'exp = TRUE' exponencia os coeficientes para apresentá-los como razões de risco (hazard ratios)

Characteristic	HR	95% CI	p-value
sex	0.59	0.42, 0.82	0.001
Abbreviations: CI = Confidence Interval, HR = Hazard Ratio

13 FIM

Fundamentals of Data Science in Healthcare Fundamentals of Data Science in Healthcare 9  Módulo 7 | Aprendizagem automática com R 7  Módulo 5 | Visualização de dados 1  Bem-vindos 2  Módulo 1 | Fundamentos 3  Módulo 2 | Data Science Collaboration tools: Github 4  Módulo 2 | Quarto 5  Módulo 4 | Manipulação de dados 6  Módulo 4 | Exploração de dados 7  Módulo 5 | Visualização de dados 8  Módulo 6 | Modelos estatísticos 9  Módulo 7 | Aprendizagem automática com R 10  Module 8 | Dashboard 11  Módulo 9 | Ferramentas intermédias de DataScience 12  Módulo 10 | Capstone 8  Módulo 6 | Modelos estatísticos

Fundamentals of Data Science in Healthcare - 8  Módulo 6 | Modelos estatísticos Fundamentals of Data Science in Healthcare - 8  Módulo 6 | Modelos estatísticos Fundamentals of Data Science in Healthcare - 8  Módulo 6 | Modelos estatísticos Fundamentals of Data Science in Healthcare