График скорректированного риска

Question

Со следующим кадром данных:

death <- c(1,0,1,0,0,0,0,1,0,0,0,1,1,0,1,0,0,1,0,1,1,1,1,0,1,1,1,0,0,1)
age<-c (70,50,60,60,55,46,50,60,59,61,62,55,66,60,59,58,61,59,58,60,60,58,60,60,61,61,59,60,61,64)
survival_time <- c(4,29,24,29,29,29,29,19,29,29,29,3,9,29,15,29,29,11,29,5,13,20,22,29,16,21,9,29,29,15)
diabetes <- c(0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0)
bloodpressure <- c(160,120,150,140,135,110,139,140,153,129,149,163,179,129,144,119,100,115,145,150,130,120,122,129,116,171,129,126,159,150)
data <- data.frame(death, survival_time, diabetes, bloodpressure,age)

Я хотел бы создать сюжет, похожий на следующий:

Переменные, которые я получил, это артериальное давление, сердечно-сосудистая смертность, время жизни, диабет, возраст. Я хотел бы приспособиться к возрасту и диабету.

Есть ли способ сделать это в R или SPSS?

Моя первая попытка не дает результата, на который я надеялся:

cfit4a <- coxph(Surv(survival_time, death) ~ age + diabetes + bloodpressure, data=data) 
surv4a <- survfit(cfit4a) 
plot(surv4a, col=c(1,2,4), mark.time=F, xscale=365.25, xlab="blood pressure", ylab="Survival")

Как создать график, который отображает значения артериального давления на оси х? Как я могу отобразить скорректированный риск вместо выживания на оси у? (Я добавил пример данных в оригинальном посте выше)

Answer 1

Вы запрашиваете график скорректированного относительного риска для артериального давления в контексте многомерной модели, которая включает в себя не только артериальное давление, но и возраст (в виде числа) и диабет (в качестве показателя 1/0). Простой подход состоит в том, чтобы использовать predict.coxph с типом «риска» и представить кадр данных с кровяным давлением, заданным в виде диапазона значений с age и diabetes, установленными на некотором заданном уровне. Для этого прогноза важно, что означают указанные уровни, потому что контрольные значения для age и diabetes будут влиять на оценки относительного риска. Я считаю, что Терно использует средние значения ковариат для определения базовой опасности. Таким образом, вместо выбора конкретного значения для этих ковариат, вы могли бы использовать mean(data$age) и mean(data$diabetes), что должно было заставить значение bloodpressure RR равным 1 при среднем значении bloodpressure.

bp.df <- data.frame(bloodpressure=seq(120,210,by=5), age=60, diabetes=0)
pred.rr.bp <- predict(cfit4a, newdata=bp.df,type="risk", se.fit=TRUE)

plot( x = bp.df$bloodpressure, y=pred.rr.bp$fit, type="l", ylim=c(0,10), main="BP Relative Risk with age=60 and diabetes=0")
lines(x = bp.df$bloodpressure, 
          y=pred.rr.bp$fit+1.96*pred.rr.bp$se.fit,lty=3)
lines(x = bp.df$bloodpressure, 
          y=pred.rr.bp$fit-1.96*pred.rr.bp$se.fit,lty=3)

Можно также использовать rms -pkg. Требуется некоторое дополнительное изучение, чтобы справиться с его особыми требованиями, но если вам действительно нужно, чтобы строки std-error соответствовали срединному значению bp, как показывает ваш пример, тогда подпрограммы построения графика предоставляют эту «особенность». Меня это всегда немного беспокоило, однако точные значения линий +/- 1,96 * s.e. действительно не так важны.