AR1_FRAILTY: Fit a Functional Regression with AutoregressIve fraiLTY (FRAILTY) model for Recurrent Event Data

Description

Jointly model longitudinal measurements and recurrent events, accommodating both scalar and functional predictors while capturing time-dependent correlations among events. The FRAILTY method employs a two-step estimation procedure. First, functional principal component analysis through conditional expectation (PACE) is applied to extract key temporal features from sparse, irregular longitudinal data. Second, the obtained scores are incorporated into a dynamic recurrent frailty model with an autoregressive structure to account for within-subject correlations across recurrent events. This function works only for univariate functional data.

Usage

AR1_FRAILTY(
  formula,
  sdat,
  fdat,
  para0 = c(0.5, 0.5),
  nbasis = 10,
  pve = 0.9,
  npc = NULL,
  makePD = FALSE,
  cov.weight.type = c("none", "counts"),
  iter.max = 50,
  eps = 1e-06
)

Value

A funsurv object containing the following components:

beta: Estimation of coefficients of scalar covariates and FPC scores. Including estimated values, standard errors, and p-values
beta_vcov: Estimated variance-covariance of the estimates of beta
eAR: Estimation of variance components ($\theta^2$ and $\rho$)
eAR_vcov: Estimated variance of estimates of $\theta^2$ and $\rho$
frailties: Estimated frailty terms (random effects)
basesurv: Estimated baseline survival probability
time: Time points associated with baseline survival probability
FPC: Functional principal components

Arguments

formula: A formula, with the response on the left of a ~ operator being a Recur object as returned by function Recur in reda, and scalar covariates on the right.
sdat: A data frame containing subject IDs, time-to-event outcomes (starting time, end point, censoring time and event status), and scalar covariates
fdat: A data frame containing subject IDs, longitudinal measurements, and the corresponding time points for each measurement.
para0: A vector of initial values for $\theta^2$ and auto-regressive coefficient $\rho$. Both default to 0.5.
nbasis: An integer, representing the number of B-spline basis functions used for estimation of the mean function and bivariate smoothing of the covariance surface. Defaults to 10 (cf. fpca.sc in refund).
pve: A numeric value between 0 and 1, the proportion of variance explained: used to choose the number of principal components. Defaults to 0.9 (cf. fpca.sc in refund).
npc: An integer, giving a prespecified value for the number of principal components. Defaults to NULL. If given, this overrides pve (cf. fpca.sc in refund).
makePD: Logical: should positive definiteness be enforced for the covariance surface estimate? Defaults to FALSE (cf. fpca.sc in refund).
cov.weight.type: The type of weighting used for the smooth covariance estimate. Defaults to "none", i.e. no weighting. Alternatively, "counts" (corresponds to fpca.sc in refund) weights the pointwise estimates of the covariance function by the number of observation points.
iter.max: Maximum number of iterations for both inner iteration and outer iteration. Defaults to 50.
eps: Tolerance criteria for a possible infinite coefficient value. Defaults to 1e-6.

Details

Model specification:

Let $T_{ij}$ denote the time of the $j$th event for subject $i$, and let $C_i$ represent the censoring time. The observed event time, accounting for right censoring, is $\widetilde{T}_{ij}=\min(T_{ij}, C_i)$, and $\delta_{ij}=I(T_{ij}\leq C_i)$ serves as an indicator of whether the $j$th event for subject $i$ is observed. The hazard function is specified as $$h(t; \boldsymbol{Z}_i, {X}_i(\cdot))=h_{0}(t-t_{i,j-1}) \exp \left(\eta_{ij}\right),$$ where $h_0(\cdot)$ is the baseline hazard function, and $\eta_{ij} = \bm{\alpha}^{\top}\boldsymbol{Z}_i +\int_{t_{i, j-1}}^{t}{X}_{i}(s)\beta(s)ds + v_{ij}$. Here, $t_{i, j-1}$ is the previous event time with $t_{i0} = 0$. $\bm{\alpha}$ is the fixed effect parameter associated with the time-invariant covariates $\boldsymbol{Z}_i$, and $\beta(t)$ is a time-varying coefficient that captures the effect of functional predictor $X_{i}(t)$ on the hazard rate of recurrent events.

Examples

Run this code

data(simDat)

fit <- AR1_FRAILTY(Recur(t_start %to% t_stop, id, status) ~ z1,
                   sdat = sdat, fdat = fdat, iter.max = 30)

summary(fit)