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


QuadApprox <- function(fun,npts,argRanges,plot=FALSE,...){

    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){
            image.plot(vals[[1]],vals[[2]],matrix(fitted(mod),npts,npts),main="Quadratic Approximation")
        }
    }

    # 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)

    qa <- QuadApprox(logpost,npts=npts,argRanges=rgs,plot=control$plotcal,surv=surv,X=X,beta=betahat,omega=omegahat,Y=Yhat,priors=priors,cov.model=cov.model,u=u,control=control)

    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)
    qa <- QuadApprox(logpost,npts=npts,argRanges=rgs,plot=control$plotcal,surv=surv,X=X,beta=betahat,omega=omegahat,Ygrid=Ygrid,priors=priors,cov.model=cov.model,u=u,control=control)

    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)
    qa <- QuadApprox(logpost,npts=npts,argRanges=rgs,plot=control$plotcal,surv=surv,X=X,beta=betahat,omega=omegahat,Ygrid=Ygrid,priors=priors,cov.model=cov.model,u=u,control=control)

    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)
    qa <- QuadApprox(logpost,npts=npts,argRanges=rgs,plot=control$plotcal,surv=surv,X=X,beta=betahat,omega=omegahat,Ygrid=Ygrid,priors=priors,cov.model=cov.model,u=u,control=control)

    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

    haz <- setupHazard(dist=control$dist,pars=omega,grad=FALSE,hess=FALSE)


    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)

}