joint_ms_hess: Computes the Hessian

In VAJointSurv: Variational Approximation for Joint Survival and Marker Models

Computes the Hessian

Description

Computes the Hessian

Usage

joint_ms_hess(
  object,
  par,
  quad_rule = object$quad_rule,
  cache_expansions = object$cache_expansions,
  eps = 1e-04,
  scale = 2,
  tol = .Machine$double.eps^(3/5),
  order = 4L,
  gh_quad_rule = object$gh_quad_rule
)

Arguments

object	a joint_ms object from joint_ms_ptr.
par	parameter vector for where the lower bound is evaluated at.
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.
eps, scale, tol, order	parameter to pass to psqn. See psqn_hess.
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 two Hessian matrices:

hessian	Hessian matrix of the model parameters with the variational parameters profiled out.
hessian_all	Hessian matrix of the model and variational parameters.

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)

# compute the Hessian
hess <- joint_ms_hess(object = model_ptr,par = fit$par)

# standard errors of the parameters
library(Matrix)
sqrt(diag(solve(hess$hessian))) 

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