Выполнение попарного теста и теста на тренд выживания

Я пытаюсь провести попарные тесты, чтобы определить, есть ли разница в выживаемости между парами групп.

используемые данные:

введите здесь описание изображения

time_Untreated<- c(20, 21, 23, 24, 24, 26, 26, 27, 28, 30)
censor_Untreated<- c(rep(1,10), rep(0,0))
censor_Untreated

time_Radiated<- c(26,28, 29, 29, 30, 30, 31, 31, 32, 35)
censor_Radiated<- c(rep(1,9), rep(0,1))
censor_Radiated

time_Radiated_BPA <- c(31, 32, 34, 35, 36, 38, 38, 39, 42, 42)
censor_Radiated_BPA <- c(rep(1,8), rep(0,2))
censor_Radiated_BPA

myData <- data.frame(time=c(time_Untreated, time_Radiated, time_Radiated_BPA),
                     status=c(censor_Untreated, censor_Radiated, censor_Radiated_BPA),
                     group= rep(1:3, each=10))

library(KMsurv)
library(survival)

Я пытался использовать функцию: pairwise_survdiff, но не смог построить на ней код.

Кроме того, я хочу провести тест на тенденцию, который проверит эту упорядоченную гипотезу (нелеченные животные будут иметь худшую выживаемость, облученные крысы будут иметь немного улучшенную выживаемость, а облученные крысы + BPA должны иметь лучшую выживаемость).

Вот что я сделал с выходными данными, но я не уверен, что такое значение хи-квадрат и значение p:
Это правильно?

KM.fit<-survfit(Surv(time,status)~group, conf.type="none", data=myData)
KM.fit

Call: survfit(formula = Surv(time, status) ~ group, data = myData, conf.type = "none") 

         n events median
group=1 10     10     25
group=2 10      9     30
group=3 10      8     37

myData.fit<-ten(Surv(time,status)~group, data=myData)
comp(myData.fit, p=0, q=0,scores =c(1,2,3))
            chiSq df     pChisq    
1          33.380  2 5.6436e-08 ***
n          30.255  2 2.6925e-07 ***
sqrtN      32.037  2 1.1046e-07 ***
S1         29.657  2 3.6307e-07 ***
S2         29.496  2 3.9349e-07 ***
FH_p=0_q=0 33.380  2 5.6436e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$tft
                   Q       Var      Z      pNorm    
1            16.0869    8.6116 5.4819 4.2081e-08 ***
n           364.0000 4741.0509 5.2864 1.2471e-07 ***
sqrtN        76.0224  196.2111 5.4272 5.7230e-08 ***
S1           11.1539    4.5558 5.2257 1.7351e-07 ***
S2           10.6871    4.2060 5.2110 1.8779e-07 ***
FH_p=0_q=0   16.0869    8.6116 5.4819 4.2081e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$scores
[1] 1 2 3

r survival
person Heidi    schedule 15.03.2020    source источник
comment
Я думаю, вам лучше преобразовать групповую переменную (1,2,3) в фактор, прежде чем подгонять объект ten. В противном случае функция будет считать, что это непрерывная переменная, что явно не так.   -  person Edward    schedule 15.03.2020

Ответы (1)


arrow_upward
3
arrow_downward

Пожалуйста, не забудьте указать все необходимые пакеты при публикации вопроса. Вы пропустили пакет survminer, необходимый для использования функции pairwise_survdiff.

library(survminer)

Следующий код работает с вашим набором данных, myData, поэтому я не уверен, какой код вы пробовали.

pairwise_survdiff(Surv(time,status)~group, data=myData)

    Pairwise comparisons using Log-Rank test 

data:  myData and group 

  1       2     
2 0.0011  -     
3 9.7e-06 0.0014

P value adjustment method: BH  # Bonferroni-Holm method of adjustment (default)

Таким образом, все три группы имеют существенно разную выживаемость.

Групповая переменная должна быть преобразована в фактор не только для целей маркировки кривых выживания, но многие функции моделирования будут предполагать, что переменная непрерывна, иначе это явно не так.

myData$group <- factor(myData$group, labels=c("Untreated","Radiated","Rad+BPA"))

KM.fit <- survfit(Surv(time,status)~group, data=myData)

ggsurvplot(KM.fit, data=myData,

       # Add median survival times (horizontal and vertical)
       surv.median.line = "hv", 

       # Legend placement, title, and labels
       legend = c(0.25, 0.75),
       legend.title = "Treatment Group",
       legend.labs = c("Untreated","Radiated","Rad+BPA"),
       # Add p-value and 95% confidence intervals
       pval = TRUE,
       conf.int = TRUE,

       # Add risk table
       risk.table = TRUE,
       tables.theme = theme_cleantable(),

       # Color palette
       palette = c("green4", "orange", "red"),
       ggtheme = theme_bw()
)

введите здесь описание изображения

myData.fit <- ten(Surv(time,status)~group, data=myData)
comp(myData.fit, p=0, q=0, scores=1:3) # Scores 1, 2, 3 will be the default 

            chiSq df     pChisq    
1          33.380  2 5.6436e-08 ***  # log-rank test
n          30.255  2 2.6925e-07 ***  # Gehan-Breslow generalized Wilcoxon
sqrtN      32.037  2 1.1046e-07 ***  # Tarone-Ware
S1         29.657  2 3.6307e-07 ***  # Peto-Peto’s modified survival estimate
S2         29.496  2 3.9349e-07 ***  # modified Peto-Peto (by Andersen)
FH_p=0_q=0 33.380  2 5.6436e-08 ***  # Fleming-Harrington

Каково значение хи-квадрата и p-значения? Все они! Функция дает вам шесть на выбор. Здесь они все существенны. Дополнительные сведения см. в документации по пакету.

В разделе $tst (не показан) представлены результаты теста на тенденции. Все сравнения и тесты на тенденции показывают, что существует статистически значимая разница в выживаемости крыс в трех группах. Нелеченные крысы имеют наихудшую выживаемость (медиана = 25 дней), за ними следуют облученные крысы (медиана = 30 дней) и облученные + BPA (медиана = 37 дней).

person Edward    schedule 15.03.2020
comment
Спасибо тебе за пояснение. Есть ли способ получить значения хи-квадрата для каждой попарно? - person Heidi; 15.03.2020
comment
Да. Вы просто выполняете обычный тест logrank: survdiff, но используете аргумент subset. Например. survdiff(Surv(time,status)~group, data=myData, subset=group!=1) сравнивает группу 2 и 3. Функция pairwise_survdiff делает именно это, но корректирует p-значение. Но будьте осторожны с выражением подмножества. Если вы факторизовали группу, используйте метки, а не исходные целочисленные коды. - person Edward; 15.03.2020
comment
Другими словами: survdiff(Surv(time,status)~group, data=myData, subset=group != "Untreated") будет сравнивать группы Radiated и Rad + BPA. - person Edward; 15.03.2020