7 Módulo 5 | Modelos estatísticos
8 SETUP
Limpar o ambiente
Definir a semente (seed)
set.seed(123)Carregar as bibliotecas
required_packages <- c("janitor",
"tidyverse",
"rio",
"data.table",
"broom",
"ggplot2",
"gt",
"gtsummary",
"survival",
"skimr"
)
for (pkg in required_packages) {
# install packages if not already present
if (!pkg %in% rownames(installed.packages())) {
install.packages(pkg)
}
# load packages to this current session
library(pkg, character.only = TRUE)
}
remove(required_packages)
remove(pkg)9 TABELAS
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 <- 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)9.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 | ||||
9.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")10 TESTES ESTATÍSTICOS
10.1 Comparar médias
10.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.
10.2 Comparar proporções
10.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.
# 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-1611 MODELOS DE REGRESSÃO
11.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.
11.1.1 Criar intuição
- 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 56 127.7607
2 45 110.5040
3 35 117.4034
4 60 129.5566
5 68 126.0226
6 70 139.2585# 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
- 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
-16.3517 -4.0741 -0.2257 3.7856 14.0862
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 98.8279 2.0553 48.09 <2e-16 ***
Age 0.5202 0.0398 13.07 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.615 on 98 degrees of freedom
Multiple R-squared: 0.6354, Adjusted R-squared: 0.6317
F-statistic: 170.8 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) 98.8 2.06 48.1 5.81e-70 94.7 103.
2 Age 0.520 0.0398 13.1 3.38e-23 0.441 0.599Podemos 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 128. 56 128. -0.199 0.0120 5.64 0.00000776 -0.0357
2 111. 45 122. -11.7 0.0111 5.52 0.0248 -2.10
3 117. 35 117. 0.368 0.0208 5.64 0.0000466 0.0662
4 130. 60 130. -0.484 0.0154 5.64 0.0000590 -0.0869
5 126. 68 134. -8.18 0.0269 5.58 0.0301 -1.48
6 139. 70 135. 4.02 0.0308 5.63 0.00838 0.726 11.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
11.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 <- clean_names(hd)Será que a idade está associada à pressão arterial? Podemos explorar esta possível associação
- Podemos fazer uma primeira exploração visual
hd |> ggplot(aes(age, resting_bp)) +
geom_point(alpha = 0.5)
11.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.622Vamos 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
11.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.
11.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.099 | 7909.143 | 7928.388 | 0.065 | 0.063 | 17.894 | 17.924model_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 | R2 (adj.) | RMSE | Sigma
-------------------------------------------------------------------------------------------------------
model2 | lm | 7907.2 (0.723) | 7907.2 (0.725) | 7921.6 (0.967) | 0.065 | 0.064 | 17.895 | 17.915
model3 | lm | 7909.1 (0.277) | 7909.1 (0.275) | 7928.4 (0.033) | 0.065 | 0.063 | 17.894 | 17.924compare_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.
11.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
11.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/
11.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")11.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.
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 011.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() 
11.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.486Ou 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) |
11.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 | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| sex | 0.59 | 0.42, 0.82 | 0.001 |
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||