Joint modeling of geometric features of longitudinal process and discrete survival time measured on nested timescales: an application to fecundity studies

In biomedical studies, information is routinely collected longitudinally on various biomarkers up to a time-to-event (usually censored), along with additional covariates. An often cited example of such data is from HIV clinical trial with the key longitudinal process of interest being CD4 counts. Many other examples from cancer studies, reproductive health etc. have been the motivation for the development of various methods in this area of research. In joint modeling, one is typically interested in (i) how to model the pattern of change of the longitudinal process, and (ii) to characterize the relationship between the survival event, the longitudinal process, and the covariates.

Most of the literature on joint modeling of longitudinal process and time-to-event have focused on modeling the mean of the longitudinal process, where the dependence between the longitudinal process and time-to-event is through the mean process with subject-specific random effect(s) (Wu:Carroll:1988, ; Tsiatis:Davidian:2004, ). Joint modeling in the context of longitudinal process and discrete time-to-events was first studied by Albert et al. (2010) albert2010approach , and later by Qiu et al. (2016)qiu2016joint in the context of shared parameter framework approach. Recently, some authors have also considered modeling the dispersions of the longitudinal process and time-to-event (gao2011joint, ; Mclain:Lum:Sundaram:2012, ). In many situations (the motivating example discussed below), it may be more of interest to study the geometric features of a cyclical longitudinal process. Another common aspect of the existing joint modeling literature is that they focus on the situations where the longitudinal process and the time-to-event are on the same time scale, while there are examples where this assumption may not be true. For instance, in our motivating example, the longitudinal process is measured on a daily scale while the time-to-pregnancy is measured in menstrual cycles (approximately monthly). In this paper, our main objective is to jointly model cyclic longitudinal process with time-to-event, measured on nested timescales, where the dependence between the cyclic longitudinal process and time to event is captured through cycle-level geometric features of the longitudinal process with subject-specific random effects.

We now present our motivating example from reproductive health studies. Reproductive hormones, like the luteinizing hormone (LH), estrogen and its metabolites (eg, estrone-3-gluconoride (e3g) etc) patterns play an important role in the study of conception, infertility and other chronic disease (Mumford:etal:2012, ; Parazzini:etal:1993, ; Terry:etal:2005, ; Whelan:etal:1994, ; Solomon:etal:2001, ; Solomon:etal:2002, ). However, data on these endogenous hormones is difficult to quantify due to complex cyclical patterns of hormones, the need for timed collection, and the cost required for multiple sample collections. So, much of research has focused on menstrual cycle characteristics such as cycle length as proxies for cumulative hormonal exposure and/or hormonal patterns as they can be easily assessed in population studies. Short and long or irregular cycles have been associated with increased risks for breast cancers, osteoporosis, type 2 diabetes mellitus, cardiovascular diseases etc. Consequently, it is important to be able to study the hormonal profiles directly. LH has important role in the luteinizing of the follicle and the functional maturation of the nucleus of the oocyte (Verpoest:etal:2000, ). Abnormalities in the LH surge may impair the development of the oocyte and consequently it’s fertilization ability. Additionally, LH surge abnormalities, such as reduced peak values, have been associated with infertility (Cohlen:etal:1993, ; Cahill:etal:1997, ), indicating that shape of the hormonal curve plays an important biological role as well. Moreover, the pattern of the LH profile is highly variable in normal menstruating women. Motivated by these issues, we are interested in modeling the LH surge and its relationship with time-to-pregnancy. A source of such data arises in the prospective pregnancy studies, where couples are followed from the time they go off contraception (with the intention of becoming pregnant) until they get pregnant. In these studies, women use ovulation kits which track the LH and e3g daily within a fixed window to precisely capture the impending day of ovulation in each menstrual cycle. This provides an opportunity to study the patterns of surge of these hormones, especially LH and its relationship with time-to-pregnancy. Motivated by scietific hypothesis, we focus on three geometric features, namely, the value of LH peak, the curvature at the LH peak and the average curvature of the LH profile within fertile window (the window of opportunity around ovulation for conception in a menstrual cycle) from the longitudinally measured LH values within each menstrual cycle.

