Modeling semi-competing risks data as a longitudinal bivariate process

The proposed framework builds on a novel representation of observed semi-competing risks outcome data. Prior to describing this, it is useful to consider how such data are typically represented. Towards this, let $L$ and $C$ denote the study entry and right censoring times, respectively. The observed outcome data for the $i^{th}$ study participant is $\mathcal{D}_{i}=(L_{i},\widetilde{T}_{1i},\delta_{1i},\widetilde{T}_{2i},\delta_{2i})$, where $\widetilde{T}_{1i}=\min(T_{1i},T_{2i},C_{i})$, $\widetilde{T}_{2i}=\min(T_{2i},C_{i})$, $\delta_{1i}=I\{\widetilde{T}_{1i}=T_{1i}\}$, and $\delta_{2i}=I\{\widetilde{T}_{2i}=T_{2i}\}$. Note, $\delta_{1i}$ and $\delta_{2i}$ indicate whether the study unit is observed to experience the non-terminal event and terminal event, respectively.

Let $\mbox{$\tau$}=\{\tau_{0},\ldots,\tau_{K}\}$ be a set of user-specified points, with $\tau_{k-1}<\tau_{k}$ for $k=1,\ldots,K$, that define a partition of the analysis time scale. Towards the study of AD using data from the ACT study, for example, one could choose age beyond 65 years as the time scale and $\mbox{$\tau$}=\{65,70,75,...,100\}$. Let $Y_{1i,k}=I\{T_{1i}\leq\tau_{k}\}$ and $Y_{2i,k}=I\{T_{2i}\leq\tau_{k}\}$ be indicators of whether participant $i$ experienced the non-terminal and terminal events by time $\tau_{k}$, respectively. Furthermore, let $k^{l}_{i}\in\{1,\ldots,K\}$ be the index of the first time point in $\tau$ after $L_{i}$ and $k^{r}_{i}\in\{1,\ldots,K\}$ the index of the first time point in $\tau$ after $\widetilde{T}_{2i}$. The observed outcome data on the partition $\tau$ can then be represented as a longitudinal bivariate process, specifically $\mbox{$Y$}_{\hskip-1.8063pti}=\{(Y_{1i,k},Y_{2i,k});k=k^{l}_{i},\ldots,k^{r}_{i}\}$. Figure A.1 in the Supplementary Materials provides a graphical representation of the observed outcome data for four hypothetical study participants.

Finally, let $\mbox{$X$}_{i}(\tau_{k})$ denote a vector of possibly time-dependent covariates measured at time $\tau_{k-1}$. Given this, let $\overline{\mbox{$X$}}_{i,k}=(\mbox{$X$}_{i}(\tau_{1}),\ldots,\mbox{$X$}_{i}(\tau_{k}))$ denote the history of all covariate information (i.e. both time-invariant and time-varying) up to time point $\tau_{k-1}$.

Let $\overline{\mbox{$X$}}_{i}\equiv\overline{\mbox{$X$}}_{i,{k^{r}_{i}}}$ is the totality of all covariate data observed for the $i^{th}$ participant. At the outset we assume that the joint distribution of the observed outcome data for the $i^{th}$ study participant, $\mbox{P}(\mbox{$Y$}_{\hskip-1.8063pti}=\mbox{$y$}_{i}|\ \overline{\mbox{$X$}}_{i})$, can be decomposed as the product:

$$\prod_{k=k^{l}_{i}}^{k^{r}_{i}}\mbox{P}(Y_{1i,k}=y_{1i,k},Y_{2i,k}=y_{2i,k}|\ Y_{1i,k-1}=y_{1i,k-1},Y_{2i,k-1}=y_{2i,k-1},\overline{\mbox{$X$}}_{i,k}),$$

(1)

where we define $Y_{1i,0}=Y_{2i,0}=0$. Underpinning this decomposition is a Markov-type assumption that the joint probability of the two events in a given interval depends on the history of the two events solely through what is known at the start of the interval, conditional on the totality of the covariate information to that point.

