## Not run:
# ###############################################################
# # A simulated data: GAFT spatial frailty model
# ###############################################################
# library(survival)
# library(spBayesSurv)
# library(coda)
# library(mvtnorm)
#
# ## True densities
# set.seed(1234567890)
# Finvsingle = function(u, F) {
# res = uniroot(function (x) F(x)-u, lower=-1000, upper=1000,
# tol=.Machine$double.eps^0.5);
# res$root
# }
# Finv = function(u, F) {sapply(u, Finvsingle, F)};
# f0 = function(x) dnorm(x, 0, 0.8);
# F0 = function(x) pnorm(x, 0, 0.8);
# shift=1
# f1 = function(x) 0.5*dnorm(x, -shift, 0.5) + 0.5*dnorm(x, shift, 0.5)
# F1 = function(x) 0.5*pnorm(x, -shift, 0.5) + 0.5*pnorm(x, shift, 0.5);
# ff = function(x, xtf=0) {
# if(xtf==0) {res=f0(x)} else res=f1(x)
# res
# }
# FF = function(x, xtf=0){
# if(xtf==0) {res=F0(x)} else res=F1(x)
# res
# }
#
# # Simulation settings;
# betaT = c(-1, 1, -0.5);
# tau2T = 0.1;
# m = 50; # blocks
# mi = 5;
# mis = rep(mi, m);
# id = rep(1:m,mis);
# n = length(id); # Total number of subjects
# # Generate symmetric adjaceny matrix, W
# wi = rep(0, m)
# while(any(wi==0)){
# W = matrix(0,m,m)
# W[upper.tri(W,diag=FALSE)]<-rbinom(m*(m-1)/2,1,.1)
# W = W+t(W)
# wi = apply(W,1,sum) # No. neighbors
# }
# # Spatial effects, v
# Wstar = matrix(0, m-1, m-1);
# Dstar = diag(wi[-m]);
# for(i in 1:(m-1)){
# for(j in 1:(m-1)){
# Wstar[i,j] = W[j,i]-W[j,m]-W[m,i]-wi[m]
# }
# }
# Qstar = Dstar-Wstar;
# covT = tau2T*solve(Qstar);
# v0 = c(rmvnorm(1,sigma=covT));
# v = c(v0,-sum(v0));
# vn = rep(v, mis);
#
# # responses
# x1 = rnorm(n, 0, 1);
# x2 = rbinom(n, 1, 0.5);
# xtf = x2; ptf = 2;
# X = cbind(1,x1,x2); pce = ncol(X);
# u = runif(n, 0, 1)
# y = rep(0, n);
# for(i in 1:n) {
# if(x2[i]==1) {
# y[i] = sum(betaT*X[i,]) + vn[i] + Finv(u[i], F1)
# }else{
# y[i] = sum(betaT*X[i,]) + vn[i] + Finv(u[i], F0)
# }
# }
#
# # generate responses
# Cen = runif(n, 0.5, 1)
# delta = (exp(y)<=Cen)+0;
# sum(delta)/n
# tTrue = exp(y);
# tobs = cbind(tTrue, tTrue);
# tobs[which(delta==0),] = cbind(Cen[which(delta==0)], NA);
# dtotal = data.frame(tleft=tobs[,1], tright=tobs[,2], x1=x1,
# x2=x2, xtf=x2, ID=id, tTrue=tTrue, censor=delta);
# ## sort the data by ID
# d = dtotal[order(dtotal$ID),];
#
# # Prior information and MCMC
# fit0 <- survival::survreg(Surv(tleft, censor)~x1+x2, dist="lognormal", data=d);
# prior = list(maxL = 4, a0 = 5, b0 = 1);
# mcmc=list(nburn=100, nsave=500, nskip=1, ndisplay=100)
#
# # Fit the model
# res = frailtyGAFT(Surv(tleft, tright, type="interval2")~x1+x2+baseline(x1, x2)+
# frailtyprior(prior="car", ID), data=d, mcmc=mcmc, prior=prior,
# Proximity=W);
# summary(res);
#
#
# ####################################################################
# ## Get curves
# ####################################################################
# par(mfrow = c(2,1))
# ygrid = seq(-3,3,0.03); tgrid = exp(ygrid);
# ngrid = length(ygrid);
# xnew = c(0, 1)
# xpred = cbind(c(1,1.5), xnew);
# xtfpred = xpred;
# nxpred = nrow(xpred);
# estimates = plot(res, xpred, ygrid, xtfpred=xtfpred, CI=0.9, PLOT=F);
# fhat = estimates$fhat;
# Shat = estimates$Shat;
# ## density in y
# xlim=c(-3,3)
# shiftpred = rep(0, nxpred);
# for(i in 1:nxpred) shiftpred[i] = (sum(c(1,xpred[i,])*betaT));
#
# ## plots
# ## density in y
# plot(ygrid, ff(ygrid-shiftpred[1],xnew[1]), "l", lwd=2,
# xlim=xlim, main="density in y")
# for(i in 1:nxpred){
# lines(ygrid, ff(ygrid-shiftpred[i],xnew[i]), "l", lwd=2,
# xlim=xlim, main="density in y")
# lines(ygrid, fhat[,i], lty=2, lwd=2, col=4);
# }
# ## survival in y
# plot(ygrid, 1-FF(ygrid-shiftpred[1],xnew[1]), "l", lwd=2,
# xlim=xlim, main="survival in y")
# for(i in 1:nxpred){
# lines(ygrid, 1-FF(ygrid-shiftpred[i],xnew[i]), "l", lwd=2,
# xlim=xlim, main="survival in y")
# lines(ygrid, Shat[,i], lty=2, lwd=2, col=4);
# }
#
# ### trace plots
# par(mfrow = c(3,1))
# traceplot(mcmc(res$beta[1,]), main="beta1");
# traceplot(mcmc(res$beta[2,]), main="beta2");
# traceplot(mcmc(res$beta[3,]), main="beta3");
# ## End(Not run)
Run the code above in your browser using DataLab