R/spm_continuous_quadratic_linear.R

In stpm: Stochastic Process Model for Analysis of Longitudinal and Time-to-Event Outcomes

#'Continuous multi-dimensional optimization with quadratic and linear terms
#'@references Yashin, A.I. et al (2007). Stochastic model for analysis of longitudinal data on aging 
#'and mortality. Mathematical Biosciences, 208(2), 538-551.<DOI:10.1016/j.mbs.2006.11.006>.
#'@param dat A data table.
#'@param a A starting value of the rate of adaptive response to any deviation of Y from f1(t).
#'@param f1 A starting value of the average age trajectories of the variables which process is forced to follow. 
#'@param Q Starting values of the quadratic hazard term.
#'@param f A starting value of the "optimal" value of variable which corresponds to the minimum of hazard rate at a respective time.
#'@param b A starting value of a diffusion coefficient representing a strength of the random disturbance from Wiener Process.
#'@param mu0 A starting value of the baseline hazard.
#'@param theta A starting value of the parameter theta (axe displacement of Gompertz function).
#'@param Q1 Q for linear term
#'@param verbose An indicator of verbosing output.
#'@param stopifbound Estimation stops if at least one parameter achieves lower or upper boundaries.
#'#'Check the NLopt website for a description of
#'the algorithms. Default: NLOPT_LN_NELDERMEAD
#'@param lb Lower bound of parameters under estimation.
#'@param ub Upper bound of parameters under estimation.
#'The program stops when the number of function evaluations exceeds maxeval. Default: 500.
#'@param pinv.tol A tolerance value for pseudo-inverse of matrix gamma (see Yashin, A.I. et al (2007). Stochastic model for analysis of longitudinal data on aging 
#'and mortality. Mathematical Biosciences, 208(2), 538-551.<DOI:10.1016/j.mbs.2006.11.006>.)
#'@param gomp A flag (FALSE by default). When it is set, then time-dependent exponential form of mu0 is used:
#' mu0 = mu0*exp(theta*t).
#'@param opts A list of options for \code{nloptr}.
#'Default value: \code{opt=list(algorithm="NLOPT_LN_NELDERMEAD", 
#'maxeval=100, ftol_rel=1e-8)}.
#'Please see \code{nloptr} documentation for more information.
#'@return A set of estimated parameters a, f1, Q, f, b, mu0, theta and
#'additional variable \code{limit} which indicates if any parameter 
#'achieved lower or upper boundary conditions (FALSE by default).
#'@return status Optimization status (see documentation for nloptr package).
#'@return LogLik A logarithm likelihood.
#'@return objective A value of objective function (given by nloptr).
#'@return message A message given by nloptr optimization function (see documentation for nloptr package).
#'@details \code{spm_continuous} runs much slower that discrete but more precise and can handle time 
#'intervals with different lengths.
#'@export
#'@examples
#'library(stpm)
#'set.seed(123)
#'#Reading the data:
#'data <- simdata_cont(N=2)
#'head(data)
#'#Parameters estimation:
#'pars <- spm_cont_quad_lin(dat=data,a=-0.05, f1=80, 
#'						           Q=2e-8, f=80, b=5, mu0=2e-5, Q1=1e-08)
#'pars
#'
spm_cont_quad_lin <- function(dat, 
                           a=-0.05, 
                           f1=80, 
                           Q=2e-8,
                           f=80,
                           b=5,
                           mu0=2e-5,
                           theta=0.08,
                           Q1=1e-08,
                           stopifbound=FALSE, 
                           lb=NULL, ub=NULL,
                           verbose=FALSE,
                           pinv.tol=0.01,
                           gomp=FALSE,
                           opts=list(algorithm="NLOPT_LN_NELDERMEAD", 
                                          maxeval=100, ftol_rel=1e-8)) {
  
  
  
    setlb <- function(k, params) 
    {
        # This function sets lower and upper boundaries for optim.
        # - k - number of dimensions
        # - params - initial parameters, a vector
        #
        # Lower boundaries:
        lower_bound <- c()
        # Setting boundaries for coefficients:
        # aH
        start=1; end=k^2
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) }))) 
        # f1H
        start=end+1; end=start+k-1
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) }))) 
        # QH
        start=end+1; end=start+k^2-1
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ ifelse(params[n] > 0, 0, params[n]+0.1*params[n]) })))
        # fH
        start=end+1; end=start+k-1
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })) )
        # bH
        start=end+1; end=start+k-1
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){params[n] + ifelse(params[n] > 0, -0.1*params[n], 0.1*params[n]) })) )
        # mu0
        start=end+1; end=start
        lower_bound <- c( lower_bound, params[start:end] - 0.1*params[start:end])
        # theta
        start=end+1; end=start
        lower_bound <- c( lower_bound, params[start:end] - 0.1*params[start:end])
        # Q1
        start=end+1; end=start+k^2-1
        lower_bound <- c(lower_bound, unlist(lapply(start:end, function(n){ ifelse(params[n] > 0, 0, params[n]+0.1*params[n]) })))
        
    
        return(lower_bound)
    }
  
  
  
    setub <- function(k, params) {
        # This function sets lower and upper boundaries for optim.
        # - k - number of dimensions
        # - params - initial parameters, a vector
        #
        #Upper boundaries:
        upper_bound <- c()
        # Setting boundaries for coefficients:
        # aH
        start=1; end=k^2
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n] + 0.1*params[n], 0*params[n]) })))
        # f1H
        start=end+1; end=start+k-1
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n] + 0.1*params[n], params[n]- 0.1*params[n]) })))
        # QH
        start=end+1; end=start+k^2-1
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n) {ifelse(params[n] > 0, params[n] + 0.1*params[n], params[n]-0.1*params[n] )})))
        # fH
        start=end+1; end=start+k-1
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n]+0.1*params[n], params[n]-0.1*params[n]) })))
        # bH
        start=end+1; end=start+k-1
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n){ifelse(params[n] > 0, params[n]+0.1*params[n], params[n]-0.1*params[n]) })))
        # mu0
        start=end+1; end=start
        upper_bound <- c( upper_bound, params[start:end] + 0.1*params[start:end])
        # theta
        start=end+1; end=start
        upper_bound <- c( upper_bound, params[start:end] + 0.1*params[start:end])
        # Q1
        start=end+1; end=start+k^2-1
        upper_bound <- c(upper_bound, unlist(lapply(start:end, function(n) {ifelse(params[n] > 0, params[n] + 0.1*params[n], params[n]-0.1*params[n] )})))
        
        upper_bound
    }
  
  
  ###=======For DEBUG========###
  #dat = dat
  #
  #a=matrix(c(-0.05,  0.001, 0.001, -0.05), nrow = 2, ncol = 2, byrow = T)
  #f1=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F))
  #Q=matrix(c(1e-06, 1e-7, 1e-7,  1e-06), nrow = 2, ncol = 2, byrow = T)
  #f=t(matrix(c(100, 200), nrow = 2, ncol = 1, byrow = F))
  #b=matrix(c(2, 5), nrow = 2, ncol = 1, byrow = F)
  #mu0=1e-4
  #theta=0.08
  #k=2
  #
  #stopifbound=FALSE
  #algorithm="NLOPT_LN_NELDERMEAD"
  #lb=NULL
  #ub=NULL
  #verbose=FALSE
  ###=========================###
  
  
  
  #avail_algorithms <- c("NLOPT_GN_DIRECT", "NLOPT_GN_DIRECT_L",
  #                      "NLOPT_GN_DIRECT_L_RAND", "NLOPT_GN_DIRECT_NOSCAL",
  #                      "NLOPT_GN_DIRECT_L_NOSCAL",
  #                      "NLOPT_GN_DIRECT_L_RAND_NOSCAL",
  #                      "NLOPT_GN_ORIG_DIRECT", "NLOPT_GN_ORIG_DIRECT_L",
  #                      "NLOPT_GD_STOGO", "NLOPT_GD_STOGO_RAND",
  #                      "NLOPT_LD_SLSQP", "NLOPT_LD_LBFGS_NOCEDAL",
  #                      "NLOPT_LD_LBFGS", "NLOPT_LN_PRAXIS", "NLOPT_LD_VAR1",
  #                      "NLOPT_LD_VAR2", "NLOPT_LD_TNEWTON",
  #                      "NLOPT_LD_TNEWTON_RESTART",
  #                      "NLOPT_LD_TNEWTON_PRECOND",
  #                      "NLOPT_LD_TNEWTON_PRECOND_RESTART",
  #                      "NLOPT_GN_CRS2_LM", "NLOPT_GN_MLSL", "NLOPT_GD_MLSL",
  #                      "NLOPT_GN_MLSL_LDS", "NLOPT_GD_MLSL_LDS",
  #                      "NLOPT_LD_MMA", "NLOPT_LN_COBYLA", "NLOPT_LN_NEWUOA",
  #                      "NLOPT_LN_NEWUOA_BOUND", "NLOPT_LN_NELDERMEAD",
  #                      "NLOPT_LN_SBPLX", "NLOPT_LN_AUGLAG", "NLOPT_LD_AUGLAG",
  #                      "NLOPT_LN_AUGLAG_EQ", "NLOPT_LD_AUGLAG_EQ",
  #                      "NLOPT_LN_BOBYQA", "NLOPT_GN_ISRES")
  #
  #if(!(algorithm %in% avail_algorithms)) {
  #  stop(cat("Provided algorithm", algorithm, "not in the list of available optimization methods."))
  #}
  
    dat <- as.matrix(dat[, 2:dim(dat)[2]])
  
    k <- dim(as.matrix(a))[1]
  
    final_res <- list()
  
    if(mu0 < 0) {mu0 <- 0}
  
    parameters <- c(t(a), f1, t(Q), f, b, mu0, theta, t(Q1))
    # Current results:
    results_tmp <- list(a=NULL, f1=NULL, Q=NULL, f=NULL, b=NULL, mu0=NULL, theta=NULL, Q1=NULL)
    e <- new.env()
    iteration <- 0
  
    bounds <- list()
  
    if(is.null(lb)) {
        bounds$lower_bound <- setlb(k, parameters)
    } else {
        bounds$lower_bound <- lb
    }
  
    if(is.null(ub)) {
        bounds$upper_bound <- setub(k, parameters)
    } else {
        bounds$upper_bound <- ub
    }
  
    L.prev <- NA
    maxlik <- function(par) {
        stopflag <- F
        # Reading parameters:
        start=1; end=k^2
        a <- matrix(par[start:end],ncol=k, , nrow=k, byrow=TRUE)
        results_tmp$a <<- a
        start=end+1; end=start+k-1
        f1 <- matrix(par[start:end],ncol=1, nrow=k, byrow=FALSE)
        results_tmp$f1 <<- f1
        start=end+1; end=start+k^2-1
        Q <- matrix(par[start:end],ncol=k, nrow=k, byrow=TRUE)
        results_tmp$Q <<- Q
        start=end+1; end=start+k-1
        f <- matrix(par[start:end],ncol=1, nrow=k, byrow=FALSE)
        results_tmp$f <<- f
        start=end+1; end=start+k-1
        b <- matrix(par[start:end],nrow=k, ncol=1, byrow=FALSE)
        results_tmp$b <<- b
        start=end+1; end=start
        mu0 <- par[start:end]
        results_tmp$mu0 <<- mu0
        start=end+1; end=start
        theta <- par[start:end]
        results_tmp$theta <<- theta
        start=end+1; end=start+k^2-1
        Q1 <- matrix(par[start:end],ncol=k, nrow=k, byrow=TRUE)
        results_tmp$Q1 <<- Q1
        # End reading parameters
    
        if(stopifbound) {
            for(i in 1:length(results_tmp)) {
                if(length(intersect(results_tmp[[i]],c(bounds$lower_bound[i], bounds$upper_bound[i]))) >= 1) {
                    cat("Parameter", names(results_tmp)[i], "achieved lower/upper bound. Process stopped.\n")
                    cat(results_tmp[[i]],"\n")
                    stopflag <- T
                    break
                }
            }
        }
    
        if(stopflag == FALSE) {
            dims <- dim(dat)
            res <- .Call("complikMD_quadratic_linear", dat, dims[1], dims[2], a, f1, Q, b, f, mu0, theta, Q1, k, pinv.tol, gomp)
            assign("results", results_tmp, envir=e)
            iteration <<- iteration + 1
            L.prev <<- res
      
            if(verbose) {
                cat("L = ", res,"\n")
                cat("Iteration: ", iteration,  "\nResults:\n") 
                print(results_tmp)
            }
        }
    
        return(as.numeric(-1*res))
    }

    # Optimization:
    if(verbose) {
        cat("Lower bound:\n")
        print(bounds$lower_bound)
        cat("Upper bound:\n")
        print(bounds$upper_bound)
    }
  
    ans <- nloptr(x0 = parameters, 
  		            eval_f = maxlik, 
  		            opts = opts,
                  lb = bounds$lower_bound, ub = bounds$upper_bound)
  
    final_results <- get("results",envir=e)
    #final_results <- results_tmp
    final_results[["status"]] <- ans$status
    final_results[["LogLik"]] <- L.prev
    final_results[["objective"]] <- ans$objective
    final_results[["message"]] <- ans$message
  
    # Check if any parameter achieved upper/lower limit and report it:
    limit <- FALSE
    for(i in 1:length(final_results)) {
        if(length(intersect(final_results[[i]],c(bounds$lower_bound[i], bounds$upper_bound[i]))) >= 1) {
            cat("Parameter", names(final_results)[i], "achieved lower/upper bound.\n")
            cat(final_results[[i]],"\n")
            limit <- TRUE
        }
    }
  
    final_results$limit <- limit
    #assign("results", final_results, envir=e)
    class(final_results) <- "spm.cont.quad.lin"
    invisible(final_results)
}