As shown in Section A.4 of the Supplementary Materials, the components of expression (1) can be written in terms of:

	$\displaystyle\pi_{1i,k}$	$\displaystyle=$	$\displaystyle\mbox{P}(Y_{1i,k}=1|\ Y_{1i,k-1}=0,Y_{2i,k-1}=0,\overline{\mbox{$X$}}_{i,k})$		(2)
	$\displaystyle\pi_{2i,k}(y_{1})$	$\displaystyle=$	$\displaystyle\mbox{P}(Y_{2i,k}=1|\ Y_{1i,k-1}=y_{1},Y_{2i,k-1}=0,\overline{\mbox{$X$}}_{i,k})$		(3)
	$\displaystyle\theta_{i,k}$	$\displaystyle=$	$\displaystyle\mbox{OR}(Y_{1i,k},Y_{2i,k}|\ Y_{1i,k-1}=0,Y_{2i,k-1}=0,\overline{\mbox{$X$}}_{i,k})$		(4)

for $y_{1}=0,1$ and where

$$\mbox{OR}(Y_{1i,k},Y_{2i,k}|\ \ldots)\ =\ \frac{\mbox{P}(Y_{1i,k}=1,Y_{2i,k}=1|\ \ldots)\times\mbox{P}(Y_{1i,k}=0,Y_{2i,k}=0|\ \ldots)}{\mbox{P}(Y_{1i,k}=0,Y_{2i,k}=1|\ \ldots)\times\mbox{P}(Y_{1i,k}=1,Y_{2i,k}=0|\ \ldots)}.$$

We emphasize that the interpretations of these three quantities are specific to the partition $\tau$. Specifically, $\pi_{1i,k}$ is the cumulative probability of experiencing the non-terminal event by time $\tau_{k}$, given that neither event had occurred by time $\tau_{k-1}$. Similarly, $\pi_{2i,k}(\cdot)$ is the cumulative probability of experiencing the terminal event by time $\tau_{k}$, given that the individual is alive at time $\tau_{k-1}$. Note, both of these correspond to quantities that are modeled in discrete time survival analyses [24, 31]. Finally, $\theta_{i,k}$ is the cross-sectional odds ratio for the 2$\times$2 table corresponding to the four possible observed ($Y_{1i,k},Y_{2i,k}$) outcome vectors at time $\tau_{k}$, given that neither the non-terminal nor the terminal event had occurred by time $\tau_{k-1}$. Given its clear importance, we return to the choice of $\tau$ in Section 3.5.

We proceed with modeling by placing structure on the components given by expressions (2)-(4) as a function of covariates. Specifically, we propose that $\pi_{1i,k}$, $\pi_{2i,k}(\cdot)$, and $\theta_{i,k}$ be modeled via the following regressions:

	$\displaystyle\pi_{1i,k}$	$\displaystyle=$	$\displaystyle g_{1}^{-1}\{\alpha_{1,k}\ +\ f_{1}(\overline{\mbox{$X$}}_{i1,k};\mbox{$\beta$}_{1})\}$		(5)
	$\displaystyle\pi_{2i,k}(y_{1})$	$\displaystyle=$	$\displaystyle g_{2}^{-1}\{\alpha_{2,k}\ +\ f_{2}(\overline{\mbox{$X$}}_{i2,k},Y_{1i,k-1};\mbox{$\beta$}_{2})\}\hskip 21.68121pt$		(6)
	$\displaystyle\theta_{i,k}$	$\displaystyle=$	$\displaystyle g_{\theta}^{-1}\{\alpha_{\theta,k}\ +\ f_{\theta}(\overline{\mbox{$X$}}_{i\theta,k};\mbox{$\beta$}_{\theta})\},$		(7)

