mmcif_sandwich {mmcif} | R Documentation |
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.
mmcif_sandwich(
object,
par,
ghq_data = object$ghq_data,
n_threads = 1L,
eps = 0.01,
scale = 2,
tol = 1e-08,
order = 3L
)
object |
an object from |
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 |
n_threads |
the number of threads to use. |
eps |
determines the step size in the numerical differentiation using
|
scale |
scaling factor in the Richardson extrapolation. Each step is
smaller by a factor |
tol |
relative convergence criteria in the extrapolation given
by |
order |
maximum number of iteration of the Richardson extrapolation. |
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.
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.
mmcif_fit
and mmcif_data
.
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()
}