0
我想了解spBayesSurv包中的函數indeptCoxph。 該函數適合貝葉斯比例風險模型。我對理解R代碼的部分以及Cox模型理論有點困惑。卡住R中的包示例代碼 - 模擬數據以適合模型
我正在處理作者的示例（下文）。 他們首先模擬了生存時間數據，而我無法按照他們的代碼來做這件事。 在我看來，他們首先模擬CDF的指數分佈F（t）= 1- exp（-λ* t） 的生存時間，不同之處在於它們的lambda值爲 exp（sum（xi * betaT） ） 而不只是一個常數。爲了模擬數據，參數betaT被賦予一個固定的常數值，這是它的真值，xi是預測數據。
問題1-由於Cox Hazard模型，這是lambda的定義/形式嗎？ 在這個例子中，作者是否對生存分佈做出了特別的假設？
問題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。他們已經生成了「z」作爲多元正態分佈（含有3個分量）的單個繪圖，而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);