where $g_{1}(\cdot)$, $g_{2}(\cdot)$ and $g_{\theta}(\cdot)$ are user-specified link functions (e.g. the logistic or log link), where $\overline{\mbox{$X$}}_{i1,k}$, $\overline{\mbox{$X$}}_{i2,k}$, $\overline{\mbox{$X$}}_{i\theta,k}$ are each subsets of $\overline{\mbox{$X$}}_{i,k}$, and where $f_{1}(\cdot)$, $f_{2}(\cdot)$ and $f_{\theta}(\cdot)$ are user-specified functions that characterize how the respective quantities depend on the covariates. For example, one may adopt a linear predictor specification with no interactions between the components of $\overline{\mbox{$X$}}_{i2,k}$ and $Y_{1i,k-1}$ in the model for $\pi_{2i,k}(y_{1})$ by setting $f_{2}(\overline{\mbox{$X$}}_{i2,k},Y_{1i,k-1};\mbox{$\beta$}_{2})\equiv\overline{\mbox{$X$}}_{i2,k}^{T}\mbox{$\beta$}_{2,X}+Y_{1i,k-1}\beta_{2,y}$.

While their precise interpretations will depend on the partition of the time axis and the chosen link functions, $\mbox{$\beta$}=(\mbox{$\beta$}_{1},\mbox{$\beta$}_{2},\mbox{$\beta$}_{\theta})$ characterize the impact of covariates; Section 3.4 considers the role and interpretation of components of $\mbox{$\beta$}_{2}$ that correspond to $Y_{1i,k-1}$. Finally, $\mbox{$\alpha$}_{1}$ = $(\alpha_{1,1},\ldots,\alpha_{1,K})$, $\mbox{$\alpha$}_{2}$ = $(\alpha_{2,1},\ldots,\alpha_{2,K})$, and $\mbox{$\alpha$}_{\theta}$ = $(\alpha_{\theta,1},\ldots,\alpha_{\theta,K})$ represent baseline time trends in $\pi_{1i,k}$, $\pi_{2i,k}(\cdot)$, and $\theta_{i,k}$, respectively, with their precise interpretation again depending on the choice of $\tau$, the link functions and covariates included in the models.

From a practical perspective, that $\alpha$ = ($\mbox{$\alpha$}_{1}$, $\mbox{$\alpha$}_{2}$, $\mbox{$\alpha$}_{\theta}$) consists of 3$K$ parameters may result in computational and/or convergence issues unless the observed data is rich (i.e. a large sample size or high event rates) or the initial partition is coarse. To mitigate such issues we propose that a B-spline structure be adopted across the $K$ components of each of $\mbox{$\alpha$}_{1}$, $\mbox{$\alpha$}_{2}$ and $\mbox{$\alpha$}_{\theta}$ [8]. Towards the specification for $\mbox{$\alpha$}_{1}$, let $\{\tilde{t}_{1}^{1},\ldots,\tilde{t}_{1}^{J_{1}}\}$ denote a collection of $J_{1}$ user-specified knots on the interior of the range ($\tau_{0},\tau_{K}$) and $q_{1}$ the user-specified degree of the local polynomial basis functions. Given these choices, we specify the $k^{th}$ component of $\mbox{$\alpha$}_{1}$ as $\alpha_{1,k}=\sum_{j=1}^{\widetilde{J}_{1}}\eta_{1,j}B_{j,k},$ where $B_{j,k}$ is the value of the $j^{th}$ B-spline at time $\tau_{k}$ and $\widetilde{J}_{1}=J_{1}-1+q_{1}$ is total number of spline terms. Thus, this specification requires estimation of $\widetilde{J}_{1}$ coefficient terms, specifically $\mbox{$\eta$}_{1}=(\eta_{1,1},\ldots,\eta_{1,\widetilde{J}_{1}})$. Applying the same strategy to $\mbox{$\alpha$}_{2}$ and $\mbox{$\alpha$}_{\theta}$ will result in $\widetilde{J}_{1}+\widetilde{J}_{2}+\widetilde{J}_{\theta}$ unknown coefficients, which we collectively denote as $\eta$ = ($\mbox{$\eta$}_{1}$, $\mbox{$\eta$}_{2}$, $\mbox{$\eta$}_{\theta}$). Finally, to distinguish between specifications, in the remainder of this paper we refer to the model based on $\alpha$ with 3$K$ unknowns as the unstructured model while that based on $\eta$ with $\widetilde{J}_{1}+\widetilde{J}_{2}+\widetilde{J}_{\theta}$ unknowns as the B-spline model.

