R/survspatNS.R

In spatsurv: Bayesian Spatial Survival Analysis with Parametric Proportional Hazards Models

##' survspatNS function
##'
##' A function to perform maximun likelihood inference for non-spatial survival data.
##'
##' @param formula the model formula in a format compatible with the function flexsurvreg from the flexsurv package
##' @param data a SpatialPointsDataFrame object containing the survival data as one of the columns
##' @param dist choice of distribution function for baseline hazard. Current options are: exponentialHaz, weibullHaz, gompertzHaz, makehamHaz, tpowHaz
##' @param control additional control parameters, see ?inference.control
##' @return an object inheriting class 'mcmcspatsurv' for which there exist methods for printing, summarising and making inference from.
##' @seealso \link{tpowHaz}, \link{exponentialHaz}, \link{gompertzHaz}, \link{makehamHaz}, \link{weibullHaz},
##' \link{covmodel}, \link{ExponentialCovFct}, \code{SpikedExponentialCovFct},
##' \link{mcmcpars}, \link{mcmcPriors}, \link{inference.control}
##' @references
##' \enumerate{
##'     \item Benjamin M. Taylor and Barry S. Rowlingson (2017). spatsurv: An R Package for Bayesian Inference with Spatial Survival Models. Journal of Statistical Software, 77(4), 1-32, doi:10.18637/jss.v077.i04.
##' }
##' @export


survspatNS <- function( formula,
                        data,
                        dist,
                        control=inference.control()){

    control$hessian <- TRUE

    responsename <- as.character(formula[[2]])
    survivaldata <- data[[responsename]]
    checkSurvivalData(survivaldata)


    # start timing,

    start <- Sys.time()

    control$dist <- dist

    ##########
    # This chunk of code borrowed from flexsurvreg
    ##########

    call <- match.call()
    indx <- match(c("formula", "data"), names(call), nomatch = 0)
    if (indx[1] == 0){
        stop("A \"formula\" argument is required")
    }
    temp <- call[c(1, indx)]
    temp[[1]] <- as.name("model.frame")
    m <- eval(temp, parent.frame())

    Terms <- attr(m, "terms")
    X <- model.matrix(Terms, m)

    ##########
    # End of borrowed code
    ##########

    X <- X[, -1, drop = FALSE]

    info <- distinfo(dist)()

    control$omegatrans <- info$trans
    control$omegaitrans <- info$itrans
    control$omegajacobian <- info$jacobian # used in computing the derivative of the log posterior with respect to the transformed omega (since it is easier to compute with respect to omega)
    control$omegahessian <- info$hessian

    control$censoringtype <- attr(survivaldata,"type")

    if(control$censoringtype=="left" | control$censoringtype=="right"){
        control$censored <- survivaldata[,"status"]==0
        control$notcensored <- !control$censored

        control$Ctest <- any(control$censored)
        control$Utest <- any(control$notcensored)

    }
    else{
        control$rightcensored <- survivaldata[,"status"] == 0
        control$notcensored <- survivaldata[,"status"] == 1
        control$leftcensored <- survivaldata[,"status"] == 2
        control$intervalcensored <- survivaldata[,"status"] == 3

        control$Rtest <- any(control$rightcensored)
        control$Utest <- any(control$notcensored)
        control$Ltest <- any(control$leftcensored)
        control$Itest <- any(control$intervalcensored)
    }

    #######

    cat("\n","Maximum likelihood using BFGS ...","\n")
    mlmod <- maxlikparamPHsurv(surv=survivaldata,X=X,control=control)
    estim <- mlmod$par
    print(mlmod)
    cat("Done.\n")

    end <- Sys.time()

    #browser()


    retlist <- list()
    retlist$formula <- formula
    retlist$dist <- dist
    retlist$control <- control

    retlist$terms <- Terms
    retlist$mlmod <- mlmod

    # construct artificial samples from which the baseline hazard can be obtained with confidence intervals
    ch <- t(chol(solve(mlmod$hessian)))
    samp <- t(mlmod$par+ch%*%matrix(rnorm(1000*ncol(ch)),ncol(ch),1000))

    betasamp <- samp[,1:ncol(X),drop=FALSE]

    ####
    #   Back transform for output
    ####
    omegasamp <- samp[,(ncol(X)+1):ncol(samp),drop=FALSE]
    if(ncol(omegasamp)>1){
        omegasamp <- t(apply(omegasamp,1,control$omegaitrans))
    }
    else{
        omegasamp <- t(t(apply(omegasamp,1,control$omegaitrans)))
    }
    colnames(omegasamp) <- info$parnames

    colnames(betasamp) <- colnames(model.matrix(formula,data))[-1] #attr(Terms,"term.labels")
    retlist$betasamp <- betasamp
    retlist$omegasamp <- omegasamp

    retlist$Ysamp <- matrix(0,1000,nrow(X))

    #retlist$loglik <- loglik

    retlist$X <- X
    retlist$survivaldata <- survivaldata
    retlist$gridded <- control$gridded
    retlist$omegatrans <- control$omegatrans
    retlist$omegaitrans <- control$omegaitrans
    retlist$control <- control
    retlist$censoringtype <- attr(survivaldata,"type")
    retlist$time.taken <- Sys.time() - start

    cat("Time taken:",retlist$time.taken,"\n")

    class(retlist) <- c("list","mlspatsurv")

    return(retlist)
}