如何在 R 中绘制 3-D 风险平面?
How to plot a 3-D risk plane in R?
我正在尝试在 R 的 3-D 图中绘制风险平面,以图形方式描述连续预测变量对其他一些连续预测变量与结果之间关联的影响修改。风险估计(HR,风险比)应该在 z 轴上,两个连续预测变量在 x 轴和 y 轴上,如下图所示:
为了说明我已经尝试过的内容,我将使用 survival
包中的 lung
数据集。
#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))
我的目的是根据
手动计算风险比
以年龄中位数和karnofsky评分(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
中的函数之后。
有谁知道这样的剧情怎么获得?
也许您可以按照这些思路使用一些东西。首先我们根据您自己的代码找到模型:
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")
我正在尝试在 R 的 3-D 图中绘制风险平面,以图形方式描述连续预测变量对其他一些连续预测变量与结果之间关联的影响修改。风险估计(HR,风险比)应该在 z 轴上,两个连续预测变量在 x 轴和 y 轴上,如下图所示:
为了说明我已经尝试过的内容,我将使用 survival
包中的 lung
数据集。
#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))
我的目的是根据
以年龄中位数和karnofsky评分(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
中的函数之后。
有谁知道这样的剧情怎么获得?
也许您可以按照这些思路使用一些东西。首先我们根据您自己的代码找到模型:
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")