R：绘制个人预测答案

【问题标题】：R: Plot Individual PredictionsR：绘制个人预测
【发布时间】：2021-06-15 10:51:37
【问题描述】：

我正在使用 R 编程语言。我正在尝试遵循本教程：https://rdrr.io/cran/randomForestSRC/man/plot.competing.risk.rfsrc.html

本教程展示了如何使用“生存随机森林”算法 - 一种用于分析生存数据的算法。在本例中，使用“follic”数据集，使用生存随机森林算法分析观察经历“状态 1”与“状态 2”的即时风险（这称为“竞争风险”）。

在下面的代码中，生存随机森林模型是在毛囊数据集上使用除最后两个观察值之外的所有观察值来训练的。然后，该模型用于预测最后两次观察的危害：

#load library
library(randomForestSRC)

#load data
data(follic, package = "randomForestSRC")

#train model on all observations except the last 2 observations
follic.obj <- rfsrc(Surv(time, status) ~ ., follic[c(1:539),], nsplit = 3, ntree = 100)

#use model to predict the last two observations
f <- predict(follic.obj, follic[540:541, ])

#plot individual curves - does not work
plot.competing.risk(f)

但是，这似乎产生了最后两次观察“状态 1 与状态 2”的平均危害。

有没有办法绘制第一次观察和第二次观察的个体风险？

谢谢

编辑1：

我知道如何为这个包中的其他功能执行此操作，例如在这里，您可以一次绘制这些曲线以进行 7 次观察：

data(veteran, package = "randomForestSRC") 
plot.survival(rfsrc(Surv(time, status)~ ., veteran), cens.model = "rfsrc")

## pbc data
data(pbc, package = "randomForestSRC") 
pbc.obj <- rfsrc(Surv(days, status) ~ ., pbc)

## use subset to focus on specific individuals
plot.survival(pbc.obj, subset = c(3, 10))

这个例子似乎同时显示了 7 个观察值的预测生存曲线（加上置信区间 - 红线是平均值）。但是我仍然不知道如何为“plot.competing.risk”功能做到这一点。

EDIT2：

我认为可能有一种间接的方法来解决这个问题 - 您可以单独预测每个观察结果：

#use model to predict the last two observations individually
f1 <- predict(follic.obj, follic[540, ])
f2 <- predict(follic.obj, follic[541, ])

#plot individual curves 
plot.competing.risk(f1)
plot.competing.risk(f2)

但我希望有一种更直接的方法来做到这一点。有人知道怎么做吗？

【问题讨论】：

标签： r plot data-visualization survival-analysis

【解决方案1】：

一种可能的方法是修改函数plot.competing.risk 用于单行，并在for 循环上绘制重叠的单行，如下所示。

#use model to predict the last three observations
f <- predict(follic.obj, follic[539:541, ])

x <- f
par(mfrow = c(2, 2))
for (k in 1:3) { #k for type of plot
    for (i in 1:dim(x$chf)[1]) { #i for all individuals in x
        #cschf <- apply(x$chf, c(2, 3), mean, na.rm = TRUE) #original group mean
        cschf = x$chf[i,,] #individual values

        #cif <- apply(x$cif, c(2, 3), mean, na.rm = TRUE) #original group mean
        cif = x$cif[i,,] #individual values

        cpc <- do.call(cbind, lapply(1:ncol(cif), function(j) {
                cif[, j]/(1 - rowSums(cif[, -j, drop = FALSE]))
            }))

        if (k==1)
            {matx = cschf
            range = range(x$chf)
            }
        if (k==2)
            {matx = cif
            range = range(x$cif)
            }
        if (k==3)
            {matx = cpc
            range = c(0,1) #manually assign, for now
            }

        ylab = c("Cause-Specific CHF","Probability (%)","Probability (%)")[k]
        matplot(x$time.interest, matx, type='l', lty=1, lwd=3, col=1:2, 
            add=ifelse(i==1,F,T), ylim=range, xlab="Time", ylab=ylab) #ADD tag for overlapping individual lines
        }
    legend <- paste(c("CSCHF","CIF","CPC")[k], 1:2, "  ")
    legend("bottomright", legend = legend, col = (1:2), lty = 1, lwd = 3)
    }

【讨论】：