# R/calibratemcmc.R In spatsurv: Bayesian Spatial Survival Analysis with Parametric Proportional Hazards Models

```##' QuadApprox function
##'
##' A function to compute the second derivative of a function (of several real variables) using a quadratic approximation  on a
##' grid of points defined by the list argRanges. Also returns the local maximum.
##'
##' @param fun a function
##' @param npts integer number of points in each direction
##' @param argRanges a list of ranges on which to construct the grid for each parameter
##' @param plot whether to plot the quadratic approximation of the posterior (for two-dimensional parameters only)
##' @param ... other arguments to be passed to fun
##' @return a 2 by 2 matrix containing the curvature at the maximum and the (x,y) value at which the maximum occurs
##' @export

npar <- length(argRanges)
vals <- lapply(argRanges,function(x){seq(x[1],x[2],length.out=npts)})
parn <- paste("x",1:npar,sep="")
parnames <- parn # later this will be used in the call to lm
paridx <- as.list(rep(NA,1+2*npar+choose(npar,2))) # this will be used later to find the function maximum via a set of simultaneous equations (intercept, single, squared and mixed terms)
paridx[(1+1:npar)] <- 1:npar
gr <- expand.grid(vals)
gr2 <- gr^2
parnames <- c(parnames,paste("x",1:npar,".2",sep=""))
paridx[(npar+2):(2*npar+1)] <- 1:npar

if(npar>1){
grcross <- matrix(NA,nrow(gr),choose(npar,2))
ct <- 1
for(i in 1:(npar-1)){
for(j in (i+1):npar){
grcross[,ct] <- gr[,i]*gr[,j]
parnames <- c(parnames,paste(parn[i],parn[j],collapse="",sep=""))
paridx[[2*npar+1+ct]] <- c(i,j)
ct <- ct + 1
}
}
partype <- c("intercept",rep("single",npar),rep("squared",npar),rep("mixed",choose(npar,2)))
dataf <- cbind(gr,gr2,grcross)
}
else{
partype <- c("intercept",rep("single",npar),rep("squared",npar))
dataf <- cbind(gr,gr2)
}

names(dataf) <- parnames

cat("Constructing quadratic approximation to posterior (this can take some time) ...\n")
if(npar>1){
dataf\$funvals <- apply(gr,1,function(params){fun(params,...)})
}
else{
dataf\$funvals <- sapply(gr[[1]],function(params){fun(params,...)})
}
cat("Done.\n")

#browser()

if(plot){
if(npar==1){
plot(vals[[1]],dataf\$funvals,main="Function")
}
if(npar==2){
image.plot(vals[[1]],vals[[2]],matrix(dataf\$funvals,npts,npts),main="Function")
}
}

form <- paste("funvals ~",paste(parnames,collapse=" + "))

#browser()
mod <- lm(form,data=dataf)
co <- coefficients(mod)

if(plot){
if(npar==1){
lines(vals[[1]],fitted(mod),type="l",main="Function")
}
if(npar==2){
}
}

# now construct matrix of second derivatives
if(npar>1){
sigmainv <- matrix(NA,npar,npar)
diag(sigmainv) <- 2 * co[which(partype=="squared")] # first the diagonal elements
idx <- which(partype=="mixed") # now the off diagonals
ct <- 1
for(i in 1:(npar-1)){
for(j in (i+1):npar){
sigmainv[i,j] <- co[idx[ct]]
sigmainv[j,i] <- co[idx[ct]]
ct <- ct + 1
}
}
}
else{
sigmainv <- 2 * co[which(partype=="squared")]
}

# lastly, create a system of simultaneous equations, Ax = b, which when solved gives the maximum
b <- (-1) * matrix(co[which(partype=="single")],npar,1)
if(npar>1){
A <- matrix(NA,npar,npar)
diag(A) <- 2 * co[which(partype=="squared")]
for(i in 1:(npar-1)){
for(j in (i+1):npar){
tst <- sapply(paridx,function(x){any(x==i)&any(x==j)})
idx <- which(tst)
A[i,j] <- co[idx]
A[j,i] <- co[idx]
}
}
etaest <- as.vector(solve(A)%*%b) # now solve the system of simultaneous equations to get an initial guess for eta
}
else{
etaest <- as.vector(b/(2 * co[which(partype=="squared")]))
}

if(npar>1){
sigmainv <- fixmatrix(sigmainv)
}

return(list(max=etaest,curvature=sigmainv,mod=mod))
}

##' fixmatrix function
##'
##' !! THIS FUNCTION IS NOT INTENDED FOR GENERAL USE !!
##'
##' A function to fix up an estimated covariance matrix using a VERY ad-hoc method.
##'
##' @param mat a matrix
##' @return the fixed matrix
##' @export

fixmatrix <- function(mat){

mat <- (-1)*mat # since mat is curvature, the negative *should* have positive eigenvalues
ev <- eigen(mat)\$values
if(all(ev>0)){
return((-1)*mat)
}
else if(all(ev<0)){
stop("Estimated covariance matrix for eta has all negative eigenvalues")
}
else{
warning("Something is wrong with the estimated covariance matrix, fixing this using a totally ad-hoc method. This will not affect ergodicity, merely the efficiency of the chain.",immediate.=TRUE)
cat("Fixing non positive definite covariance matrix for eta ...\n")

diag(mat) <- abs(diag(mat)) # hmmmm ....

if(all(dim(mat)==2)){
fun <- function(x){
tmp <- mat
tmp[1,2] <- tmp[2,1] <- mat[1,2] / x
posev <- abs(ev)
ev1 <- eigen(tmp)\$values
if(!all(ev1>0)){
return(.Machine\$double.xmax)
}
else{
df1 <- (posev[1]-ev1[1])/posev[1]
df2 <- (posev[2]-ev1[2])/posev[2]
return(df1^2+df2^2)
}
}
op <- suppressWarnings(try(optimise(fun,interval=c(0,10))))
if(inherits(op,"try-error")){
stop("Failed to fix negative definite matix")
}
ans <- mat
ans[1,2] <- ans[2,1] <- mat[1,2] / op\$minimum

}
else{
#browser()
fun1 <- function(pars){
tmp <- mat
tmp[lower.tri(tmp)] <- tmp[lower.tri(tmp)] / pars
tmp[upper.tri(tmp)] <- tmp[upper.tri(tmp)] / pars
posev <- abs(ev)
ev1 <- eigen(tmp)\$values
if(!all(ev1>0)){
return(.Machine\$double.xmax)
}
else{
dff <- sum(((posev-ev1)/posev)^2)
return(dff)
}
}
op <- suppressWarnings(try(optim(par=rep(1,ncol(mat)),fn=fun1)))
if(inherits(op,"try-error")){
stop("Failed to fix negative definite matix")
}

ans <- mat
ans[lower.tri(ans)] <- ans[lower.tri(ans)] / op\$par
ans[upper.tri(ans)] <- ans[upper.tri(ans)] / op\$par
}

ct <- nrow(ans)
if(!all(eigen(ans)\$values>0)){
while(!all(eigen(ans)\$values>0) & ct>0){
ans[ct,1:(ct-1)] <- 0
ans[1:(ct-1),ct] <- 0
ct <- ct - 1
}
}

if(!all(eigen(ans)\$values>0)){
stop("Failed to fix negative definite matix")
}

ans <- (-1)*ans

return(ans)
}
}

##' proposalVariance function
##'
##' A function to compute an approximate scaling matrix for the MCMC algorithm. Not intended for general use.
##'
##' @param X the design matrix, containing covariate information
##' @param surv an object of class Surv
##' @param betahat an estimate of beta
##' @param omegahat an estimate of omega
##' @param Yhat an estimate of Y
##' @param priors the priors
##' @param cov.model the spatial covariance model
##' @param u a vector of pairwise distances
##' @param control a list containg various control parameters for the MCMC and post-processing routines
##' @return an estimate of eta and also an approximate scaling matrix for the MCMC
##' @export

proposalVariance <- function(X,surv,betahat,omegahat,Yhat,priors,cov.model,u,control){

n <- nrow(X)
lenbeta <- length(betahat)
lenomega <- length(omegahat)
leneta <- getleneta(cov.model)
lenY <- length(Yhat)
npars <- lenbeta + lenomega + leneta + lenY

sigma <- matrix(0,npars,npars)

# eta
logpost <- function(eta,surv,X,beta,omega,Y,priors,cov.model,u,control){
etapars <- cov.model\$itrans(eta)
sigma <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etapars),n,n)
cholsigma <- t(chol(sigma))
cholsigmainv <- solve(cholsigma)
MU <- -etapars[control\$sigmaidx]^2/2
gamma <- cholsigmainv%*%(Y-MU)

logpost <- logPosterior(surv=surv,X=X,beta=beta,omega=omega,eta=eta,gamma=gamma,priors=priors,cov.model=cov.model,u=u,control=control)

return(logpost)
}

npts <- 20
if(leneta>=3){
npts <- 10
}
rgs <- getparranges(priors=priors,leneta=leneta)

matr <- qa\$curvature
etahat <- qa\$max

# entry for eta in proposal covariance
sigma[(lenbeta+lenomega+1):(lenbeta+lenomega+leneta),(lenbeta+lenomega+1):(lenbeta+lenomega+leneta)] <- matr

#estimate of gamma
etahatpars <- cov.model\$itrans(etahat)
ssigma <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etahatpars),n,n)
cholssigma <- t(chol(ssigma))
MU <- -etahatpars[control\$sigmaidx]^2/2
gammahat <- solve(cholssigma)%*%(Yhat-MU)

hessian <- logPosterior(surv=surv,X=X,beta=betahat,omega=omegahat,eta=etahat,gamma=gammahat,priors=priors,cov.model=cov.model,u=u,control=control,hessian=TRUE)

# beta and omega
sigma[1:lenbeta,1:lenbeta] <- hessian\$hess_beta
sigma[(lenbeta+1):(lenbeta+lenomega),(lenbeta+1):(lenbeta+lenomega)] <- hessian\$hess_omega
sigma[(lenbeta+1):(lenbeta+lenomega),(1:lenbeta)] <- hessian\$hess_omega_beta
sigma[(1:lenbeta),(lenbeta+1):(lenbeta+lenomega)] <- t(hessian\$hess_omega_beta)
# gamma
diag(sigma)[(lenbeta+lenomega+leneta+1):npars] <- hessian\$hess_gamma

return(list(etahat=etahat,sigma=solve(-sigma)))
}

##' proposalVariance_gridded function
##'
##' A function to compute an approximate scaling matrix for the MCMC algorithm. Not intended for general use.
##'
##' @param X the design matrix, containing covariate information
##' @param surv an object of class Surv
##' @param betahat an estimate of beta
##' @param omegahat an estimate of omega
##' @param Yhat an estimate of Y
##' @param priors the priors
##' @param cov.model the spatial covariance model
##' @param u a vector of pairwise distances
##' @param control a list containg various control parameters for the MCMC and post-processing routines
##' @return an estimate of eta and also an approximate scaling matrix for the MCMC
##' @export

proposalVariance_gridded <- function(X,surv,betahat,omegahat,Yhat,priors,cov.model,u,control){

Ygrid <- gridY(Y=Yhat,control=control)

n <- nrow(X)
lenbeta <- length(betahat)
lenomega <- length(omegahat)
leneta <- getleneta(cov.model)
lenY <- length(Ygrid)
npars <- lenbeta + lenomega + leneta + lenY

# eta
logpost <- function(eta,surv,X,beta,omega,Ygrid,priors,cov.model,u,control){

etapars <- cov.model\$itrans(eta)
covbase <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etapars),control\$Mext,control\$Next)

rootQeigs <- sqrt(1/Re(fft(covbase)))

ymean <- -etapars[control\$sigmaidx]^2/2
gamma <- GammafromY(Ygrid,rootQeigs=rootQeigs,mu=ymean)

logpost <- logPosterior_gridded(surv=surv,X=X,beta=beta,omega=omega,eta=eta,gamma=gamma,priors=priors,cov.model=cov.model,u=u,control=control)

return(logpost)
}

npts <- 20
if(leneta>=3){
npts <- 10
}
rgs <- getparranges(priors=priors,leneta=leneta)

matr <- qa\$curvature
etahat <- qa\$max

sigma <- matrix(0,lenbeta + lenomega + leneta,lenbeta + lenomega + leneta)

# entry for eta in proposal covariance
sigma[(lenbeta+lenomega+1):(lenbeta+lenomega+leneta),(lenbeta+lenomega+1):(lenbeta+lenomega+leneta)] <- matr

#estimate of gamma
etahatpars <- cov.model\$itrans(etahat)
covbase <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etahatpars),control\$Mext,control\$Next)
rootQeigs <- sqrt(1/Re(fft(covbase)))
ymean <- -etahatpars[control\$sigmaidx]^2/2
gammahat <- GammafromY(Ygrid,rootQeigs=rootQeigs,mu=ymean)

hessian <- logPosterior_gridded(surv=surv,X=X,beta=betahat,omega=omegahat,eta=etahat,gamma=gammahat,priors=priors,cov.model=cov.model,u=u,control=control,hessian=TRUE)

# beta and omega
sigma[1:lenbeta,1:lenbeta] <- hessian\$hess_beta
sigma[(lenbeta+1):(lenbeta+lenomega),(lenbeta+1):(lenbeta+lenomega)] <- hessian\$hess_omega
sigma[(lenbeta+1):(lenbeta+lenomega),(1:lenbeta)] <- hessian\$hess_omega_beta
sigma[(1:lenbeta),(lenbeta+1):(lenbeta+lenomega)] <- t(hessian\$hess_omega_beta)

#sigma[1:(lenbeta+lenomega),1:(lenbeta+lenomega)] <- as.matrix(nearPD(sigma[1:(lenbeta+lenomega),1:(lenbeta+lenomega)])\$mat)

# gamma
hess_gam <- hessian\$hess_gamma

sigma <- (-1) * sigma # variance is inverse of observed information

matidx <- (lenbeta+lenomega+leneta+1):npars
matidx <- matrix(matidx,nrow=length(matidx),ncol=2)

#browser()
#sigma[1:11,1:11] <- 0
#sigma[1:2,1:2] <- diag(1/1e-4,2)
#sigma[3:9,3:9] <- diag(1/1e-4,7)
#sigma[10:11,10:11] <- diag(1/1e-4,2)

sigmaret <- Matrix(0,npars,npars)
sigmaret[1:(lenbeta+lenomega+leneta),1:(lenbeta+lenomega+leneta)] <- solve(sigma)
sigmaret[matidx] <- -1/hess_gam

return(list(etahat=etahat,sigma=sigmaret))
}

##' proposalVariance_polygonal function
##'
##' A function to compute an approximate scaling matrix for the MCMC algorithm. Not intended for general use.
##'
##' @param X the design matrix, containing covariate information
##' @param surv an object of class Surv
##' @param betahat an estimate of beta
##' @param omegahat an estimate of omega
##' @param Yhat an estimate of Y
##' @param priors the priors
##' @param cov.model the spatial covariance model
##' @param u a vector of pairwise distances
##' @param control a list containg various control parameters for the MCMC and post-processing routines
##' @return an estimate of eta and also an approximate scaling matrix for the MCMC
##' @export

proposalVariance_polygonal <- function(X,surv,betahat,omegahat,Yhat,priors,cov.model,u,control){

Ygrid <- gridY_polygonal(Y=Yhat,control=control)

lenbeta <- length(betahat)
lenomega <- length(omegahat)
leneta <- getleneta(cov.model)
lenY <- length(Ygrid)
npars <- lenbeta + lenomega + leneta + lenY

sigma <- matrix(0,npars,npars)

# eta
logpost <- function(eta,surv,X,beta,omega,Ygrid,priors,cov.model,u,control){

etapars <- cov.model\$itrans(eta)
sigma <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etapars),control\$n,control\$n)
cholsigma <- t(chol(sigma))
cholsigmainv <- solve(cholsigma)
MU <- -etapars[control\$sigmaidx]^2/2
gamma <- cholsigmainv%*%(Ygrid-MU)

logpost <- logPosterior_polygonal(surv=surv,X=X,beta=beta,omega=omega,eta=eta,gamma=gamma,priors=priors,cov.model=cov.model,u=u,control=control)

return(logpost)
}

npts <- 20
if(leneta>=3){
npts <- 10
}
rgs <- getparranges(priors=priors,leneta=leneta)

matr <- qa\$curvature
etahat <- qa\$max

# entry for eta in proposal covariance
sigma[(lenbeta+lenomega+1):(lenbeta+lenomega+leneta),(lenbeta+lenomega+1):(lenbeta+lenomega+leneta)] <- matr

#estimate of gamma
etahatpars <- cov.model\$itrans(etahat)
ssigma <- matrix(EvalCov(cov.model=cov.model,u=u,parameters=etahatpars),control\$n,control\$n)
cholssigma <- t(chol(ssigma))
MU <- -etahatpars[control\$sigmaidx]^2/2
gammahat <- solve(cholssigma)%*%(Ygrid-MU)

hessian <- logPosterior_polygonal(surv=surv,X=X,beta=betahat,omega=omegahat,eta=etahat,gamma=gammahat,priors=priors,cov.model=cov.model,u=u,control=control,hessian=TRUE)

# if(any(eigen(hessian\$hess_beta)\$values<0)){
#     cat("Problem with calibrating beta component, fixing using ad-hoc method...\n")
#     hessian\$hess_beta <- diag(abs(hessian\$hess_beta))
#     hessian\$hess_omega_beta <- 0
# }
# if(any(eigen(hessian\$hess_omega)\$values<0)){
#     cat("Problem with calibrating omega component, fixing using ad-hoc method ...\n")
#     hessian\$hess_omega <- diag(abs(hessian\$hess_omega))
#     hessian\$hess_omega_beta <- 0
# }

# beta and omega
sigma[1:lenbeta,1:lenbeta] <- hessian\$hess_beta
sigma[(lenbeta+1):(lenbeta+lenomega),(lenbeta+1):(lenbeta+lenomega)] <- hessian\$hess_omega
sigma[(lenbeta+1):(lenbeta+lenomega),(1:lenbeta)] <- hessian\$hess_omega_beta
sigma[(1:lenbeta),(lenbeta+1):(lenbeta+lenomega)] <- t(hessian\$hess_omega_beta)
# gamma
diag(sigma)[(lenbeta+lenomega+leneta+1):npars] <- hessian\$hess_gamma

#browser()

return(list(etahat=etahat,sigma=solve(-sigma)))
}

##' proposalVariance_SPDE function
##'
##' A function to compute an approximate scaling matrix for the MCMC algorithm. Not intended for general use.
##'
##' @param X the design matrix, containing covariate information
##' @param surv an object of class Surv
##' @param betahat an estimate of beta
##' @param omegahat an estimate of omega
##' @param Yhat an estimate of Y
##' @param priors the priors
##' @param cov.model the spatial covariance model
##' @param u a vector of pairwise distances
##' @param control a list containg various control parameters for the MCMC and post-processing routines
##' @return an estimate of eta and also an approximate scaling matrix for the MCMC
##' @export

proposalVariance_SPDE <- function(X,surv,betahat,omegahat,Yhat,priors,cov.model,u,control){

Ygrid <- gridY_polygonal(Y=Yhat,control=control)

lenbeta <- length(betahat)
lenomega <- length(omegahat)
leneta <- getleneta(cov.model)
lenY <- length(Ygrid)
npars <- lenbeta + lenomega + leneta + lenY

sigma <- matrix(0,npars,npars)

# eta
logpost <- function(eta,surv,X,beta,omega,Ygrid,priors,cov.model,u,control){

etapars <- cov.model\$itrans(eta)
prec <- (1/etapars[1])*control\$precmat(SPDEprec(etapars[2],cov.model\$order))
U <- Matrix::chol(prec)

if(cov.model\$order>1){
margvar <- etapars[1]/(4*pi*(cov.model\$order-1)*(etapars[2]-4)^(cov.model\$order-1)) # marginal var of Y = psi*variance_of_GMRF
}
else{
margvar <- etapars[1]*etapars[2]/(4*pi)
}
MU <- rep(-margvar/2,control\$n)
gamma <- GammaFromY_SPDE(Ygrid,U=U,mu=MU)

# y <- YFromGamma_SPDE(gamma=gamma,U=U,mu=MU)
# hist(Ygrid-y)
# browser()

logpost <- logPosterior_SPDE(surv=surv,X=X,beta=beta,omega=omega,eta=eta,gamma=gamma,priors=priors,cov.model=cov.model,u=u,control=control)

return(logpost)
}

npts <- 20
if(leneta>=3){
npts <- 10
}
rgs <- getparranges(priors=priors,leneta=leneta)

matr <- qa\$curvature
etahat <- qa\$max

# entry for eta in proposal covariance
sigma[(lenbeta+lenomega+1):(lenbeta+lenomega+leneta),(lenbeta+lenomega+1):(lenbeta+lenomega+leneta)] <- matr

#estimate of gamma
etahatpars <- cov.model\$itrans(etahat)
prec <- (1/etahatpars[1])*control\$precmat(SPDEprec(etahatpars[2],cov.model\$order))
# cholprec <- t(chol(prec))
# margvar <- etahatpars[1]/(4*pi*(cov.model\$order-1)*(etahatpars[2]-4)^(cov.model\$order-1)) # marginal var of Y = psi*variance_of_GMRF
# MU <- rep(-margvar/2,control\$n)
# gammahat <- GammaFromY_SPDE(Ygrid,P=prec,L=cholprec,mu=MU)
U <- chol(prec)
if(cov.model\$order>1){
margvar <- etahatpars[1]/(4*pi*(cov.model\$order-1)*(etahatpars[2]-4)^(cov.model\$order-1)) # marginal var of Y = psi*variance_of_GMRF
}
else{
margvar <- etahatpars[1]*etahatpars[2]/(4*pi)
}
MU <- rep(-margvar/2,control\$n)
gammahat <- GammaFromY_SPDE(Ygrid,U=U,mu=MU)

hessian <- logPosterior_SPDE(surv=surv,X=X,beta=betahat,omega=omegahat,eta=etahat,gamma=gammahat,priors=priors,cov.model=cov.model,u=u,control=control,hessian=TRUE)

# beta and omega
sigma[1:lenbeta,1:lenbeta] <- hessian\$hess_beta
sigma[(lenbeta+1):(lenbeta+lenomega),(lenbeta+1):(lenbeta+lenomega)] <- hessian\$hess_omega
sigma[(lenbeta+1):(lenbeta+lenomega),(1:lenbeta)] <- hessian\$hess_omega_beta
sigma[(1:lenbeta),(lenbeta+1):(lenbeta+lenomega)] <- t(hessian\$hess_omega_beta)
# gamma
diag(sigma)[(lenbeta+lenomega+leneta+1):npars] <- hessian\$hess_gamma

return(list(etahat=etahat,sigma=solve(-sigma)))
}

##' estimateY function
##'
##' A function to get an initial estimate of Y, to be used in calibrating the MCMC. Not for general use
##'
##' @param X the design matrix containing covariate information
##' @param betahat an estimate of beta
##' @param omegahat an estimate of omega
##' @param surv an object of class Surv
##' @param control a list containg various control parameters for the MCMC and post-processing routines
##' @return an estimate of Y, to be used in calibrating the MCMC
##' @export

estimateY <- function(X,betahat,omegahat,surv,control){

omega <- control\$omegaitrans(omegahat) # this is omega on the correct scale

tsubs <- guess_t(surv)
#browser()

Y <- -X%*%betahat - log(haz\$H(tsubs)) # greedy estimate of Y (maximise individual contributions to log-likelihood) ... note log(delta) is now omitted

return(Y)
}

##' guess_t function
##'
##' A function to get an initial guess of the failure time t, to be used in calibrating the MCMC. Not for general use
##'
##' @param surv an object of class Surv
##' @return a guess at the failure times
##' @export

guess_t <- function(surv){

n <- nrow(surv)

censoringtype <- attr(surv,"type")

if(censoringtype=="left" | censoringtype=="right"){
notcensored <- surv[,"status"]==1
}
else{
rightcensored <- surv[,"status"] == 0
notcensored <- surv[,"status"] == 1
leftcensored <- surv[,"status"] == 2
intervalcensored <- surv[,"status"] == 3
}

# setup function J=exp(X%*%beta + Y)*H_0(t)
if(censoringtype=="left" | censoringtype=="right"){
tsubs <- surv[,"time"]
}
else{ # else interval censored
tsubs <- surv[,"time1"]
}

for(i in 1:n){

if(notcensored[i]){
next
}

if(censoringtype=="left" | censoringtype=="right"){
if(censoringtype=="right"){
tpot <- tsubs[notcensored][tsubs[notcensored]>tsubs[i]] # potential t
}
else{ # censoringtype=="left"
tpot <- tsubs[notcensored][tsubs[notcensored]<tsubs[i]] # potential t
}
}
else{
if(rightcensored[i]){
tpot <- tsubs[notcensored][tsubs[notcensored]>tsubs[i]] # potential t
}
if(leftcensored[i]){
tpot <- tsubs[notcensored][tsubs[notcensored]<tsubs[i]] # potential t
}
if(intervalcensored[i]){
tpot <- surv[,"time1"][i] + 0.5*(surv[,"time2"][i]-surv[,"time1"][i]) # mid point of interval
}
}

if(length(tpot)==0){
next # leave tsubs[i] alone
}
else{
if(length(tpot)==1){
tpot <- c(tpot,tpot)
}
tsubs[i] <- sample(tpot,1) # ignoring covariates, sample from empirical distribution of times exceeding (right censored), or less than (left censored) the observed time
}

}

return(tsubs)

}
```

## Try the spatsurv package in your browser

Any scripts or data that you put into this service are public.

spatsurv documentation built on May 22, 2018, 5:09 p.m.