###############################################################
# A simulated data: mixture of two normals
###############################################################
rm(list=ls())
library(survival)
library(spBayesSurv)
library(coda)
library(MASS)
## True parameters
betaT = cbind(c(3.5, 0.5), c(2.5, -1));
wT = c(0.4, 0.6);
sig2T = c(1^2, 0.5^2);
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);
X = cbind(rep(1,ntot), x);
p = ncol(X); # number of covariates + 1
## Generate transformed log of survival times
z = mvrnorm(1, rep(0, ntot), corT);
## The pdf of Ti:
fi = function(y, xi, w=wT){
nw = length(w);
ny = length(y);
res = matrix(0, ny, nw);
Xi = c(1,xi);
for (k in 1:nw){
res[,k] = w[k]*dnorm(y, sum(Xi*betaT[,k]), sqrt(sig2T[k]) )
}
apply(res, 1, sum)
}
## true plot
xx = seq(-2, 7, 0.01)
#plot(xx, fi(xx, -1), "l", lwd=2, col=2)
#lines(xx, fi(xx, 1), "l", lwd=2, col=3)
## The CDF of Ti:
Fi = function(y, xi, w=wT){
nw = length(w);
ny = length(y);
res = matrix(0, ny, nw);
Xi = c(1,xi);
for (k in 1:nw){
res[,k] = w[k]*pnorm(y, sum(Xi*betaT[,k]), sqrt(sig2T[k]) )
}
apply(res, 1, sum)
}
## The inverse for CDF of Ti
Finvsingle = function(u, xi) {
res = uniroot(function (x) Fi(x, xi)-u, lower=-500, upper=500);
res$root
}
Finv = function(u, xi) {sapply(u, Finvsingle, xi)};
## Generate log of survival times y
u = pnorm(z);
y = rep(0, ntot);
for (i in 1:ntot){
y[i] = Finv(u[i], x[i]);
}
#plot(x,y);
yTrue = y;
## Censoring scheme
Centime = runif(ntot, 3.5,5);
Centime = 10000;
delta = (y<=Centime) +0 ;
sum(delta)/ntot;
cen = which(delta==0);
y[cen] = Centime[cen];
## make a data frame
dtotal = data.frame(s1=s1, s2=s2, y=y, x=x, delta=delta, yTrue=yTrue);
## Hold out npred=30 for prediction purpose
predindex = sample(1:ntot, npred);
dpred = dtotal[predindex,];
dtrain = dtotal[-predindex,];
# rename the variables
d = dtrain;n=nrow(d); n;
s = cbind(d$s1, d$s2);
y = d$y;
x = d$x;
delta =d$delta;
# Prediction settings
xpred = dpred$x;
s0 = cbind( dpred$s1, dpred$s2 );
prediction = list(spred=s0, xpred=xpred);
###############################################################
# ANOVA DDP
###############################################################
# Prior information
prior = list(N = 10,
a0 = 2, b0 = 2);
# current state values
state <- NULL
# MCMC parameters
nburn <- 500
nsave <- 500
nskip <- 0
ndisplay <- 500
mcmc <- list(nburn=nburn,
nsave=nsave,
nskip=nskip,
ndisplay=ndisplay)
# Note larger nburn, nsave and nskip should be used in practice.
# Fit model
ptm<-proc.time()
res1 = anovaDDP( y = y,
delta =delta,
x = x,
prediction=prediction,
prior=prior,
mcmc=mcmc,
state=state);
systime1=proc.time()-ptm; systime1;
par(mfrow = c(2,2))
w.save = res1$w;
Kindex = which.max(rowMeans(w.save));
traceplot(mcmc(w.save[Kindex,]), main="w")
sig2.save = res1$sigma2;
traceplot(mcmc(sig2.save[Kindex,]), main="sig2")
beta.save = res1$beta;
traceplot(mcmc(beta.save[2,Kindex,]), main="beta")
alpha.save = res1$alpha;
traceplot(mcmc(alpha.save), main="alpha")
## LPML
#cpo=CPOanovaDDP(y, delta, X, beta.save, sig2.save, w.save )$cpo;
LPML1 = sum(log(res1$cpo)); LPML1;
## MSPE
mean((dpred$yTrue-apply(res1$Ypred, 1, median))^2);
## number of non-negligible components
quantile(colSums(res1$w>0.05))
## plots
par(mfrow = c(2,2))
xnew = c(-1, 1)
xpred = cbind(xnew);
nxpred = nrow(xpred);
ygrid = seq(0,6.0,0.2); tgrid = exp(ygrid);
ngrid = length(ygrid);
estimates = GetCurves(res1, xpred, ygrid, CI=c(0.05, 0.95));
fhat = estimates$fhat;
Shat = estimates$Shat;
## density in y
plot(ygrid, fi(ygrid, xnew[1]), "l", lwd=2, ylim=c(0, 0.8),
xlim=c(0,6), main="density in y")
for(i in 1:nxpred){
lines(ygrid, fi(ygrid, xnew[i]), lwd=2)
lines(ygrid, fhat[,i], lty=2, lwd=2, col=4);
lines(ygrid, estimates$fhatup[,i], lty=2, lwd=1, col=4);
lines(ygrid, estimates$fhatlow[,i], lty=2, lwd=1, col=4);
}
## survival in y
plot(ygrid, 1-Fi(ygrid, xnew[1]), "l", lwd=2, ylim=c(0, 1),
xlim=c(0,6), main="survival in y")
for(i in 1:nxpred){
lines(ygrid, 1-Fi(ygrid, xnew[i]), lwd=2)
lines(ygrid, Shat[,i], lty=2, lwd=2, col=4);
lines(ygrid, estimates$Shatup[,i], lty=2, lwd=1, col=4);
lines(ygrid, estimates$Shatlow[,i], lty=2, lwd=1, col=4);
}
## density in t
plot(tgrid, fi(ygrid, xnew[1])/tgrid, "l", lwd=2, ylim=c(0, 0.15),
xlim=c(0,100), main="density in t")
for(i in 1:nxpred){
lines(tgrid, fi(ygrid, xnew[i])/tgrid, lwd=2)
lines(tgrid, fhat[,i]/tgrid, lty=2, lwd=2, col=4);
lines(tgrid, estimates$fhatup[,i]/tgrid, lty=2, lwd=1, col=4);
lines(tgrid, estimates$fhatlow[,i]/tgrid, lty=2, lwd=1, col=4);
}
## survival in t
plot(tgrid, 1-Fi(ygrid, xnew[1]), "l", lwd=2, ylim=c(0, 1),
xlim=c(0,100), main="survival in t")
for(i in 1:nxpred){
lines(tgrid, 1-Fi(ygrid, xnew[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);
}
Run the code above in your browser using DataLab