R/spm_time-dependent.R

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

#' A function for the model with time-dependent model parameters.
#'
#'@param x Input data table.
#'@param start A list of starting parameters, default: 
#'\code{start=list(a=-0.5, f1=80, Q=2e-8, f=80, b=5, mu0=1e-5)}.
#'@param frm A list of formulas that define age (time) - dependency. 
#'Default: \code{frm=list(at="a", f1t="f1", Qt="Q", ft="f", bt="b", 
#'mu0t="mu0")}.
#'@param stopifbound Estimation stops if at least one parameter 
#'achieves lower or upper boundaries. Default: \code{FALSE}.
#'@param lb Lower bound of parameters under estimation.
#'@param ub Upper bound of parameters under estimation.
#'@param verbose Turns on verbosing output.
#'@param opts A list of options for \code{nloptr}.
#'Default value: \code{opt=list(algorithm="NLOPT_LN_NELDERMEAD", 
#'maxeval=100, ftol_rel=1e-8)}.
#'@param lrtest Indicates should Likelihood-Ratio test be performed.
#'Possible values: \code{TRUE}, \code{H01}, \code{H02}, \code{H03},
#'\code{H04}, \code{H05} (see package Vignette for details)
#'Default value: \code{FALSE}.
#'Please see \code{nloptr} documentation for more information.
#'@return A set of estimates of \code{a}, \code{f1}, \code{Q}, 
#'\code{f}, \code{b}, \code{mu0}.
#'@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 \code{nloptr} optimization function 
#'(see documentation for nloptr package).
#'@references Yashin, A. et al (2007), Health decline, 
#'aging and mortality: how are they related? 
#'Biogerontology, 8(3), 291-302.<DOI:10.1007/s10522-006-9073-3>.
#'@export
#'@examples
#'library(stpm)
#'set.seed(123)
#'#Data preparation:
#'n <- 5
#'data <- simdata_time_dep(N=n)
#'# Estimation:
#'opt.par <- spm_time_dep(data)
#'opt.par
spm_time_dep <- function(x, 
                         start=list(a=-0.05, f1=80, Q=2e-8, f=80, b=5, mu0=1e-3),
                         frm=list(at="a", f1t="f1", Qt="Q", ft="f", 
                                  bt="b", mu0t="mu0"), 
                         stopifbound=FALSE, 
                         lb=NULL, ub=NULL,
                         verbose=FALSE, 
                         opts=NULL,
                         lrtest=FALSE) {
    
    if(is.null(opts) | !is.list(opts)) {
        opts=list(algorithm="NLOPT_LN_NELDERMEAD")
    } else if(is.null(opts$algorithm)) {
        opts$algorithm="NLOPT_LN_NELDERMEAD"
    }
    
    hypotheses <- c("H01", "H02", "H03", "H04", "H05")
    
    if(is.null(start)) {
        start <- list(a=-0.05, f1=80, Q=2e-8, f=80, b=5, mu0=1e-3)
    }
  
    res <- NA
    if(lrtest == FALSE) {
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                     verbose = verbose, opts=opts)
        res.null <- NA
    } else if(lrtest=="H01" | lrtest==TRUE) { # a_q = 0, b_q = 0
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                   verbose = verbose, opts=opts)
        frm[["Qt"]] <- "0"
        start[["Q"]] <- NULL
        lb <- lb[-which(names(lb) == "Q")]
        ub <- ub[-which(names(ub) == "Q")]
        
        res.null <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                     verbose = verbose, opts=opts)
        
        lr.test.pval <- LRTest(res$LogLik, res.null$LogLik)
        res[["lr.test.pval"]] <- lr.test.pval
    } else if(lrtest=="H02") { # a_q = const, b_q = 0
        const.q <- frm[["Qt"]]
        frm[["Qt"]] <- "Q"
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                   verbose = verbose, opts=opts)
        
        frm[["Qt"]] <- const.q
        start[["Q"]] <- NULL
        lb <- lb[-which(names(lb) == "Q")]
        ub <- ub[-which(names(ub) == "Q")]
      
        res.null <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                        verbose = verbose, opts=opts)
        res.null[["Q"]] <- as.numeric(const.q)
        lr.test.pval <- LRTest(res$LogLik, res.null$LogLik)
        res[["lr.test.pval"]] <- lr.test.pval
    } else if(lrtest=="H03") { # a_f1 = 0, b_f1 = 0
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                   verbose = verbose, opts=opts)
        frm[["f1t"]] <- "0"
        start[["f1"]] <- NULL
        lb <- lb[-which(names(lb) == "f1")]
        ub <- ub[-which(names(ub) == "f1")]
      
        res.null <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                        verbose = verbose, opts=opts)
      
        lr.test.pval <- LRTest(res$LogLik, res.null$LogLik)
        res[["lr.test.pval"]] <- lr.test.pval
    } else if(lrtest=="H04") { # a_f1 = const, b_f1 = 0
        const.q <- frm[["f1t"]]
        frm[["f1t"]] <- "f1"
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                   verbose = verbose, opts=opts)
        frm[["f1t"]] <- const.q
        start[["f1"]] <- NULL
        lb <- lb[-which(names(lb) == "Q")]
        ub <- ub[-which(names(ub) == "Q")]
      
        res.null <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                        verbose = verbose, opts=opts)
      
        lr.test.pval <- LRTest(res$LogLik, res.null$LogLik)
        res[["lr.test.pval"]] <- lr.test.pval
    } else if(lrtest=="H05") { # a_y = const, b_y = 0
        const.q <- frm[["at"]]
        frm[["at"]] <- "a"
        res <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                   verbose = verbose, opts=opts)
      
        frm[["at"]] <- const.q
        start[["a"]] <- NULL
        lb <- lb[-which(names(lb) == "a")]
        ub <- ub[-which(names(ub) == "a")]
      
        res.null <- spm_time_dep_internal(x=x, start=start, frm=frm, stopifbound=stopifbound, lb=lb, ub=ub,
                                        verbose = verbose, opts=opts)
        res.null[["a"]] <- as.numeric(const.q)
        lr.test.pval <- LRTest(res$LogLik, res.null$LogLik)
        res[["lr.test.pval"]] <- lr.test.pval
    }
    
    if((lrtest == TRUE) | (lrtest %in% hypotheses)) {
        return(list(alternative=res, null=res.null))
    } else {
        return(res)
    }
}


