Как построить трехмерную плоскость риска в R?

Я пытаюсь построить плоскость риска на трехмерном графике в языке R, чтобы графически изобразить изменение влияния непрерывного предсказателя на связь между некоторым другим непрерывным предсказателем и результатом. Оценки риска (HR, отношение рисков) должны быть на оси z, а две непрерывные переменные-предикторы - на осях x и y, как на графике ниже:

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

Чтобы проиллюстрировать то, что я уже пробовал, я воспользуюсь набором данных lung из пакета survival.

#install.packages("survival")
#install.packages("rgl")

library(survival); library(rgl)

#Remove missing values with listwise deletion
I1 <- is.na(lung$age) | is.na(lung$ph.karno)
lung <- lung[!I1,]

m1 <- coxph(Surv(time, status==2) ~ age*ph.karno, data = lung)
m1

z <- outer(lung$age, lung$ph.karno, FUN=function(x=lung$age, y=lung$ph.karno, model=m1){
  ref.x <- median(x)
  ref.y <- median(y)
  for(i in 1:length(x)){
    exp(summary(model)$coef[1,1]*(x[i]-ref.x)+summary(model)$coef[2,1]*(y[i]-ref.y)+
               summary(model)$coef[3,1]*(x[i]-ref.x)*(y[i]-ref.y))
  }
})

persp3d(x=lung$age, y=lung$ph.karno, z=z)

В exp(summary(model)$coef[1,1]*(x[i]-ref.x)+summary(model)$coef[2,1]*(y[i]-ref.y)+summary(model)$coef[3,1]*(x[i]-ref.x)*(y[i]-ref.y)) я стремился вручную рассчитать коэффициент опасности в соответствии с

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

со средним возрастом и оценкой Карновского (ph.karno), установленными в качестве соответствующих эталонов. Однако, когда я запускаю этот код, я сталкиваюсь со следующими двумя ошибками: Error in dim(robj) <- c(dX, dY) : attempt to set an attribute on NULL после выполнения функции в пределах outer() и Error in persp3d.default(x = lung$age, y = lung$ph.karno, z = z) : Increasing 'x' and 'y' values expected.

Кто-нибудь знает, как получить такой сюжет?


r 3d plot matrix survival-analysis
person Dion    schedule 01.08.2019    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Может быть, ты сможешь использовать что-нибудь в этом роде. Сначала мы находим модель по вашему собственному коду:

library("survival"); library("rgl")

#Remove missing values with listwise deletion
I1 <- is.na(lung$age) | is.na(lung$ph.karno)
lung <- lung[!I1,]

m1 <- coxph(Surv(time, status==2) ~ age*ph.karno, data = lung)
m1

Затем используйте функцию predict(), чтобы вычислить риск по модели. Поскольку модель включает взаимодействие, оно автоматически включается в прогноз. В качестве входных данных мы используем соответствующие интервалы значений в наблюдаемых диапазонах age и ph.karno.

age.range <- seq(min(lung$age), max(lung$age), 5)
ph.range <- seq(min(lung$ph.karno), max(lung$ph.karno), 5)

z <- outer(age.range, ph.range, FUN=function(x, y) {
  predict(m1, newdata = data.frame(age=x, ph.karno=y), type="risk")
  })

rgl::persp3d(age.range, ph.range, z, col="lightblue")
person ekstroem    schedule 05.08.2019
comment
Спасибо за ответ! Мне все еще интересно, как будет выглядеть конструкция, когда взаимодействие будет скорректировано для дополнительных переменных, например: m2 <- coxph(Surv(time, status==2) ~ age*ph.karno + sex + wt.loss, data = lung) - person Dion; 06.08.2019
comment
Пока вы рассматриваете взаимодействия только между двумя предикторами, другие переменные по существу не имеют значения (для иллюстрации взаимодействия). Они просто увеличат или уменьшат общий риск. Для более чем двух предикторов вам придется разделить данные на или график в 4D: o) - person ekstroem; 06.08.2019