joint_ms_profile {VAJointSurv}

Approximate Likelihood Ratio based Confidence Intervals

Description

Approximate Likelihood Ratio based Confidence Intervals

Usage

joint_ms_profile(
  object,
  opt_out,
  which_prof,
  delta,
  level = 0.95,
  max_step = 15L,
  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,
  verbose = TRUE,
  mask = integer(),
  cache_expansions = object$cache_expansions,
  gr_tol = -1,
  hess = NULL
)

Arguments

object

a joint_ms object from joint_ms_ptr.

opt_out

maximum lower bound estimator from joint_ms_opt.

which_prof

index of the parameter to profile.

delta

numeric scalar greater than zero for the initial step size. Steps are made of size 2^(iteration - 1) * delta. A guess of the standard deviation is a good value.

level

confidence level.

max_step

maximum number of steps to take in each direction when constructing the approximate profile likelihood curve.

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.

verbose

logical for whether to print output during the construction of the approximate profile likelihood curve.

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.

hess

the Hessian from joint_ms_hess. It is used to get better starting values along the profile likelihood curve. Use NULL if it is not passed.

Value

A list with the following elements:

confs

profile likelihood based confidence interval.

xs

the value of the parameter at which the profile likelihood is evaluated at.

p_log_Lik

numeric scalar with the profile log-likelihood.

data

list of lists of the output of each point where the profile likelihood is evaluated with the optimal parameter values of the other parameters given the constrained value of the parameter that is being profiled and the optimal value of the lower bound.

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)

# compute the standard errors
library(Matrix)
se <- sqrt(diag(solve(hess$hessian)))

# find index for the first association parameter
which_prof <- model_ptr$indices$survival[[1]]$associations[1]

# initial step size for finding the confidence interval limits
delta <- 2*se[which_prof]

# compute profile likelihood based confidence interval
# for the first association parameter
profile_CI <- joint_ms_profile(
  object = model_ptr, opt_out = fit, which_prof = which_prof,
  delta= delta, gr_tol = .01)

# comparison of CIs
profile_CI$confs
fit$par[which_prof]+c(-1,1)*qnorm(0.975)*se[which_prof]