spm_time_dep_internal <- function(x, start, frm, stopifbound, lb, ub, verbose, opts) {
    
    #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_GN_CRS2_LM", 
    #                      "NLOPT_LN_COBYLA", "NLOPT_LN_NEWUOA",
    #                      "NLOPT_LN_NEWUOA_BOUND", "NLOPT_LN_NELDERMEAD",
    #                      "NLOPT_LN_SBPLX",
    #                      "NLOPT_LN_BOBYQA", "NLOPT_GN_ISRES")
    
    e <- new.env()
    
    #--------------Begin of optimize function-------------------#
    optimize <- function(data, starting_params,  formulas, verbose, 
                         lb, ub, 
                         stopifbound,
                         opts) {
    
        final_res <- list()
    
        # Current results:
        results <- list()
    
        N <- dim(data)[1]
        at <- NULL
        f1t <- NULL
        Qt <- NULL
        ft <- NULL
        bt <- NULL
        mu0t <- NULL
        variables <- c()
        
        # Assigning parameters:
        p.const.ind <- c()
        p.coeff.ind <- c()
        
        model.coef <- c("at", "f1t", "Qt", "ft", "bt", "mu0t")
        
        for(c in model.coef) {
            parsed <- parse_parameters(formulas, c, p.const.ind, p.coeff.ind, variables)
            p.const.ind <- parsed$p.const.ind
            p.coeff.ind <- parsed$p.coeff.ind
            variables <- parsed$variables
        }
        
        ########=============
        p.const.ind <- unique(p.const.ind)
        p.coeff.ind <- unique(p.coeff.ind)
        variables <- suppressWarnings(variables[which(is.na(as.numeric(variables)))])
        
        
        for(name in names(starting_params)) {
            if(!(name %in% names(lb))) {
                lb[[name]] <- ifelse(starting_params[[name]] < 0, 
                                 starting_params[[name]] + 0.25*starting_params[[name]], 
                                 starting_params[[name]] - 0.25*starting_params[[name]])
            }
          
            if(!(name %in% names(ub))) {
                ub[[name]] <- ifelse(starting_params[[name]] < 0, 
                                 starting_params[[name]] - 0.25*starting_params[[name]],
                                 starting_params[[name]] + 0.25*starting_params[[name]])
            }
        }
        
        
        
        # Starting parameters, lower and upper boundaries:
        stpar <- c()
        lower_bound <- c()
        upper_bound <- c()
        for(name in names(starting_params)) {
            stpar <- c(stpar, starting_params[[name]])
            lower_bound <- c(lower_bound, lb[[name]])
            upper_bound <- c(upper_bound, ub[[name]])
        }
        names(stpar) <- variables
    
        for(p in names(stpar)) {
            results[[p]] <- stpar[p]
            #assign(x=p, value=stpar[p], envir=.GlobalEnv)
            assign_to_global(what=p,value=stpar[p])
        }
    
    
        comp_func_params <- function(astring, f1string, qstring, fstring, bstring, mu0string) {
            #at <<- eval(parse(text = paste('at <<- function(t) { return(' , astring , ')}', sep='')))
            #f1t <<- eval(parse(text = paste('f1t <<- function(t) { return(' , f1string , ')}', sep='')))
            #Qt <<- eval(parse(text = paste('Qt <<- function(t) { return(' , qstring , ')}', sep='')))
            #ft <<- eval(parse(text = paste('ft <<- function(t) { return(' , fstring , ')}', sep='')))
            #bt <<- eval(parse(text = paste('bt <<- function(t) { return(' , fstring , ')}', sep='')))
            #mu0t <<- eval(parse(text = paste('mu0t <<- function(t) { return(' , mu0string , ')}', sep='')))
            
            at <<- eval(bquote(function(t) .(parse(text = astring)[[1]])))
            f1t <<- eval(bquote(function(t) .(parse(text = f1string)[[1]]))) 
            Qt <<- eval(bquote(function(t) .(parse(text = qstring)[[1]])))
            ft <<- eval(bquote(function(t) .(parse(text = fstring)[[1]])))
            bt <<- eval(bquote(function(t) .(parse(text = bstring)[[1]])))
            mu0t <<- eval(bquote(function(t) .(parse(text = mu0string)[[1]])))
        }
    
        iteration <- 0
        L.prev <- 0
    
        maxlik_t <- function(params) {
            stopflag <- FALSE
      
            if(verbose)
                cat("Iteration: ", iteration, "\n")
      
            names(params) <- names(stpar)
            
      
            for(p in names(stpar)) {
                #assign(p, params[[p]], envir=e)
                #assign(p, params[[p]], envir=.GlobalEnv)
                assign_to_global(what=p,value=params[[p]])
                results[[p]] <<- params[[p]]
                if(verbose)
                    cat(paste(p, results[[p]]), " ")
            }
            
            
            sigma_sq <- function(t1, t2) {
                # t2 = t_{j}, t1 = t_{j-1}
                ans <- bt(t1)*(t2-t1)
                ans
            }
      
            m <- function(y, t1, t2) {
                # y = y_{j-1}, t1 = t_{j-1}, t2 = t_{j}
                ans <- y + at(t1)*(y - f1t(t1))*(t2 - t1)
                ans
            }
      
            mu <- function(y,t) {
                ans <- mu0t(t) + (y - ft(t))^2*Qt(t)
                ans
            }
      
            if(stopifbound) {
                for(i in 1:length(results)) {
                    if(length(intersect(results[[i]],c(lower_bound[i], upper_bound[i]))) >= 2) {
                        cat("Parameter", names(results)[i], "achieved lower/upper bound. Process stopped.\n")
                        cat(results[[i]],"\n")
                        stopflag <- TRUE
                        break
                    }
                }
            }
      
            L <- 0
      
            if(stopflag == FALSE) {
        
                for(i in 1:N) {
                    delta <- data[i,1]
                    t1 <- data[i, 2]; t2 <- data[i, 3]
                    if(t1 < t2) {
                        ind <- ifelse(is.na(data[i, 5]), 0, 1)
                        S <- exp(-1*mu(data[i, 4],t1)*(t2-t1))
                        if(S <= 1e-5) {
                            S <- 1e-5
                        }
                        if(ind == 0) {
                            L <- L + (1 - delta)*log(S) + delta*log(1-S)
                            #if(is.nan(L)) {
                            #  print(mu(data[i, 4],t1)*(t2-t1))
                            #  print(S)
                            #  print(data[i,])
                            #  stop()
                            #}
                        } else {
                            yj <- data[i,5]
                            mj <- m(data[i,4], t1, t2)
                            pn <- -1*log(sqrt(2*pi)*sqrt(sigma_sq(t1, t2))) - (yj - mj)^2/(2*sigma_sq(t1, t2))
                            L <- L + pn + (1 - delta)*log(S) + delta*log(1-S)
                        }
            
                    }
                }
        
                #assign("results", results, envir=e)
                #assign("results", results, envir=.GlobalEnv)
                assign_to_global(what="results",value=results)
                
                iteration <<- iteration + 1
                L.prev <<- L
        
        
                if(verbose) {
                    cat("\n")
                    cat(paste("L", L.prev), "\n")
                }
        
            } else {
        
                cat("Optimization stopped. Parametes achieved lower or upper bound.\nPerhaps you need more data or these returned parameters might be enough.\n")
                print("###########################################################")
                L <- NA
        
            }
      
            L <- -1*L
            L
        }
    
        comp_func_params(formulas$at, formulas$f1t, formulas$Qt, formulas$ft, formulas$bt, formulas$mu0t)
    
        #if( setequal(names(starting_params), variables) == FALSE) {
        #  stop("Provided set of function parameters is not equal to that one provided in starting list or vise-versa.")
        #}
    
        if(verbose == TRUE) {
            cat("Variables:\n")
            cat(variables,"\n")
            cat("Functions:\n")
            print(at)
            print(f1t)
            print(Qt)
            print(ft)
            print(bt)
            print(mu0t)
        }
    
        # Optimization:
        if(verbose) {
            cat("Lower bound:\n")
            print(lower_bound)
            cat("Upper bound:\n")
            print(upper_bound)
      
        }
    
    
        tryCatch({
            ans <- nloptr(x0 = unlist(stpar), 
                            eval_f = maxlik_t, opts = opts,
                            lb = lower_bound, ub = upper_bound)
            i <- 1
            
            for(p in names(stpar)) {
                results[[p]] <<- ans$solution[i]
                i <- i + 1
                if(verbose)
                    cat(paste(p, get(p)), " ")
            }
            results[["status"]] <- ans$status
            results[["LogLik"]] <- L.prev
            results[["objective"]] <- ans$objective
            results[["message"]] <- ans$message
            
        },  error=function(e) {
            if(verbose  == TRUE) {print(e)}
        }, finally=NA)
    
        final_res <- results
        return(final_res)
    }
    #---------------------End of optimize---------------------------#
  
    formulas <- frm
    data <- as.matrix(x)
    data <- data[, 2:dim(data)[2]]
    formulas.work = list(at="a", f1t="f1", Qt="Q", ft="f", bt="b", mu0t="mu0")
    for(item in names(formulas)) {
        formulas.work[[item]] <- formulas[[item]]
    }
    # Optimization:
    res <- optimize(data, start, formulas.work, verbose, lb, ub, stopifbound, opts)
    
    invisible(res)
}


