JC train strep_tb
Nikita Burlov (Ebm_base)
2023-09-22
Тренируем разные варианты анализа для данных “strep_tb” в рамках обучения програмирования на R в Журнальном клубе
EDA (exploratory data analysis)
Сначала загружаем необходимые библиотеки
library(tidyverse)
library(gtsummary)Загружаем необходимые данные и смотрим общую информацию о них и первые строчки
df <- medicaldata::strep_tb
summary(df)## patient_id arm dose_strep_g dose_PAS_g gender
## Length:107 Streptomycin:55 Min. :0.000 Min. :0 F:59
## Class :character Control :52 1st Qu.:0.000 1st Qu.:0 M:48
## Mode :character Median :2.000 Median :0
## Mean :1.028 Mean :0
## 3rd Qu.:2.000 3rd Qu.:0
## Max. :2.000 Max. :0
## baseline_condition baseline_temp baseline_esr baseline_cavitation
## 1_Good:16 1_98-98.9F : 7 1_0-10 : 0 no :45
## 2_Fair:37 2_99-99.9F :25 2_11-20: 5 yes:62
## 3_Poor:54 3_100-100.9F:32 3_21-50:36
## 4_100F+ :43 4_51+ :65
## NA's : 1
##
## strep_resistance radiologic_6m rad_num
## 1_sens_0-8 :65 6_Considerable_improvement :32 Min. :1.000
## 2_mod_8-99 : 8 5_Moderate_improvement :23 1st Qu.:2.000
## 3_resist_100+:34 4_No_change : 5 Median :5.000
## 3_Moderate_deterioration :17 Mean :3.925
## 2_Considerable_deterioration:12 3rd Qu.:6.000
## 1_Death :18 Max. :6.000
## improved
## Mode :logical
## FALSE:52
## TRUE :55
##
##
## head (df)## # A tibble: 6 × 13
## patient_id arm dose_strep_g dose_PAS_g gender baseline_condition
## <chr> <fct> <dbl> <dbl> <fct> <fct>
## 1 0001 Control 0 0 M 1_Good
## 2 0002 Control 0 0 F 1_Good
## 3 0003 Control 0 0 F 1_Good
## 4 0004 Control 0 0 M 1_Good
## 5 0005 Control 0 0 F 1_Good
## 6 0006 Control 0 0 M 1_Good
## # ℹ 7 more variables: baseline_temp <fct>, baseline_esr <fct>,
## # baseline_cavitation <fct>, strep_resistance <fct>, radiologic_6m <fct>,
## # rad_num <dbl>, improved <lgl>Исходы, которые нас интересуют, находятся в переменных “radiologic_6m” (радиологическое изменение через 6 месяцев), “rad_num” (тоже самое, но в цифрах) и “improved” (есть улучшение или нет по результатам Rg). Мы можем посмотреть посмотреть основные описательные характеристики, чтобы очень упрощенно понять, что из себя представляют данные.
df |>
dplyr::select(arm, gender, baseline_condition, baseline_cavitation, baseline_temp, baseline_esr, radiologic_6m, rad_num, improved) |>
tbl_summary(missing = "no")| Characteristic | N = 1071 |
|---|---|
| arm | |
| Streptomycin | 55 (51%) |
| Control | 52 (49%) |
| gender | |
| F | 59 (55%) |
| M | 48 (45%) |
| baseline_condition | |
| 1_Good | 16 (15%) |
| 2_Fair | 37 (35%) |
| 3_Poor | 54 (50%) |
| baseline_cavitation | 62 (58%) |
| baseline_temp | |
| 1_98-98.9F | 7 (6.5%) |
| 2_99-99.9F | 25 (23%) |
| 3_100-100.9F | 32 (30%) |
| 4_100F+ | 43 (40%) |
| baseline_esr | |
| 1_0-10 | 0 (0%) |
| 2_11-20 | 5 (4.7%) |
| 3_21-50 | 36 (34%) |
| 4_51+ | 65 (61%) |
| radiologic_6m | |
| 6_Considerable_improvement | 32 (30%) |
| 5_Moderate_improvement | 23 (21%) |
| 4_No_change | 5 (4.7%) |
| 3_Moderate_deterioration | 17 (16%) |
| 2_Considerable_deterioration | 12 (11%) |
| 1_Death | 18 (17%) |
| rad_num | |
| 1 | 18 (17%) |
| 2 | 12 (11%) |
| 3 | 17 (16%) |
| 4 | 5 (4.7%) |
| 5 | 23 (21%) |
| 6 | 32 (30%) |
| improved | 55 (51%) |
| 1 n (%) | |
Как мы видим, переменные “radiologic_6m” и “rad_num” являются категориальными (точнее порядковыми), а “improved” бинарной. Также можно заметить, что есть переменная “arm”, где закодировано в какую группу был рандомизирован пациент.
Давайте уберем пока показатели исходов и посмотрим, как сбалансированы группы между собой по остальным переменным (можете отдельно почитать здесь, что они означают)
df |>
dplyr::select(arm, gender, baseline_condition, baseline_cavitation, baseline_temp, baseline_esr) |>
tbl_summary(by = arm, missing = "no")| Characteristic | Streptomycin, N = 551 | Control, N = 521 |
|---|---|---|
| gender | ||
| F | 31 (56%) | 28 (54%) |
| M | 24 (44%) | 24 (46%) |
| baseline_condition | ||
| 1_Good | 8 (15%) | 8 (15%) |
| 2_Fair | 17 (31%) | 20 (38%) |
| 3_Poor | 30 (55%) | 24 (46%) |
| baseline_cavitation | 32 (58%) | 30 (58%) |
| baseline_temp | ||
| 1_98-98.9F | 3 (5.5%) | 4 (7.7%) |
| 2_99-99.9F | 13 (24%) | 12 (23%) |
| 3_100-100.9F | 15 (27%) | 17 (33%) |
| 4_100F+ | 24 (44%) | 19 (37%) |
| baseline_esr | ||
| 1_0-10 | 0 (0%) | 0 (0%) |
| 2_11-20 | 3 (5.5%) | 2 (3.9%) |
| 3_21-50 | 16 (29%) | 20 (39%) |
| 4_51+ | 36 (65%) | 29 (57%) |
| 1 n (%) | ||
Сильных различий между группами незаметно. Значит мы можем сделать вывод, что наш EDA завершен, и можно приступать к оценке исходов между группами.
Testing hypotesis
Так как у нас качественные переменные, то наша нулевая гипотеза H0: ps = pc, где “ps” - частота исхода (любого) в группе стрептомицина, “pc” - в группе контроля
Давайте сначала посмотрим, как различаются исходы между группами на графиках.
df |>
count(rad_num = rad_num, arm = arm) %>%
group_by(arm) %>%
mutate(Freq = prop.table(n)) |>
ggplot() +
geom_col(aes(x = rad_num, fill = arm, y = Freq), position = position_dodge(width = 0.7),
width = 0.65) +
theme_test() +
scale_fill_brewer(palette = "Set1") +
theme(text = element_text(family = "Arial", size = 12))
library(ggmosaic)
df |>
ggplot() +
geom_mosaic(aes(x = product(improved, arm), fill=improved), show.legend = F) +
theme_test() +
scale_fill_brewer(palette = "Set1")
Думаю тут уже можно заметить, что различия есть, и они скорее всего окажутся значимыми. Время это проверить через стандартный тест \(X^{2}\) Пирсона.
chisq.test(df$arm, df$rad_num)##
## Pearson's Chi-squared test
##
## data: df$arm and df$rad_num
## X-squared = 26.966, df = 5, p-value = 5.791e-05chisq.test(df$arm, df$improved)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: df$arm and df$improved
## X-squared = 12.756, df = 1, p-value = 0.0003548df |>
dplyr::select(arm, rad_num) |>
tbl_summary(by = arm) |>
add_p(test = everything() ~ chisq.test)## Warning for variable 'rad_num':
## simpleWarning in stats::chisq.test(x = c(6, 5, 5, 5, 5, 6, 5, 5, 5, 5, 6, 5, 5, : Chi-squared approximation may be incorrect| Characteristic | Streptomycin, N = 551 | Control, N = 521 | p-value2 |
|---|---|---|---|
| rad_num | <0.001 | ||
| 1 | 4 (7.3%) | 14 (27%) | |
| 2 | 6 (11%) | 6 (12%) | |
| 3 | 5 (9.1%) | 12 (23%) | |
| 4 | 2 (3.6%) | 3 (5.8%) | |
| 5 | 10 (18%) | 13 (25%) | |
| 6 | 28 (51%) | 4 (7.7%) | |
| 1 n (%) | |||
| 2 Pearson’s Chi-squared test | |||
df |>
dplyr::select(arm, improved) |>
tbl_summary(by = arm) |>
add_p(test = everything() ~ chisq.test)| Characteristic | Streptomycin, N = 551 | Control, N = 521 | p-value2 |
|---|---|---|---|
| improved | 38 (69%) | 17 (33%) | <0.001 |
| 1 n (%) | |||
| 2 Pearson’s Chi-squared test | |||
Бинго! Мы получили р<0,05, значит можем отклонить нулевую гипотезу (различий нет), следовательно наблюдаем значимые различия в исходах между группами.
И на этом можно было закончить, но…
Models
Мы можем посмотреть более интересные данные, например, использовать логистическую регрессию для бинарного показателя “improved”
df$arm <- factor(df$arm, levels = c("Control", "Streptomycin"))
tbl_regression(glm(improved ~ arm, df, family = binomial(link = "logit")), exponentiate = T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| arm | |||
| Control | — | — | |
| Streptomycin | 4.60 | 2.07, 10.6 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
Получаем, что шансы получить улучшение в группе стрептомицина в 4,6 (95% ДИ 2,1 - 10,6) раз выше, чем в группе контроля, что также значимо. Но т.к. это РКИ мы можем еще и посчитать относительный риск.
tbl_regression(glm(improved ~ arm, df, family = binomial(link = "log")), exponentiate = T)| Characteristic | RR1 | 95% CI1 | p-value |
|---|---|---|---|
| arm | |||
| Control | — | — | |
| Streptomycin | 2.11 | 1.42, 3.38 | <0.001 |
| 1 RR = Relative Risk, CI = Confidence Interval | |||
Бинго! Риск получить улучшение в группе стрептомицина также значимо выше, в 2,1 (95% ДИ 1,4 - 3,4), чем в группе контроля.
Complex models
Но это мы оценивали только бинарный показатель, а у нас есть еще порядковый “rad_num”. Как оценить его? Тут мы переходим на поле сложных моделей. Одна из таких - это порядковая регрессия (ordinal regression).
library(rms)fit <- orm(rad_num ~ arm, df, x = T, y = T)
print(fit)## Logistic (Proportional Odds) Ordinal Regression Model
##
## orm(formula = rad_num ~ arm, data = df, x = T, y = T)
##
##
## Frequencies of Responses
##
## 1 2 3 4 5 6
## 18 12 17 5 23 32
##
## Model Likelihood Discrimination Rank Discrim.
## Ratio Test Indexes Indexes
## Obs 107 LR chi2 21.96 R2 0.192 rho 0.442
## Distinct Y 6 d.f. 1 R2(1,107) 0.178
## Median Y 5 Pr(> chi2) <0.0001 R2(1,102) 0.186
## max |deriv| 0.0004 Score chi2 21.44 |Pr(Y>=median)-0.5| 0.200
## Pr(> chi2) <0.0001
##
## Coef S.E. Wald Z Pr(>|Z|)
## y>=2 0.9679 0.2926 3.31 0.0009
## y>=3 0.2681 0.2623 1.02 0.3068
## y>=4 -0.5010 0.2599 -1.93 0.0539
## y>=5 -0.7173 0.2645 -2.71 0.0067
## y>=6 -1.8059 0.3148 -5.74 <0.0001
## arm=Streptomycin 1.6928 0.3751 4.51 <0.0001Здесь нам интересны коэфициенты и уровни р. Нам показаны как меняются шансы перехода из одного исхода “Y = 1” в другой “Y ≥ 2” в логарифимческой шкале. Можем их посмотреть в виде отношения шансов (OR).
data.frame(exp(fit$coefficients))## exp.fit.coefficients.
## y>=2 2.6324648
## y>=3 1.3074419
## y>=4 0.6059108
## y>=5 0.4880833
## y>=6 0.1643315
## arm=Streptomycin 5.4344360Мы получаем уже много интересных данных: например, шанс перехода в группе контроля из “Y = 1” в “Y ≥ 2” достаточно высокий, что логично, т.к. в “Y = 1” меньше случаев, чем в других значениях суммарно.
При этом также видно, как шансы постепенно снижаются перейти в статус “улучшение” в контрольной группе. А добавление стрептомицина достаточно сильно эти шансы повышают (но тут “сложная” математика, где надо подключать формулы).
Попроще это можно отразить в виде предсказания вероятностей получить конкретный уровень в разных группах.
library(MASS)fit <- polr(factor(rad_num) ~ arm, df, Hess = T)
predict(fit, data.frame(arm = c("Control", "Streptomycin")), type="p")## 1 2 3 4 5 6
## 1 0.27529088 0.15808604 0.1893226 0.04930645 0.1868585 0.1411355
## 2 0.06533106 0.05804195 0.1095712 0.04083921 0.2544634 0.4717531Где по горизонтали группы (1 - контроль, 2 - стрептомицин), по вертикали уровни исходов (как в переменной “rad_num”). Как видно вероятность получить конкретный исход в каждой группе разный, больше вероятностьт получить улучшение в группе стрептомицина, а ухудшение в группе контроля.
Эти вероятность хорошо согласуются с нашим первичным графиком (только на графике цифры чуть немного отличаются, потому что это исходные частоты), где мы смотрели разницу исходов между группами.
df |>
count(rad_num = rad_num, arm = arm) %>%
group_by(arm) %>%
mutate(Freq = prop.table(n)) |>
ggplot() +
geom_col(aes(x = rad_num, fill = arm, y = Freq), position = position_dodge(width = 0.7),
width = 0.65) +
theme_test() +
scale_fill_brewer(palette = "Set1") +
theme(text = element_text(family = "Arial", size = 12))
Summary
Вот такими разными методами можно было оценить разницу в исходах между группами. И это еще не предел, т.к. здесь мы рассмотрели только методы частотнического подхода.