mmcif_sandwich {mmcif}

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()
}