Нахождение таблицы ANOVA и 95% ДИ в R

Question

Я пытаюсь запустить ANOVA в исследовании по пересадке костного мозга в Университете штата Огайо.

Вот данные с использованием кода из ( Модель пропорциональной опасности Кокса ).

time_Allo_NHL<- c(28,32,49,84,357,933,1078,1183,1560,2114,2144)
censor_Allo_NHL<- c(rep(1,5), rep(0,6))
df_Allo_NHL <- data.frame(group = "Allo NHL", 
                          time = time_Allo_NHL,
                          censor = censor_Allo_NHL,
                          Z1 = c(90,30,40,60,70,90,100,90,80,80,90),
                          Z2 = c(24,7,8,10,42,9,16,16,20,27,5))

time_Auto_NHL<- c(42,53,57,63,81,140,176,210,252,476,524,1037)
censor_Auto_NHL<- c(rep(1,7), rep(0,1), rep(1,1), rep(0,1), rep(1,1), rep(0,1))
df_Auto_NHL <- data.frame(group = "Auto NHL", 
                          time = time_Auto_NHL, 
                          censor = censor_Auto_NHL,
                          Z1 = c(80,90,30,60,50,100,80,90,90,90,90,90),
                          Z2 = c(19,17,9,13,12,11,38,16,21,24,39,84))

time_Allo_HOD<- c(2,4,72,77,79)
censor_Allo_HOD<- c(rep(1,5))
df_Allo_HOD <- data.frame(group = "Allo HOD", 
                          time = time_Allo_HOD, 
                          censor = censor_Allo_HOD,
                          Z1 = c(20,50,80,60,70),
                          Z2 = c(34,28,59,102,71))

time_Auto_HOD<- c(30,36,41,52,62,108,132,180,307,406,446,484,748,1290,1345)
censor_Auto_HOD<- c(rep(1,7), rep(0,8))
df_Auto_HOD <- data.frame(group = "Auto HOD", 
                          time = time_Auto_HOD, 
                          censor = censor_Auto_HOD,
                          Z1 = c(90,80,70,60,90,70,60,100,100,100,100,90,90,90,80),
                          Z2 = c(73,61,34,18,40,65,17,61,24,48,52,84,171,20,98))

myData <- Reduce(rbind, list(df_Allo_NHL, df_Auto_NHL, df_Allo_HOD, df_Auto_HOD)) 

Затем я использовал телеканал:

myData<- Reduce(rbind, list(df_Allo_NHL, df_Auto_NHL, df_Allo_HOD, df_Auto_HOD))
myData

library(survival)

for(i in 1:43){
  if (myData$group[i]=="Auto NHL")
    myData$Z1[i]<-1
  else myData$Z1[i]<-0
}

for(i in 1:43){
  if (myData$group[i]=="Allo HOD")
    myData$Z2[i]<-1
  else myData$Z2[i]<-0
}

for(i in 1:43){
  if (myData$group[i]=="Auto HOD")
    myData$Z3[i]<-1
  else myData$Z3[i]<-0
}

myData

Coxfit<-coxph(Surv(time,censor)~Z1+Z2+Z3, data = myData)
summary(Coxfit)

Затем, 1 - Я хочу получить ANOVA. Как я могу это сделать, включая (DF, оценки параметров, стандартную ошибку, хи-квадрат и p-значение) в результатах?

Значения должны быть следующими:

2- Как я могу получить ANOVA типа заболевания путем взаимодействия трансплантата?

Взаимодействие было получено следующим образом из ( Модель пропорционального риска Кокса-взаимодействие ):

library(tidyr)
myData <- separate(myData, col=group, into=c("disease","transpl"))
head(myData)



Coxfit.W<-coxph(Surv(time,censor)~transplant_type*disease_type, data = myData)
summary(Coxfit.W)

Значения должны быть следующими:

3- Как получить точечные оценки и 95% ДИ для Относительный риск смерти пациента с НХЛ-аутотрансплантацией по сравнению с пациентом с трансплантацией НХЛ Алло?

Результаты должны быть: 1,94 при 95% ДИ (0,64, 5,87)

Answer 1

Для 1 параметры находятся в coefficients(Coxfit), а значение anova можно получить с помощью anova(Coxfit), вы можете попробовать:

aovres = anova(Coxfit)
cbind(paramter = coefficients(Coxfit),aovres[-1,c(3,2,4)])

    paramter Df      Chisq  Pr(>|Chi|)
Z1 0.6422532  1 0.38223478 0.536409640
Z2 1.8212306  1 7.31476152 0.006839047
Z3 0.1584653  1 0.07297854 0.787048385

Для 2 вы можете изменить коэффициенты:

splitData <- separate(myData, col=group, into=c("disease","transpl"))
splitData$transpl = factor(splitData$transpl,levels=c("NHL","HOD"))
Coxfit.W<-coxph(Surv(time,censor)~transpl*disease, data = splitData)

Если вы проверяете сводку, она дает вам необходимые коэффициенты (как в предполагаемом соотношении шансов) в anova:

summary(Coxfit.W)

                           coef exp(coef) se(coef)      z Pr(>|z|)   
transplHOD              1.82123   6.17946  0.67473  2.699  0.00695 **
diseaseAuto             0.64225   1.90076  0.56415  1.138  0.25493   
transplHOD:diseaseAuto -2.30502   0.09976  0.84938 -2.714  0.00665 **

cbind(paramter = coefficients(Coxfit.W),anova(Coxfit.W)[-1,c(3,2,4)])

Для 3. это может быть проще использовать ваш оригинальный data.frame с «group», поэтому вы устанавливаете опорный уровень как «NHL allo», и отношение всех шансов будет относительно этого:

myData$group = factor(myData$group,levels=c("Allo NHL","Auto NHL","Allo HOD","Auto HOD"))

summary(coxph(Surv(time,censor)~group, data = myData))
Call:
coxph(formula = Surv(time, censor) ~ group, data = myData)

  n= 43, number of events= 26 

                coef exp(coef) se(coef)     z Pr(>|z|)   
groupAuto NHL 0.6423    1.9008   0.5641 1.138  0.25493   
groupAllo HOD 1.8212    6.1795   0.6747 2.699  0.00695 **
groupAuto HOD 0.1585    1.1717   0.5889 0.269  0.78785   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
groupAuto NHL     1.901     0.5261    0.6291     5.743
groupAllo HOD     6.179     0.1618    1.6467    23.189
groupAuto HOD     1.172     0.8535    0.3695     3.716

У вас есть отношение шансов в «groupAuto» НХЛ ".. (не очень понятно, почему он отличается от того, что вы ожидаете)