Motivated by the example described above, we consider joint modeling of the cyclic longitudinal process (“hormonal profile”) and a discrete survival time (“time-to-pregnancy”), where one is interested in modeling various geometric features of the longitudinal process at the cycle-level (e.g., value at the peak, curvature at the peak, average curvature within fertile window of a cycle). In the next section, we introduce our data and the modeling framework. In Section 3, we provide the estimation approach for the parameters of interest and also discuss the prediction approach for the time-to-event distribution of a new subject given its longitudinal measurement history. We assess the performance of the proposed estimates through simulation studies in Section 4. In Section 5, we present detailed analysis of the Stress and Time-to-pregnancy, a sub-component of the Oxford Conception Study (louis2011stress, ; McLain:Sundaram:Louis:2015, ). .

For woman $i$, $i=1,\dots,n$, let $T_{i}$ denote the time-to-pregnancy, i.e., the number of menstrual cycles it took to get pregnant. As is typical in time-to-event studies, $T_{i}$ is subject to right censoring and one observes $X_{i}=\min(T_{i},\tau_{i})$ and $\delta_{i}=I(T_{i}\leq\tau_{i})$, where $\tau_{i}$ denotes the censoring time and $I(\cdot)$ denotes the indicator function. Throughout this article, we assume that the censoring time $\tau_{i}$ is independent of $T_{i}$.

In prospective pregnancy studies, ovulation kits are used by each woman participant to identify the day of ovulation. These kits typically require testing to be done on day $6$ through day $25$ of every menstrual cycle for measuring the hormone levels to identify the day of ovulation within each cycle. So, essentially the test results in underlying hormonal data on LH from day $6$ through day $25$ in each cycle for every woman. We denote the hormonal data in a cycle as $h(t)$, $t\in[L,M]$ where $L,M$ are days from the start of cycle where hormone measurements are collected within each menstrual cycle and $\hat{t}$ is the time point from the start of the cycle where the hormonal profile reaches the highest level (peak). Also, we define the fertile window “the window of opportunity for conception around ovulation” in a cycle as $[\hat{t}-5,\ \hat{t}+1].$ The three geometric features of hormonal profile that are of interest here are: (1) cycle-specific curvature at the hormone peak, $k(\hat{t})=|h^{\prime\prime}(\hat{t})|\{1+h^{\prime}(\hat{t})\}^{-3/2}$, which measures the sharpness of the hormonal profile at the peak time; (2) cycle-specific hormone peak value, $h(\hat{t})$; (3) the average curvature of the hormone profile within fertile window, $\sum_{t=\hat{t}-5}^{\hat{t}+1}k(t)/6$; within each cycle.

Let $\tilde{Y}_{ij}$ be the true geometric feature of hormonal profile of interest for $i$-th woman in menstrual cycle $j$ and let $Z_{ij}$ be the vector of covariates for $\tilde{Y}_{ij}$. Denote by $Y_{ij}$ the observed geometric feature of hormonal profile for $i$-th woman in menstrual cycle $j$. The observed geometric features for each cycle are then calculated from the observed hormone data, first by using some smoothing technique, e.g. B-splines and then calculating the feature based on the formula mentioned in the previous paragraph.

We relate the true geometric feature $\tilde{Y}_{ij}$ to the covariates $Z_{ij}$ through a linear model with subject-specific random intercept $b_{Y,i}$,

where $\bm{\beta}$ is the vector of regression coefficients corresponding to $Z_{ij}$. The observed hormonal geometric features $Y_{ij}$ are then modeled by

where $\epsilon_{ij}$ are all independent and identically distributed and follow $N(0,\sigma^{2})$ distribution.

We use the discrete survival model by Sundaram et al. (2012) Sundaram2:McLain:Louis:2012 for TTP where the hazard of discrete survival time is related linearly to the covariates when transformed by a complementary log-log function. It also accounts for the fact that, the hazard for conception in a cycle is zero if the couple does not have any intercourse in the fertile window of that cycle. The model is given by

where for subject $i$ in cycle $j$, $U_{ij}$ is the vector of covariates for TTP (which could have overlap with $Z_{ij}$, the covariates for the geometric features), $A_{ij}$ is the indicator of intercourse within fertile window, $\rho_{j}$ is the baseline effect for cycle $j$, $\bm{\gamma}$ represents the regression coefficients of $U_{ij}$ and $b_{T,i}$ is a subject-specific random effect. Recall that the fertile window refers to the days in a menstrual cycle around the day of ovulation when a couple having intercourse can potentially conceive; $A_{ij}=0$ means that couple $i$ did not have intercourse during the fertile window of cycle $j$, which implies that there is no risk of pregnancy in that cycle.

