joint_ms_opt: Optimizes the Lower Bound

In VAJointSurv: Variational Approximation for Joint Survival and Marker Models

Optimizes the Lower Bound

Description

Optimizes the Lower Bound

Usage

joint_ms_opt(
  object,
  par = object$start_val,
  rel_eps = 1e-08,
  max_it = 1000L,
  n_threads = object$max_threads,
  c1 = 1e-04,
  c2 = 0.9,
  use_bfgs = TRUE,
  trace = 0L,
  cg_tol = 0.5,
  strong_wolfe = TRUE,
  max_cg = 0L,
  pre_method = 3L,
  quad_rule = object$quad_rule,
  mask = integer(),
  cache_expansions = object$cache_expansions,
  gr_tol = -1,
  gh_quad_rule = object$gh_quad_rule
)

Arguments

object	a joint_ms object from joint_ms_ptr.
par	starting value.
rel_eps, max_it, c1, c2, use_bfgs, trace, cg_tol, strong_wolfe, max_cg, pre_method, mask, gr_tol	arguments to pass to the C++ version of psqn.
n_threads	number of threads to use. This is not supported on Windows.
quad_rule	list with nodes and weights for a quadrature rule for the integral from zero to one.
cache_expansions	TRUE if the expansions in the numerical integration in the survival parts of the lower bound should be cached (not recomputed). This requires more memory and may be an advantage particularly with expansions that take longer to compute (like ns_term and bs_term). The computation time may be worse particularly if you use more threads as the CPU cache is not well utilized.
gh_quad_rule	list with two numeric vectors called node and weight with Gauss–Hermite quadrature nodes and weights to handle delayed entry. A low number of quadrature nodes and weights is used when NULL is passed. This seems to work well when delayed entry happens at time with large marginal survival probabilities. The nodes and weights can be obtained e.g. from fastGHQuad::gaussHermiteData.

Value

A list with the following elements:

par	numeric vector of estimated model parameters.
value	numeric scalar with the value of optimized lower bound.
counts	integer vector with the function counts and the number of conjugate gradient iterations. See psqn.
convergence	logical for whether the optimization converged.

Examples

# load in the data
library(survival)
data(pbc, package = "survival")

# re-scale by year
pbcseq <- transform(pbcseq, day_use = day / 365.25)
pbc <- transform(pbc, time_use = time / 365.25)

# create the marker terms
m1 <- marker_term(
  log(bili) ~ 1, id = id, data = pbcseq,
  time_fixef = bs_term(day_use, df = 5L),
  time_rng = poly_term(day_use, degree = 1L, raw = TRUE, intercept = TRUE))
m2 <- marker_term(
  albumin ~ 1, id = id, data = pbcseq,
  time_fixef = bs_term(day_use, df = 5L),
  time_rng = poly_term(day_use, degree = 1L, raw = TRUE, intercept = TRUE))

# base knots on observed event times
bs_term_knots <-
  with(pbc, quantile(time_use[status == 2], probs = seq(0, 1, by = .2)))

boundary <- c(bs_term_knots[ c(1, length(bs_term_knots))])
interior <- c(bs_term_knots[-c(1, length(bs_term_knots))])

# create the survival term
s_term <- surv_term(
  Surv(time_use, status == 2) ~ 1, id = id, data = pbc,
  time_fixef = bs_term(time_use, Boundary.knots = boundary, knots = interior))

# create the C++ object to do the fitting
model_ptr <- joint_ms_ptr(
  markers = list(m1, m2), survival_terms = s_term,
  max_threads = 2L, ders = list(0L, c(0L, -1L)))

# find the starting values
start_vals <- joint_ms_start_val(model_ptr)

# optimize lower bound
fit <- joint_ms_opt(object = model_ptr, par = start_vals, gr_tol = .01)

# formatted maximum likelihood estimators
joint_ms_format(model_ptr, fit$par)

VAJointSurv documentation built on Aug. 14, 2022, 9:05 a.m.