Joint modeling of multivariate longitudinal outcomes and informative observation times with multivariate skew-t random effects.
SkewTIOT is an R package implementing a Monte-Carlo ECM (MCECM) estimator for
a joint model in which
- a multivariate longitudinal outcome (dimension
K) is observed at irregular, subject-specific visit times, and - those visit times are informative — they follow a recurrent-event Cox intensity whose frailty is the same subject-level random effect that drives the longitudinal mean.
The shared random effect is given a multivariate skew-t distribution, so
the model accommodates heavy tails and asymmetry in the latent subject effects.
The E-step uses a per-subject Gibbs sampler written in C++
(Rcpp + RcppArmadillo + OpenMP); the M-step updates are in R.
This repository is the companion code for the manuscript "Variance Estimation and Validation for a Joint Multivariate-Longitudinal / Informative-Observation- Time Model with Skew-t Random Effects."
The package contains compiled C++ and requires a working toolchain
(Xcode command-line tools on macOS, build-essential on Linux, Rtools on
Windows) plus the R packages listed under Imports in DESCRIPTION.
# install.packages("remotes")
remotes::install_github("dayusun/SkewTIOT")Or from a local clone:
R CMD INSTALL .A self-contained demo simulates a small data set (200 subjects, two correlated outcomes, informative visit times) and fits the model:
Rscript demo/run_demo.RIt prints the recovered beta, eta, gamma, and skewness Del next to the
true values (about a minute single-core). The script also shows, in comments,
how to switch to stiot_fnu_est() for profile-likelihood standard errors,
z-scores, AIC, and BIC.
| Function | Purpose |
|---|---|
stiot_fnu_est_only() |
MCECM point estimation at a fixed degrees-of-freedom nu (fastest). |
stiot_fnu_est() |
Point estimation plus profile-likelihood standard errors, z-scores, observed log-likelihood, AIC/BIC. |
stiot_fnu_cv_loss() |
Held-out marginal log-likelihood of a fitted model, for choosing nu by cross-validation. |
stiot_est_init() |
Warm-start initial values via sn::selm and reReg::reReg. |
stiot_tr_est() |
Variant that applies a user-supplied transformation to Y before the skew-t likelihood. |
All estimators expect data grouped contiguously by subject:
ID— length-all_nvector repeating the subject id per observation.Y—all_n × Kmatrix of longitudinal outcomes.X—all_n × p_Xcovariate matrix for the longitudinal mean. An intercept column is prepended internally whenifintercetp = TRUE(the default); passifintercetp = FALSEifXalready has one.Z— subject-level covariate matrix, one row per unique subject (nrow(Z) == length(unique(ID))), used for the visit-intensity Cox model.Time— length-all_nobservation times.censor— length-nper-subject end of follow-up (defaults to each subject's lastTimeifNULL).
A fitted object always carries beta, Sige, eta, gamma, lambda, Gam, Del, nu;
stiot_fnu_est() adds SE, zres, Covmat, logl_obs, AIC, BIC.
The per-subject E-step parallelizes over subjects. Pass a parallel PSOCK
cluster via the paracl argument to reuse it across calls, or leave
paracl = NULL (with ifpar = TRUE) to have the estimator create and tear down
its own cluster. OpenMP threading inside the C++ Gibbs sampler is controlled by
OMP_NUM_THREADS.
GPL (>= 2). See DESCRIPTION.