To study the association between a woman’s hormonal profile and TTP, we take into account the cycle-level geometric feature of a woman’s hormonal profile, $\tilde{Y}_{ij}$, in (2). Recalling that $\tilde{Y}_{ij}=Z_{ij}^{\prime}\bm{\beta}+b_{Y,i}$, we propose the following model

where $\psi_{\mu}$ is the regression coefficient of $\tilde{Y}_{ij}$, and assume that $\mathbf{b}_{i}\equiv(b_{Y,i},b_{T,i})$ follows multivariate normal distribution with mean $\mathbf{0}$ and variance covariance matrix $D$, MVN$(\mathbf{0},D)$, with

and $\mathbf{b}_{i}$’s are independent and identically distributed (iid) and independent of $\epsilon_{ij}$. The association between the hormonal profile and TTP is taken into account not only through the fact that the cycle-level geometric feature is included in the model for TTP as a predictor, but also through the possible correlation between subject specific random effects $b_{Y,i}$ and $b_{T_{i}}$.

In this section, we conducted simulation studies to investigate the performance of the proposed estimates using likelihood-based approach under various settings. For simplicity, we assumed covariates $U_{ij}=U_{i}$ and generated $U_{i}$ from $\textrm{N}(2,1)$. Let $Z_{ij}=Z_{i}=(1,U_{i})$. $\mathbf{b}_{i}=(b_{Y,i},b_{T,i})$ were generated from MVN$(\mathbf{0},D)$,

The random error $\epsilon_{ij}$ were generated from $\textrm{N}(0,0.9^{2})$ and $Y_{ij}$ were generated using (1). Similarly, TTP for $i$-th subject is simulated using (2). True values of $\bm{\beta},\bm{\gamma},\psi,\sigma_{1},\sigma_{2}$, and $\zeta$ are listed in Table 2. Subjects who had not experienced an event till $j=6$ were censored.

We have considered 9 different simulation settings with varying sample sizes ($n$= 300, 400, and 500) as well as, with varying intercourse success probabilities ($p(A_{ij})$= 0.95,0.90 and 0.85). These choices are made in such a way that a particular case ($n=400$, $p(A_{ij})=0.95$) closely mimics the real data in terms of distribution of TTP while the other cases either overestimate or under-estimate the number of cycles and sample size in comparison to real data. For each simulation scenario, $m=$1000 replicates are generated. We split each replicate data into a training set (randomly selected 2/3rd) and a test set, as is originally only done for real data. Table 1 provides a summary of TTP distribution (both overall as well as in the training set). The shaded row indicates the simulation scenario that mimics the real data in terms of sample size ($337$ in real data), # of cycles ($1023$ in real data), intercourse probability ($\approx 94.7\%$ in real data) and # of cycles in training data ($686$ in real data). In Table 1 we have provided avg. TTP, proportion of right censoring ($\hat{p}(\delta=0)$), # of cycles in overall as well as in training sets, each averaged (mean(Mn)/median(Med)) over 1000 replicates (provided with interquartile range (IQR)).

Estimate accuracy with the training sets: For each replicate data, training set is used to fit the model, while test set is used to estimate predictive accuracy. Based on training set, we have reported estimation bias (Bias), standard deviation (SD) and coverage probability (CP) based on 95% CI of the parameter estimate with varying samples sizes ($n$= 300, 400, and 500) for a given intercourse success probability, p($A_{ij}$)= 0.95, in Table 2. Same set of metrics are reported for p($A_{ij}$)= 0.90 and 0.85 in Appendix A, Table 7, and Table 8 respectively. $n_{tr}$ represents the training set size, which is 2/3rd of the sampled data. Estimation used Gaussian quadrature with 50 nodes and the simulation was conducted in software R/4.2.

The simulation results show that the proposed estimation approach works well in all settings considered here with reasonably small bias and coverage probabilities close to the nominal level. For $\rho_{3}$, $\sigma_{1}$, and $\zeta$, CP falls slightly below 95% relative to other parameters, which may be attributable to numerical approximation of marginal likelihood.