The central innovation of the framework in expressions (5)-(7) is that dependence between the non-terminal and terminal events is quantified in two distinct and yet complimentary ways. The first is captured through $\theta_{i,k}$ which can be viewed as a measure of local dependence in that it quantifies the risk of co-occurrence of the two events, via the odds-ratio [31, 27, 4, 35, 24, 29], during the $(\tau_{k-1},\tau_{k}]$ interval; a large positive value of $\theta_{i,k}$ indicates that if the non-terminal event (e.g. AD) occurs during $(\tau_{k-1},\tau_{k}]$ then the terminal event (e.g. death) is likely to subsequently occur during the same interval. Because of the novel parameterization, we emphasize that the local dependence can vary over the time scale or as a function of covariates. As such, one could investigate whether local dependence between AD and death is weaker at younger ages relative to older ages. Furthermore, one could investigate whether the probability that the two events co-occur changes in response to key life events such as the death of a partner.

The second quantification of dependence is through the components of $\mbox{$\beta$}_{2}$ in expression (6) that correspond to how the status of the non-terminal event, $Y_{1i,k-1}$, influences $\pi_{2i,k}$. Labelling these components as $\mbox{$\beta$}_{2,y}$, analogous to how one interprets covariates effects in regression models, these parameters capture the extent to which whether the non-terminal event has occurred is associated with a change in risk of the terminal event; for this reason, we refer to $\mbox{$\beta$}_{2,y}$ as capturing global dependence. Note, global dependence is conceptually similar to the explanatory hazard ratio (EHR) and cross-quantile residual ratio (CQRR) [41] in that it concerns the role that the non-terminal event plays in modifying the future occurrence of the terminal event (see Section A.1 of the Supplementary Materials).

The partition, $\tau$, plays a critical role in the proposed framework in that it provides the foundation for being able to distinguish between local and global dependence, as we have conceptualized them, and for being able to investigate the role that covariates play. In Section A.4 of the Supplementary Materials we show that a distribution for $(T_{1},T_{2})$ induces the quantities (5)–(7) for any $\tau$. Thus, regardless of the choice of $\tau$, the components of the proposed model are well-defined mathematical objects and are, therefore, valid targets for estimation and inference. An important consequence of this is that one cannot say that any given partition corresponds to the ‘truth’. As indicated in Section 3.2, however, the choice of $\tau$ dictates the numerical values and interpretation of the quantities given by (2)-(4), and correspondingly (in part, at least) the numerical values and interpretation of the parameters in the regression structure given by expressions (5)-(7). As such the choice of $\tau$ is a critical challenge that requires careful consideration by the study investigators.

In principle, one could approach choosing $\tau$ through consideration of the clinical condition under investigation such as the pace at which the disease progresses. In the ACT study, for example, the choice to schedule follow-up biennially was made in part for logistical considerations (it would be challenging to follow-up approximately 4,000 participants each year) but also because AD is a slowly-developing condition. Alternatively, one may pursue a data-driven strategy where, for example, some goodness-of-fit criterion is specified and then optimized as a function of $\tau$. Our perspective is that the decision should be based primarily, if not exclusively, on clinical considerations. Central to this position is that, in addition to their interpretation, the numerical values of (2)-(4) change with $\tau$. To see this, we again note that the probabilities in each of (2)-(4) speak to the cumulative incidence of events during the interval $(\tau_{k-1},\tau_{k}]$. If the length of the interval is decreased then the incidence will necessarily decrease. The corresponding new interpretation and numerical value will not be ‘wrong’, however, but just different. Put another way, the change in the interpretation and the numerical values of the model parameters that result from, say, adopting a finer partition of the time scale should, arguably, be viewed as a change in the question that is being answered. Thus, in our view, purely data-driven approaches, while they may have some initial intuitive appeal, should be avoided.

