The function example is conducted under the framework of spatial copula models (i.e. the function spCopulaCoxph). Briefly speaking, you just need to ignore the spred=s0 in the prediction settings, that is, prediction=list(xpred=xpred) is sufficient. To be more clear, a new example is attached at the end.
Alternatively, the newly developed function survregbayes (https://rdrr.io/cran/spBayesSurv/man/survregbayes.html) is more user-friendly to use, which fits three popular semiparametric survival models (either non-, iid-, CAR-, or GRF-frailties): proportional hazards, accelerated failure time, and proportional odds.
###############################################################
# A simulated data: Cox PH
###############################################################
rm(list=ls())
library(survival)
library(spBayesSurv)
library(coda)
library(MASS)
## True parameters
betaT = c(1,1);
n=500; npred=30; ntot=n+npred;
## Baseline Survival
f0oft = function(t) 0.5*dlnorm(t, -1, 0.5)+0.5*dlnorm(t,1,0.5);
S0oft = function(t) (0.5*plnorm(t, -1, 0.5, lower.tail=FALSE)+
0.5*plnorm(t, 1, 0.5, lower.tail=FALSE))
## The Survival function:
Sioft = function(t,x) exp( log(S0oft(t))*exp(sum(x*betaT)) ) ;
fioft = function(t,x) exp(sum(x*betaT))*f0oft(t)/S0oft(t)*Sioft(t,x);
Fioft = function(t,x) 1-Sioft(t,x);
## The inverse for Fioft
Finv = function(u, x) uniroot(function (t) Fioft(t,x)-u, lower=1e-100,
upper=1e100, extendInt ="yes", tol=1e-6)$root
## generate x
x1 = rbinom(ntot, 1, 0.5); x2 = rnorm(ntot, 0, 1); X = cbind(x1, x2);
## generate survival times
u = runif(ntot);
tT = rep(0, ntot);
for (i in 1:ntot){
tT[i] = Finv(u[i], X[i,]);
}
## right censoring
t_obs=tT
Centime = runif(ntot, 2, 6);
delta = (tT<=Centime) +0 ;
length(which(delta==0))/ntot; # censoring rate
rcen = which(delta==0);
t_obs[rcen] = Centime[rcen]; ## observed time
## make a data frame
dtotal = data.frame(t_obs=t_obs, x1=x1, x2=x2, delta=delta,
tT=tT);
## Hold out npred=30 for prediction purpose
predindex = sample(1:ntot, npred);
dpred = dtotal[predindex,];
dtrain = dtotal[-predindex,];
# Prediction settings
xpred = cbind(dpred$x1,dpred$x2);
prediction = list(xpred=xpred);
###############################################################
# Independent Cox PH
###############################################################
# MCMC parameters
nburn=1000; nsave=1000; nskip=0;
# Note larger nburn, nsave and nskip should be used in practice.
mcmc=list(nburn=nburn, nsave=nsave, nskip=nskip, ndisplay=1000);
prior = list(M=10);
state <- NULL;
# Fit the Cox PH model
res1 = indeptCoxph( y = dtrain$t_obs, delta =dtrain$delta,
x = cbind(dtrain$x1, dtrain$x2),RandomIntervals=FALSE,
prediction=prediction, prior=prior, mcmc=mcmc, state=state);
save.beta = res1$beta; row.names(save.beta)=c("x1","x2")
apply(save.beta, 1, mean); # coefficient estimates
apply(save.beta, 1, sd); # standard errors
apply(save.beta, 1, function(x) quantile(x, probs=c(0.025, 0.975))) # 95% CI
## traceplot
par(mfrow = c(2,1))
traceplot(mcmc(save.beta[1,]), main="beta1")
traceplot(mcmc(save.beta[2,]), main="beta2")
res1$ratebeta; # adaptive MH acceptance rate
## LPML
LPML1 = sum(log(res1$cpo)); LPML1;
## MSPE
mean((dpred$tT-apply(res1$Tpred, 1, median))^2);
## plots
par(mfrow = c(2,1))
x1new = c(0, 0);
x2new = c(0, 1)
xpred = cbind(x1new, x2new);
nxpred = nrow(xpred);
tgrid = seq(1e-10, 4, 0.03);
ngrid = length(tgrid);
estimates = GetCurves(res1, xpred, log(tgrid), CI=c(0.05, 0.95));
fhat = estimates$fhat;
Shat = estimates$Shat;
## density in t
plot(tgrid, fioft(tgrid, xpred[1,]), "l", lwd=2, ylim=c(0,3), main="density")
for(i in 1:nxpred){
lines(tgrid, fioft(tgrid, xpred[i,]), lwd=2)
lines(tgrid, fhat[,i], lty=2, lwd=2, col=4);
}
## survival in t
plot(tgrid, Sioft(tgrid, xpred[1,]), "l", lwd=2, ylim=c(0,1), main="survival")
for(i in 1:nxpred){
lines(tgrid, Sioft(tgrid, xpred[i,]), lwd=2)
lines(tgrid, Shat[,i], lty=2, lwd=2, col=4);
lines(tgrid, estimates$Shatup[,i], lty=2, lwd=1, col=4);
lines(tgrid, estimates$Shatlow[,i], lty=2, lwd=1, col=4);
}
Best Answer
What I ended up doing with (some success, but I'm still fiddling to improve) was use a cure rate model where:
SurvivalObsevered(t) = TreatedPortion(t) + (1-TreatedPortion(t))*SurvivalUntreated(t)
In effect what I'm doing is considering treatment to be a "cure" for untreated mortality. This allows me to model and ultimately solve for the quantity of interest which is untreated survival versus time which is not directly observed.