We illustrate our proposed method by analyzing the Oxford Conception Study, Stress and Time-to-Pregnancy component (Pyper:etal:2006, ), which is a prospective cohort study with preconception enrollment of women aged 18-40 years who were attempting to become pregnant. The women in the Oxford Conception Study (hereafter OCS) provided daily level information on reproductive hormones, in particular, luteinizing hormone during mid-cycle, intercourse acts, and host of other lifestyle variables, along with couple level baseline covariates. The lutenizing hormone was observed through the ovulation kits used by the woman to track her daily fertility level to identify ovulation, the values provided by the monitor were monotonically transformed values of the lutenizing hormones (manufacturer only provides these transformed values for research purpose). These women were prospectively followed for the number of menstrual cycles it took them to become pregnant (i.e., human-chorionic gonadotropin confirmed pregnancy on the day of expected menses), i.e., TTP, or a maximum of $6$ menstrual cycles. Other examples of prospective pregnancy cohort studies’ include the Longitudinal Investigation of Fertility and Environment (Louis:etal:2011, ), Fertili (Colombo:Masarotto:2000, ), and Billings (Colombo:etal:2006, ).

The data consisted of 337 women with LH measurements, resulting in a total of 1023 menstrual cycles. We randomly selected 225 (about two-thirds) of the women with 686 cycles in the training set and the rest of the women were taken to be the prediction set. The hormonal measurements were taken daily within a fixed window in each menstrual cycle. For each cycle of every woman, we used B-spline functions to smooth the observed hormone data and then calculate the geometric features of interest. As mentioned in previous sections, we considered three geometric features of the hormone profiles: curvature at the LH peak which measures the sharpness of the LH profile at peak, LH peak value and the average LH profile curvature within fertile window. Figure 2 gives an example of the LH measurements for 4 randomly selected women in the data set. Note the number of varying cycles worth of information based on TTP for each woman in the figure. Motivated by the underlying biological hypothesis discussed in detail in the introduction, we fit three different joint models corresponding to each of the geometric feature. In each joint model, we model TTP with each one of the cycle-level geometric features using the training dataset and evaluate their prediction abilities, respectively. In our analysis, we focus on the following covariates: female age (in model for hormonal data), couple average age and difference between female and male age (in model for TTP), and female’s body mass index (BMI) (categorized as underweight or normal weight if BMI$<25$; overweight if $25\leq$BMI$<30$; obese if $30\leq$BMI)), female’s smoking category (smoke.n if not smoking; smoke.m if smoking but average cigarette smoked per day $\leq 10$; smoke.h if average cigarette smoked per day $>10$), stress level measured by the salivary biomarker alpha amylase and parity (nulliparous or multiparous; nulliparous women are those who haven’t had live birth before) in both models. The three age related covariates as well as alpha amylase level were scaled by suitable constants in the analysis and were denoted as Age*, avg.age*, dif.age* and Alp*, respectively. The geometric features of hormonal profiles were also scaled by adequate constants in the analysis. Subjects with BMI smaller than 25 who are non-smokers and nulliparous were considered as the reference group.

Table 4 presents the estimates and $95\%$ confidence intervals for the covariate effects on the curvature at LH peak, covariate effects on TTP, cycle-specific baseline effect on TTP, effect of curvature at LH peak on TTP as well as the variance covariance parameters for the random effects. Notice that all covariates seem non-significant for the curvature at LH peak. A couple’s average age and female-male age difference seem to have significant negative effect on TTP. Women who smoke but smoke no more than 10 cigarettes per day have significantly longer TTP compared to non-smoking women while the difference between the heavy smoking women ($>10$ cigarettes per day) and non-smoking women was not significant. Multiparous women seem to have significantly shorter TTP than nulliparous women. We also found that the curvature at LH peak has significant positive effect on TTP. This implies that the sharper a woman’s LH peak is, the shorter her TTP tends to be.

We also present the estimation results for jointly modeling LH peak value and TTP in Table 5. The conclusions of covariate effects were generally consistent with that in Table 4, except that obesity seems to have protective effect on TTP and heavy smoking women also seem to have significantly longer TTP than non-smoking women.

The estimation results for joint modeling of average curvature of LH profile within fertile window and TTP are shown in Table 6 The overall trend of covariate effects is consistent with that in Table 4 for curvature at LH peak, except that overweight women seem to have significant longer TTP than underweight/normal weight women, women with higher stress level (alpha amylase) seem to have longer TTP, while the effect of smoking on TTP does not seem to be significant.