Nevertheless, we do acknowledge that it may not be easy to elicit a single partition, from the literature or collaborators, on which to base the analyses and conclusions. If it is the case that there is no clear choice, analysts may opt to perform a range of analyses over different partitions, both in terms of how fine the partition is and in terms of where the cut-points are for a given (common) interval length. We pursue this strategy in conducting the analyses of the data from ACT in Section 6.

Building on the notation developed in Section 3, the first two columns of Table 1 provide the six possible outcome data scenarios in the $k^{th}$ interval of the partition given by $\tau$ as a function of the outcome vector in the previous interval (i.e. as a function of ($Y_{1i,k-1},Y_{2i,k-1}$)). The third column provides the corresponding likelihood contributions, that is the interval-specific components in the decomposition given by expression (1), with

	$\displaystyle\pi_{12i,k}$	$\displaystyle=$	$\displaystyle\mbox{P}(Y_{1i,k}=1,Y_{2i,k}=1|Y_{1i,k-1}=0,Y_{2i,k-1}=0,\overline{\mbox{$X$}}_{i,k})$
		$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{cc}\pi_{1i,k}\pi_{2i,k}(0)&\hskip 14.45377pt\theta_{i,k}=1\\ \frac{1}{{2(\theta_{i,k}-1)}}\times\left[1+a_{ij}-\sqrt{(1+a_{i,k})^{2}-4\theta_{i,k}(\theta_{i,k}-1)\pi_{1i,k}\pi_{2i,k}(0)}\right]&\hskip 14.45377pt\theta_{i,k}\neq 1\end{array}\right.,$

where $a_{i,k}=(\pi_{1i,k}+\pi_{2i,k}(0))(\theta_{i,k}-1)$ [35].

Let $\boldsymbol{\phi}$ denote the collection of unknown parameters in the specification of model (5)-(7). Note, if the unstructured form of the model is fit then $\boldsymbol{\phi}\equiv\boldsymbol{\phi}^{\boldsymbol{\alpha}}=(\boldsymbol{\alpha},\boldsymbol{\beta})$ while $\boldsymbol{\phi}\equiv\boldsymbol{\phi}^{\mbox{$\eta$}}=(\mbox{$\eta$},\boldsymbol{\beta})$ if the B-spline model is fit. For either specification, the observed data likelihood for a random sample of $N$ study participants from the population of interest, $\mathcal{L}(\boldsymbol{\phi})$ is the product of $N$ terms, each of the form:

	$\displaystyle\mathcal{L}_{i}(\boldsymbol{\phi})\$	$\displaystyle=\ \prod_{k=k^{l}_{i}}^{k^{r}_{i}}\mbox{P}(Y_{1i,k}=y_{1i,k},Y_{2i,k}=y_{2i,k}|\ Y_{1i,k-1}=y_{1i,k-1},Y_{2i,k-1}=y_{2i,k-1},\overline{\mbox{$X$}}_{i,k}),$
		$\displaystyle=\ \prod_{k=k^{l}_{i}}^{k^{r}_{i}}\Big{\{}[\pi_{12i,k}]^{y_{1i,k}y_{2i,k}}[\pi_{1i,k}-\pi_{12i,k}]^{y_{1i,k}(1-y_{2i,k})}[\pi_{2i,k}(0)-\pi_{12i,k}]^{(1-y_{1i,k})y_{2i,k}}$
		$\displaystyle\hskip 36.135pt\times\ [1-\pi_{1i,k}-\pi_{2i,k}(0)+\pi_{12i,k}]^{(1-y_{1i,k})(1-y_{2i,k})}\Big{\}}^{(1-y_{1i,k-1})(1-y_{2i,k-1})}$
		$\displaystyle\hskip 36.135pt\times\Big{\{}[\pi_{2i,k}(1)]^{y_{2i,k}}[1-\pi_{2i,k}(1)]^{1-y_{2i,k}}\Big{\}}^{y_{1i,k-1}(1-y_{2i,k-1})}.$

