坚持使用 R 中的包示例代码 - 模拟数据以适应模型

Question

我想了解 spBayesSurv 包中的函数 indeptCoxph。该函数符合贝叶斯比例风险模型。我对部分 R 代码和 Cox 模型理论的理解有点卡住了。

我正在处理作者的示例（如下）。他们首先模拟了生存时间数据，但我无法按照他们的代码执行此操作。在我看来，首先他们是用 CDF 的指数分布来模拟生存时间 F(t) = 1- exp(-lambda*t) 除了它们对 lambda 的价值是 exp(sum(xi * betaT)) 而不仅仅是一个常数。为了模拟数据，参数 betaT 被赋予一个固定的常数值，即它的真实值，xi 是预测数据。

问题 1 - 这 definition/form 的 lambda 是由 Cox 风险模型引起的吗？在这个例子中，作者是否对生存分布做出了特殊假设？

问题2-我一直在理解以下生成生存时间数据的关键代码（当然它依赖于最后给出的早期代码）：

## Generate survival times t

u = pnorm(z);
t = rep(0, ntot);
for (i in 1:ntot){
t[i] = Finv(u[i], x[i]);
}
tTrue = t; #plot(x,t);

函数Finv(u,xi)获取满足F(t) = u的生存时间t的值，这里我认为xi是预测变量。我真的不明白为什么你必须来自普通的 CDF。他们从多变量正态分布（具有 3 个分量）中生成 "z" 作为单次抽取，并且 u 是正态 CDF 值 u = pnorm(z) 的向量。同样，不确定为什么必须以这种方式生成 "u" - 如果可以澄清 u,z,t 和 lambda 之间的关系，那将非常有帮助。 "z" 的协方差矩阵也是由作者从代码中的两个行向量 s1 和 s2 生成的——但如果我只是用生存时间数据拟合模型，s1、s2 的作用会令人困惑"t" 和预测变量 "x".

作者代码：

###############################################################
# A simulated data: Cox PH
###############################################################

rm(list=ls())
library(survival)
library(spBayesSurv)
library(coda)
library(MASS)
## True parameters
betaT = c(-1);
theta1 = 0.98; theta2 = 100000;
## generate coordinates:
## npred is the # of locations for prediction
n = 100; npred = 30; ntot = n + npred;
ldist = 100; wdist = 40;
s1 = runif(ntot, 0, wdist); s2 = runif(ntot, 0, ldist);
s = rbind(s1,s2); #plot(s[1,], s[2,]);
## Covariance matrix
corT = matrix(1, ntot, ntot);
for (i in 1:(ntot-1)){
for (j in (i+1):ntot){
dij = sqrt(sum( (s[,i]-s[,j])^2 ));
corT[i,j] = theta1*exp(-theta2*dij);
corT[j,i] = theta1*exp(-theta2*dij);
}
}
## Generate x
x = runif(ntot,-1.5,1.5);
## Generate transformed log of survival times
z = mvrnorm(1, rep(0, ntot), corT);
## The CDF of Ti: Lambda(t) = t;
Fi = function(t, xi){
res = 1-exp(-t*exp(sum(xi*betaT)));
res[which(t<0)] = 0;
res
}
## The pdf of Ti:
fi = function(t, xi){
res=(1-Fi(t,xi))*exp(sum(xi*betaT));
res[which(t<0)] = 0;
res
}
#integrate(function(x) fi(x, 0), -Inf, Inf)
## true plot
xx = seq(0, 10, 0.1)
#plot(xx, fi(xx, -1), "l", lwd=2, col=2)
#lines(xx, fi(xx, 1), "l", lwd=2, col=3)

## The inverse for CDF of Ti
Finvsingle = function(u, xi) {
res = uniroot(function (x) Fi(x, xi)-u, lower=0, upper=5000);
res$root
}
Finv = function(u, xi) {sapply(u, Finvsingle, xi)};

## Generate survival times t
u = pnorm(z);
t = rep(0, ntot);
for (i in 1:ntot){
t[i] = Finv(u[i], x[i]);
}
tTrue = t; #plot(x,t);

Answer 1

实际上，数据是在空间copula Cox PH模型的框架下生成的。阅读 the supplemental material of Zhou et al. (2015). As you are fitting non-spatial PH model, the data generating procedure can be sampled without the use of s1 and s2; see the new example at https://stats.stackexchange.com/questions/253368/bayesian-survival-analysis 的第 4.1 节很有帮助。

在这个新示例中，f0oft(t) 和 S0oft(t) 分别是基线生存函数。给定具有协变量 x 的受试者，Sioft(t,x) 和 fioft(t,x) 是该受试者的存活率和密度。 Finv(u,x)是Fioft(t,x)=1-Sioft(t,x)的反函数，即Finv(u,x)是Fioft(t,x)=uw.r.tt的解。

要生成生存数据，我们可以先生成协变量：

    x1 = rbinom(ntot, 1, 0.5); x2 = rnorm(ntot, 0, 1); X = cbind(x1, x2);

给定每个协变量向量 X，真实生存时间 tT 可以生成为

    u = runif(ntot);
    tT = rep(0, ntot);
    for (i in 1:ntot){
      tT[i] = Finv(u[i], X[i,]);
    }

这里背后的基本原理是，如果 T|x ~ F(t,x)，则 F(T,x) ~ Uniform(0,1)。

坚持使用 R 中的包示例代码 - 模拟数据以适应模型

Stuck with package example code in R - simulating data to fit a model

r

survival-analysis