mmcif_sandwich: Computes the Sandwich Estimator

In mmcif: Mixed Multivariate Cumulative Incidence Functions

Computes the Sandwich Estimator

Description

Computes the sandwich estimator of the covariance matrix. The parameter that is passed is using the log Cholesky decomposition. The Hessian is computed using numerical differentiation with Richardson extrapolation to refine the estimate.

Usage

mmcif_sandwich(
  object,
  par,
  ghq_data = object$ghq_data,
  n_threads = 1L,
  eps = 0.01,
  scale = 2,
  tol = 1e-08,
  order = 3L
)

Arguments

object	an object from mmcif_data.
par	numeric vector with the parameters to compute the sandwich estimator at.
ghq_data	the Gauss-Hermite quadrature nodes and weights to use. It should be a list with two elements called "node" and "weight". A default is provided if NULL is passed.
n_threads	the number of threads to use.
eps	determines the step size in the numerical differentiation using max(sqrt(.Machine$double.eps), |par[i]| * eps) for each parameter i.
scale	scaling factor in the Richardson extrapolation. Each step is smaller by a factor scale.
tol	relative convergence criteria in the extrapolation given by max(tol, |g[i]| * tol) with g being the gradient and for each parameter i.
order	maximum number of iteration of the Richardson extrapolation.

Value

The sandwich estimator along attributes called

• "meat" for the "meat" of the sandwich estimator.

• "hessian" for the Hessian of the log composite likelihood.

• "res vcov" which is the sandwich estimator where the last elements are the upper triangle of the covariance matrix of the random effects rather than the log Cholesky decomposition of the matrix.

References

Cederkvist, L., Holst, K. K., Andersen, K. K., & Scheike, T. H. (2019). Modeling the cumulative incidence function of multivariate competing risks data allowing for within-cluster dependence of risk and timing. Biostatistics, Apr 1, 20(2), 199-217.

See Also

mmcif_fit and mmcif_data.

Examples

if(require(mets)){
  # prepare the data
  data(prt)

  # truncate the time
  max_time <- 90
  prt <- within(prt, {
    status[time >= max_time] <- 0
    time <- pmin(time, max_time)
  })

  # select the DZ twins and re-code the status
  prt_use <- subset(prt, zyg == "DZ") |>
    transform(status = ifelse(status == 0, 3L, status))

  # randomly sub-sample
  set.seed(1)
  prt_use <- subset(
    prt_use, id %in% sample(unique(id), length(unique(id)) %/% 10L))

  n_threads <- 2L
  mmcif_obj <- mmcif_data(
    ~ country - 1, prt_use, status, time, id, max_time,
    2L, strata = country)

  # get the staring values
  start_vals <- mmcif_start_values(mmcif_obj, n_threads = n_threads)

  # estimate the parameters
  ests <- mmcif_fit(start_vals$upper, mmcif_obj, n_threads = n_threads)

  # get the sandwich estimator
  vcov_est <- mmcif_sandwich(
    mmcif_obj, ests$par, n_threads = n_threads, order = 2L)

  # show the parameter estimates along with the standard errors
  rbind(Estimate = ests$par,
        SE = sqrt(diag(vcov_est))) |>
    print()

  # show the upper triangle of the covariance matrix and the SEs
  rbind(`Estimate (vcov)` = tail(ests$par, 10) |> log_chol_inv() |>
          (\(x) x[upper.tri(x, TRUE)])() ,
        SE = attr(vcov_est, "res vcov") |> diag() |> sqrt() |> tail(10)) |>
    print()
}

mmcif documentation built on July 18, 2022, 1:06 a.m.