If the unstructured form of model (5)-(7) is adopted, estimation can proceed in a straightforward manner by finding the maximizer of $\ell(\boldsymbol{\phi}^{\boldsymbol{\alpha}})$ = $\log\mathcal{L}(\boldsymbol{\phi}^{\boldsymbol{\alpha}})$, denoted by $\widehat{\boldsymbol{\phi}}^{\boldsymbol{\alpha}}$. Under the proposed B-spline model, however, care is needed to avoid overfitting. To resolve this, we maximize the observed data likelihood subject to a penalty that imposes smoothness in resulting estimated functions. One such penalty is the integrated squared second derivative which, following Eilers et al. [8], can be approximated by penalizing coefficient differences. Let $\Delta$ be the difference operator so that, for example, $\Delta\eta_{1,j}=\eta_{1,j}-\eta_{1,j-1}$. Then, letting $\boldsymbol{\lambda}=(\lambda_{1},\lambda_{2},\lambda_{\theta})$, consider the penalized likelihood:

$$\mathcal{L}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})\ =\ \mathcal{L}(\boldsymbol{\phi}^{\mbox{$\eta$}})\ -\ \lambda_{1}\sum_{j=m+1}^{\widetilde{J}_{1}}(\Delta^{m}\eta_{1,j})^{2}\ -\ \lambda_{2}\sum_{j=m+1}^{\widetilde{J}_{2}}(\Delta^{m}\eta_{2,j})^{2}\ -\ \lambda_{\theta}\sum_{j=m+1}^{\widetilde{J}_{\theta}}(\Delta^{m}\eta_{\theta,j})^{2},$$

where $\Delta^{m}$ corresponds to the difference operator having been applied $m$ times; for example, $\Delta^{2}\eta_{1,j}=\Delta(\Delta\eta_{1,j})=\eta_{1,j}-2\eta_{1,j-1}+\eta_{1,j-2}$. Note, letting $\mbox{$D$}_{m}$ denote the matrix representation of $\Delta^{m}$, such that, for example, $\sum_{j=m+1}^{\widetilde{J}_{1}}(\Delta^{m}\eta_{1,j})^{2}$ can be written as $\mbox{$\eta$}_{1}^{T}\mbox{$D$}_{m}^{T}\mbox{$D$}_{m}\mbox{$\eta$}_{1}$ (Section A.4 of the Supplementary Materials), the penalized likelihood can be written as:

$$\mathcal{L}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})\ =\ \mathcal{L}(\boldsymbol{\phi}^{\mbox{$\eta$}})\ -\ \lambda_{1}\mbox{$\eta$}_{1}^{T}\mbox{$D$}_{m}^{T}\mbox{$D$}_{m}\mbox{$\eta$}_{1}\ -\ \lambda_{2}\mbox{$\eta$}_{2}^{T}\mbox{$D$}_{m}^{T}\mbox{$D$}_{m}\mbox{$\eta$}_{2}\ -\ \lambda_{\theta}\mbox{$\eta$}_{\theta}^{T}\mbox{$D$}_{m}^{T}\mbox{$D$}_{m}\mbox{$\eta$}_{\theta}.$$

(9)

Let $\ell(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})=\log\mathcal{L}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})$ be the log-penalized likelihood. Furthermore, let $\boldsymbol{U}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})=\nabla_{\boldsymbol{\phi}^{\mbox{$\eta$}}}\ell(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})$ be the gradient of $\ell(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})$ with respect to $\boldsymbol{\phi}^{\mbox{$\eta$}}$ and $\mathcal{H}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})=\nabla_{\boldsymbol{\phi}\boldsymbol{\phi}^{\mbox{$\eta$}}}\ell(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})$ the corresponding matrix of second partial derivatives (i.e. the Hessian matrix). For a given value of $\boldsymbol{\lambda}$, the penalized maximum likelihood estimator, which we denote by $\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}}_{\boldsymbol{\lambda}}$, is the solution to $\boldsymbol{U}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})=\mbox{$0$}$. Note, since the penalty is quadratic in $\mbox{$\eta$}_{1}$, $\mbox{$\eta$}_{2}$, and $\mbox{$\eta$}_{\theta}$, gradient-based maximization is relatively straightforward to carry out.