parse_parameters <- function(formulas, parameter, p.const.ind, p.coeff.ind, variables) {
  #---
  parameters <- trim(unlist(strsplit(formulas[[parameter]],"[\\+\\*\\(\\)]",fixed=F)))
  parameters <- parameters[which(!(parameters %in% c("t","exp")))]
  ee <- new.env()
  #for(p in parameters) {assign(p,NULL, envir = .GlobalEnv); variables <- c(variables, p);}
  for(p in parameters) {assign(p,NULL, envir=ee); variables <- c(variables, p);}
  variables <- unique(variables)
  p.constants <- c()
  p.coeffs <- c()
  p.temp <- trim(unlist(strsplit(formulas$at,"[\\+\\-\\(\\)]",fixed=F)))
  for(i in 1:length(p.temp)) {
    if(grepl('t', p.temp[i]) == FALSE) {
      p.constants <- c(p.constants, p.temp[i])
    } else {
      p.temp.coeff <- trim(unlist(strsplit(p.temp[i],"[\\*]",fixed=F)))
      p.coeffs <- c(p.coeffs, p.temp.coeff[which(!(p.temp.coeff %in% c("t","exp")))])
    }
  }
  
  for(i in 1:length(p.constants)) {
    p.const.ind <- c(p.const.ind, which(variables == p.constants[i]))
  }
  for(i in 1:length(p.coeffs)) {
    p.coeff.ind <- c(p.coeff.ind, which(variables == p.coeffs[i]))
  }
  
  
  return(list(p.const.ind=p.const.ind, p.coeff.ind=p.coeff.ind, variables=variables))
  
}