Now that we have fitted the joint models using the training dataset, we are interested in predicting the probability of subfertility (i.e., TTP $>6$ cycles) given survival past one cycle. That is, for all women in the prediction set with $T_{i}>1$, we wish to classify $I(T_{i}>6)$ by the conditional survival probability $\pi_{i}(6\mid 1)$. To measure the classification rate, we empirically estimated the sensitivity: $P(\hat{\pi}_{i}(6\mid 1)>c\mid T_{i}>6)$, and specificity $P(\hat{\pi}_{i}(6\mid 1)\leq c\mid T_{i}\leq 6)$, where $\hat{\pi}_{i}(6\mid 1)$ is the mean of the Monte Carlo sample $\{\pi_{i}^{(l)},l=1,\dots,L\}$ obtained following the steps in Section 2.1. Of the 79 women available for prediction, 20 were censored between 1 and 6 cycles and removed from the prediction analysis. The classification measures using ROC for the models with curvature at LH peak, LH peak value and average curvature of LH profile within fertile window, are displayed in Figure 3(a),(b),(c), respectively, for all $c\in[0,1]$. The AUC is 0.650 for the model which includes the curvature at LH peak, 0.646 for the model which includes the LH peak value and 0.614 for the model which includes the average curvature of LH profile within fertile window. This indicates that the prediction ability of the three models are moderately good.

As mentioned in Section 2.1, one could use the fitted joint model and hormonal measurement history up to cycle 1, to predict the conditional survival probability $\pi_{i}(j\mid 1)$, for any $j>1$ if that is of interest. The data used in the above analysis are available from the corresponding author upon reasonable request.

We have proposed a joint modeling approach to assess the association between cycle-level geometric features of a woman’s hormonal profile and her fecundity measured by TTP. A likelihood based approach with the use of Gaussian quadrature approximation was proposed for estimation of the unknown parameters. Simulation studies have demonstrated that the proposed estimation approach works reasonable well in the situation similar to the Oxford data, with reasonably small bias and coverage probabilities close to the nominal level. With the estimates from the joint models, we also derived the approach to predict individual characteristics of TTP given a set of longitudinal measurements up to a certain cycle. The prediction of the probability that a woman was subinfertile given her past one menstrual cycle behavior without getting pregnant, was moderately accurate for model with any one of the three geometric features of LH profile, especially with curvature at hormonal profile peak (AUC around 0.65).

The analysis of Oxford data found that couple average age and difference between female and male age are signifcantly associated with TTP in the sense that older couples and couples with larger female minus male age difference have significantly longer TTP. The association between BMI and TTP seem to be marginal with moderate evidence that overweight women have a slower rate of pregnancy while obese women have a faster rate of pregnancy. Multiparous women were found to have significantly shorter TTP than nulliparous women. Furthermore, we found that women with sharper LH peaks, higher LH peaks and overall more curved LH profiles within fertile window tend to have significantly shorter TTP.

We focused on cycle-varying geometric features of the hormonal profile, motivated by biology as well as due to the availability of hormonal data only around ovulation. In the scenario where hormonal profile is available throughout the menstrual cycle, a more general approach of considerable interest is to look at the full hormonal profile and assess patterns in the framework of joint modeling of functional data and time to event. Such approach have been studied by Li et al.(2017,2019,2022)(li2017functional, ; li2019bayesian, ; li2022joint, ) and may be extended in this context of cyclical hormonal profile and TTP to shed light on meaningful patterns on TTP, as well as its predictive ability on infertility. Another interesting approach along this line is the recent development of the functional model proposed for the cyclical longitudinal process (ji2020semiparametric, ), which can be extended in the context of joint modeling.

Finally, our approach though motivated by reproductive epidemiology may be applicable to examples arising in other disciplines. For example, one may be interested in understanding the association between the high blood pressure and risk for heart attack and not just on their average blood pressure measurements. In conclusion, this paper studies a novel model for joint models of longitudinal and survival data where one is interested in the longitudinally varying geometric features with survival data.

Acknowledgments

This research was supported in part by the Intramural Research Program of the National Institutes of Health, Eunice Kennedy Shriver National Institute of Child Health and Human Development. The authors would like to acknowledge that this study utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, MD (http://biowulf.nih.gov). The research of Animikh Biswas was supported in part by the NSF grant DMS 1517027.