Я пытаюсь понять функцию indeptCoxph в пакете spBayesSurv. Эта функция соответствует байесовской модели пропорциональных опасностей. Я немного зациклился на понимании частей кода R, а также на теории моделей Кокса.Застрял с примером кода пакета в R - имитируя данные для соответствия модели
Я работаю над примером авторов (см. Ниже). У них есть первые симулированные данные о времени выживания, и у меня возникают проблемы с их кодом для этого. Мне кажется, что сначала они моделируют времена выживания из экспоненциального распределения с CDF F (t) = 1- exp (-lambda * t) , за исключением того, что их значение для лямбда составляет exp (sum (xi * betaT)) , а не просто константа. Для моделирования данных параметру betaT присваивается фиксированное постоянное значение, которое является его истинным значением, а xi - данными предсказателя.
Вопрос 1-Это определение/форма лямбда из-за модели Cox Hazard? В этом примере авторы делают особые предположения о распределении выживаемости?
Вопрос 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) получает значение времени выживания t, которое удовлетворяет F (t) = u, где, я думаю, xi - предикторная переменная. Я действительно не понимаю, почему у вас должно быть нормальное CDF. Они генерировали «z» в виде одиночной ничьей из многомерного нормального распределения (с 3 компонентами), а u - вектор нормальных значений CDF u = pnorm (z). Опять же, не уверен, почему «u» должно быть сгенерировано таким образом - было бы действительно полезно, если бы связь между u, z, t и лямбдой могла быть выяснена. Ковариационная матрица для «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);