Finally, towards choosing the value of $\boldsymbol{\lambda}$ on which the final results are to be based, say $\boldsymbol{\lambda}^{*}$, one can proceed using any standard model-selection criteria such as the Akaike Information Criterion (AIC) or cross validation [11, 8]. In our implementation of the methods, we follow Gray [12] by taking the trace of $\mathcal{H}(\boldsymbol{\phi}^{\mbox{$\eta$}};\mbox{$0$})\mathcal{H}^{-1}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})$ as the effective degrees of freedom when calculating the AIC.

For the unstructured model, under standard regularity conditions, the asymptotic distribution for the maximum likelihood estimator $\widehat{\boldsymbol{\phi}}^{\boldsymbol{\alpha}}$ follows from standard likelihood theory. That is, $\widehat{\boldsymbol{\phi}}^{\boldsymbol{\alpha}}$ is consistent and, letting $\widetilde{\boldsymbol{\phi}}^{\boldsymbol{\alpha}}$ denote the value of ${\boldsymbol{\phi}}^{\boldsymbol{\alpha}}$ induced from the partition, $N^{1/2}(\widehat{\boldsymbol{\phi}}^{\boldsymbol{\alpha}}-\widetilde{\boldsymbol{\phi}}^{\boldsymbol{\alpha}})$ is asymptotically normal with the variance being the usual inverse of Fisher information matrix.

For the B-spline model, careful consideration of the asymptotic properties of $\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}}$ requires specification of how the number of knots, $\widetilde{\mbox{$J$}}=(\widetilde{J}_{1},\widetilde{J}_{2},\widetilde{J}_{\theta})$, and the penalty, $\boldsymbol{\lambda}^{*}$, change with the sample size $N$. This has been the focus of a large body of work, much of which is concerned with error estimation for the non-parametric function(s) [e.g., 26, 18, 5]. Here, however, primary interest lies with statistical inference for the parameters $\boldsymbol{\phi}^{\mbox{$\eta$}}$. Given this, and as recommended by number of authors for development of practical inference in a range of setting [11, 12, 39, 42], we consider the behavior of $\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}}$ in settings where both $\widetilde{\mbox{$J$}}$ and $\boldsymbol{\lambda}$ are fixed. With this, henceforth, we denote the estimator as $\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}}_{\boldsymbol{\lambda}}$.

To this end, we show in Section A.5 of the Supplementary Materials, that under standard regularity assumptions $N^{1/2}(\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}}_{\boldsymbol{\lambda}}-\widetilde{\boldsymbol{\phi}}^{\mbox{$\eta$}}_{\boldsymbol{\lambda}})$ is asymptotically normal with variance that can be consistently estimated by

$$\widehat{\mathcal{V}}\ =\ N\mathcal{H}^{-1}(\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}};\boldsymbol{\lambda})\left(\frac{1}{N}\sum_{i=1}^{N}\boldsymbol{U}_{i}(\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}};\boldsymbol{\lambda})\boldsymbol{U}_{i}^{T}(\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}};\boldsymbol{\lambda})\right)N\mathcal{H}^{-1}(\widehat{\boldsymbol{\phi}}^{\mbox{$\eta$}};\boldsymbol{\lambda}),$$

(10)

where $\widetilde{\boldsymbol{\phi}}^{\mbox{$\eta$}}_{\boldsymbol{\lambda}}$ is the solution of $E[\boldsymbol{U}(\boldsymbol{\phi}^{\mbox{$\eta$}};\boldsymbol{\lambda})]=\mbox{$0$}$ and where $\boldsymbol{U}_{i}(\cdot;\boldsymbol{\lambda})$ is the gradient of the log penalized likelihood of a single observation $i$. For individual parameters, such as the coefficient of APOE in the model for $\theta$, the appropriate entry in the diagonal of (10) is used for estimating the variance. For the baseline time trends, let $\widehat{\mathcal{V}}_{\mbox{$\eta$}_{1}}$ be the sub-matrix corresponding to the estimated variance matrix of $\widehat{\mbox{$\eta$}}_{1}$. Note that $\widehat{\boldsymbol{\alpha}}_{1}$ can be written as $\widehat{\boldsymbol{\alpha}}_{1}=\widehat{\mbox{$\eta$}}^{T}_{1}\mathcal{B}$ with $\mathcal{B}$ being a $\widetilde{J}_{1}\times K$ matrix with elements $\mathcal{B}_{jk}$=$B_{j,k}$. Therefore, to construct a 95% confidence interval for the AD time-varying reference probability, $\pi_{1i,k}$, one can use $g^{-1}_{1}\{\widehat{\mbox{$\eta$}}^{1}\pm 1.96[diag(\mathcal{B}^{T}\widehat{\mathcal{V}}_{\mbox{$\eta$}_{1}}\mathcal{B})]^{1/2}\}$ where here $g^{-1}$ to be understood as operating entrywise on the vector, and $diag(\cdot)$ returns the diagonal of a matrix.

Prior to presenting results based on the proposed framework, we consider a series of analyses based on existing methods described in Section A.1 of the Supplementary Materials. Results are given in Figure 2, and in Section A.6 of the Supplementary Materials.

Finally, we report results from a series of analyses using the proposed longitudinal bivariate modeling framework. As emphasized in Section 3.5, the choice of the partition $\tau$ is important. Towards investigating the role of the choice, and recalling from Section 2 that participants underwent biennial visits, we considered two partitions of the time scale $[65,100)$: $\mbox{$\tau$}^{2.5}=\{65.0,67.5,70.0,\ldots,97.5,100.0\}$, for which $K$=14, and $\mbox{$\tau$}^{5.0}=\{65.0,70.0,75.0\ldots,95.0,100.0\}$, for which $K$=7. Tables A.2 and A.3 in the Supplementary Materials provide the $2\times 2$ outcome tables for each interval, under the two participations.

Towards specification of the models in expressions (5)-(7), we used logit links for $\pi_{1i,k}$ and $\pi_{2i,k}(t_{1})$, so that $g_{1}^{-1}(\cdot)=g_{2}^{-1}(\cdot)=\mbox{expit}(\cdot)$. For the cross-sectional odds ratio, $\theta_{i,k}$, we used a log link, so that $g_{\theta}^{-1}(\cdot)=\mbox{exp}(\cdot)$. Based on these, we considered unstructured and B-spline specifications for the three sets of baseline parameters, $\mbox{$\alpha$}_{1}$, $\mbox{$\alpha$}_{2}$ and $\mbox{$\alpha$}_{\theta}$. For the latter, we used 10 knots for $\mbox{$\tau$}^{2.5}$ and 5 knots for $\mbox{$\tau$}^{5}$, with a cubic spline and a second order difference penalty (i.e. $m=2$ in $\Delta^{m}$), and considered varying degrees of penalty, specifically setting $\lambda_{1}$ = $\lambda_{2}$ = $\lambda_{\theta}$ = $\lambda$, with $\lambda\in\{0.0,0.1,0.5,1.0,2.5,5.0\}$. Table A.17 in the Supplementary Materials reports AIC from the fitted models, across the different $\lambda$.

Finally, $\overline{\mbox{$X$}}_{i1,k}$, $\overline{\mbox{$X$}}_{i2,k}$, and $\overline{\mbox{$X$}}_{i\theta,k}$ were specified with the overarching goal of assessing the joint impact of having $\geq$ 1 APOE $\epsilon$4 allele and gender on the joint risk of AD and death. As such, the models for $\pi_{1i,k}$ and $\theta_{i,k}$ included main effects for having $\geq$ 1 APOE $\epsilon$4 allele and gender, and their interaction, while the model for $\pi_{2i,k}(\cdot)$ included main effects, two-way interactions and the three-way interaction between AD diagnosis, having $\geq$ 1 APOE $\epsilon$4 allele and gender. For $\overline{\mbox{$X$}}_{i1,k}$, $\overline{\mbox{$X$}}_{i2,k}$, we additionally included all other variables available in dataset (i.e. race/ethnicity, marital status, education and depression).