In longitudinal biomedical studies, each patient’s process of disease progression is often measured over clinical visits; however, the frequency of and duration between visits can substantially differ between patients and within an individual patient over time. As a result, it becomes very unlikely that any two patients will have the same observed visiting times [19]. Moreover, a given patient might have irregular but dense visit patterns when experiencing a poorer health condition and have regular visits at other periods. On the other hand, patients could forego visits during periods of worsening health. The resulting missingness is not negligible and therefore must be modeled as part of the analysis. In addition, clinical patients might miss part of visits during the observational period due to a health condition, which is not negligible in the analysis. Ignoring the impact of the visiting process can yield biased results in estimating the effect of important clinical factors on the longitudinal measurements in analysis [14].

To overcome this challenge in longitudinal data analysis, two types of approaches have often been applied in the literature. The first type of approach is based on the idea of weighted generalized estimating equations (GEE), which is also widely used in addressing missing data when the longitudinal outcomes are missing not at random [18, 3]. The method relies on inverse-probability weighting (IPW) to estimate the outcomes through GEE models, which provide robust inference at the population level but lack individual-level estimation. A second common approach, which requires additional assumptions to gain individual-level estimation, is based on the shared-parameter joint model framework [21], which refers to jointly modeling the longitudinal measurements and the visiting process by sharing information between two outcomes. The underlying assumption in these joint models is, conditional on the shared information, the longitudinal measurements are independent of the visiting process, and the impact of irregular visits may be estimated by the shared information. Lipsitz and colleagues proposed a joint model linking a linear mixed-effects model for a longitudinal response and another model of the counting process for the number of visits [13]. The model assumed the study observed at random visit time and that the number of visits depends on the previous longitudinal response. Later work by Sun et al. focused on the marginal mean of longitudinal measurements; they proposed a joint model, which used GEE to fit the longitudinal process and a nonhomogeneous Poisson process for the visiting process [28]. Broadly assessing strengths and weaknesses of approaches, Neuhaus et al. (2018) compared a single generalized linear mixed-effects model with random intercept and slope, a GEE with inverse intensity rate ratio and a joint model which includes the generalized linear model and a gamma-distributed visiting process for the cases of the binary and normal responses at designed times of visits [17]. They concluded that mixed-effects models have better performance in estimating the covariate effects. Gasparini et al. accounted for the informative visiting process by fitting a proportional hazard jointly with a linear mixed-effect model using a random intercept [5]. They showed that the joint model had a better fit than the GEE model with normalized inverse intensity weighting.

The aforementioned analysis of continuous longitudinal measurements has shown that mixed-effect models have a smaller bias in scenarios with informative outcomes [17]. Thus, we investigate mixed-effects models with different random effect structures and introduce a Gaussian process (GP) to explain the fluctuation and correlation over time within a subject. A parametric Weibull proportional hazard model is used to fit the visiting process. A log-normal distributed frailty term is estimated to account for the variation in visit behavior between subjects. Previous work proposed a joint model which inserts the exponential frailty term into the longitudinal mixed-effects model to share the component and link the two outcomes [5]. We employ a similar idea, assuming these two outcomes are associated through a hazard ratio (HR) of having a visit at the corresponding time point. To facilitate estimation, we utilize a Bayesian approach. Our main interest is on the longitudinal measurements and covariate:outcome associations in both the mixed-effect and survival models.

This work is motivated by studying the acquisition of lung function measurements at clinical encounters from pediatric patients with cystic fibrosis (CF) lung disease. CF is a life-threatening genetic disease which mainly affects lung function with persistent coughing and frequent lung infections. Based on the most recent reports from CF registries (2016-2018), an estimated 106,000 people live with CF worldwide [23]. In the clinical setting, monitoring lung function via the marker forced expiratory volume in 1 s of % predicted (hereafter, FEV1) is essential to survival in CF, as the primary cause of death is respiratory failure [4]. Attenuated decreases in lung function relative to patient- and/or center-level norms, clinically termed rapid decline, typically manifest during adolescence and early adulthood, but can persist throughout the lifespan [31]. Hence, the US Cystic Fibrosis Foundation guidelines recommend quarterly visits to assess pulmonary function in patients aged 6 years and older [15]; however, patients may complete clinic visits at different frequencies for various reasons. Attendance also varies within an individual patient over time. It’s plausible that increased visit frequency is related to disease severity, but most CF observational studies of lung function aggregate FEV1 data. For example, past studies have aggregated subject-specific FEV1 data across visits into quarterly values prior to modeling [10, 30]. More recent work examining visit frequency as a time-varying covariate in longitudinal FEV1 modeling suggests increased visit intensity is associated with improved lung function [31], while other analyses with conventional mixed-effects modeling have found no evidence for this association [29].

Accurate estimation of CF patients’ lung function can shed light on disease progression and enhance physician efforts to monitor a patient’s pulmonary condition, but the extent to which estimation is confounded by visit frequency has yet to be examined. Figure 1 (a) shows observed FEV1 with highlighted lines marking individual trajectories for 6 representative subjects that were included in our analysis cohort, which is subsequently described. There is substantial between- and within-subject variability in FEV1. The patients also have different lengths of follow-up and distinct visit frequencies during the collection period. Figure 1 (b) is the histogram plot of the gap/waiting time between visits. The gap time is the duration of time (in years) between two visits. The typical (median) between-visit time is around 0.19 years or 2.3 months, shown as the dashed vertical line in Figure 1 (b). The majority of subjects have 0.5 years or 6 months between visits, and the maximum gap time is around 1.5 years. In this paper, we propose a Bayesian joint model for the longitudinal lung function and patients’ visiting processes to determine the extent to which irregular visit patterns may impact lung function decline. We also fit a single mixed-effects model adjusted by the number of visits within the prior year as a time-varying covariate, to account for the impact of irregular visits, which is recommended by Sun and colleagues [28] and has been used in previous CF studies [31]. We compare the joint model with this single mixed-effects model by different Bayesian criteria, and we compare model performance in estimating the effect of clinically relevant factors, such as the coefficients of covariates, after accounting for the visiting process. We also examine the impact of model assumptions through simulation studies.

The paper is organized into the following sections. The notation and model set-ups with different random effects structures are introduced in Section 2. Bayesian priors, model estimation and comparison criteria are described in Section 3. An overview of the CF cohort data and its application, including model comparison, estimation, and inference, are presented in Section 4. Simulation studies are reported in Section 5. We conclude with a discussion of the approach and findings in Section 6. Implementation was performed in R with code available as supplemental files.

We conduct simulation studies in order to examine the performance of the proposed model and the ability of the proposed model to recover the true values of parameters. We investigate four scenarios combining the two factors, the frequency (sparse or dense) of visits and the extent of outcome-dependent or independent visits. The data with dense visiting mimics patients with a higher number of visits but shorter gap times between visits, compared to patients with sparse visits within the same length of time. The data with outcome-dependent visits mimics the outcome-dependent visiting process, compared to outcome-independent visit. For each simulated dataset, we generate the observations for N = 100 subjects. Then we fit the data with the 9 models defined in Section 4.2, which include the single LME models with or without ‘nvisit’ and joint models under different random effects structures. Given the computational intensity, we replicated M = 50 simulated datasets under each scenario. All HMC chains passed convergence diagnostics.

In each simulation, we generate each patient’s visit times using a Weibull distribution, assuming shape parameter of 1.85 and scale parameter of 1 for the scenario of sparse visits (i.e., wider gap time), and scale parameter 0.1 for the scenario of dense visits (i.e., narrower gap time). Average gap time between visits was either 1 or 0.3 years; these annual and quarterly visit scenarios correspond to less frequent vs. guidelines-based visiting processes in the CF example. We simulate two variables, one is a time-varying variable, representing age at visits; initial age ranges from 6 to 20 years old generated from U n i f o r m ( 6 , 20 ); the other is a static variable, representing gender, which is generated from a B e r n o u l l i ( 0.5 ). The true regression parameters of age and gender in the continuous longitudinal response y ( i , j ) are assigned as ( α 1 , α 2 ) ′ = ( − 0.5 , 1 ) ′ and in the survival model as ( β 1 , β 2 ) ′ = ( 0.1 , − 0.2 ) ′ . The intercept α 0 is 85. We generate the random intercept b ( 0 , i ) from the distribution N ( 0 , 1 2 ). The GP terms ψ ( i , j ) are generated from a G P ( 0 , Σ ( 4 , 0.5 ) ). The log of the frailty term u i is generated from a N ( 0 , ( 0.5 ) 2 ). We set the association parameter equal to −2 for the outcome-dependent scenario and equal to 0 for the outcome-independent scenario. We then generate the longitudinal continuous response y ( i , j ) from a normal distribution with mean response function: α 0 + α 1 Age + α 2 Gender + b ( 0 , i ) + ψ ( i , j ) ( t ) + γ exp ( β 1 Age + β 2 Gender + u i ) and the variance σ 2 is set to 1.

After M = 50 replications of each simulation, we summarize the results by the percentage of best performances based on each of the comparison criteria LPML, DIC, WAIC, RMSE, and MAD among all the models compared. A model with a larger value of this percentage means that this model has more chances to perform better in the simulated scenario. To check the ability to recover the true parameters, we used three commonly used criteria, including the coverage probability of 95% credible intervals (CrIs), C 95 , which is calculated as C 95 = ∑ m = 1 M I m ( 95 % CrI contains true parameter value ) M , the bias, which is calculated as Bias = ∑ m = 1 M ( θ ˆ m − θ ) M , with θ ˆ m denoting the posterior mean estimate of the unknown parameter θ at the m th replication, and RMSE, which is calculated as RMSE = ∑ m = 1 M ( θ ˆ m − θ ) 2 M 1 / 2 . A parameter with a larger C 95 , smaller bias, and smaller RMSE indicates a more accurate estimation of the true parameter in the model.

The simulation results for the scenarios of patients with dense visits and longitudinal measurements dependent on the vising process are shown in Table 4, and the independent case is reported in Table 5. Model selection results corresponding to the dependent and independent cases are reported in Tables S3 and S4, respectively. Results for the scenarios of patients with sparse visits are reported in Tables S5-S6.

We first check the models by the order ‘_ n u l l’, ‘_ n v i s i t’, and ‘_ J M’ in all the scenarios. When the responses are dependent on patients’ visits, the ‘_ n u l l’ models across different random effect structures do not provide good estimation for ‘Intercept’ and time-varying covariates, such as ‘Age’, and the estimation for static covariates, such as ‘Gender’ is slightly better but still not satisfactory. When the responses are independent of patients’ visits, the ‘_ n u l l’ models perform better and provide reasonable estimates of covariate effects. The ‘_ n v i s i t’ models are similar to ‘_ n u l l’ estimating ‘Age’ and ‘Gender’ and slightly better in estimating ‘Intercept’, especially in the models only with random intercept. We notice that, in the scenarios of outcome-dependent responses shown in Table 4 and Table S5, although the ‘_ n v i s i t’ model does not recover the true parameters well, the bias and RMSE in the ‘_ n v i s i t’ models are smaller than the ‘_ n u l l’ models. This indicates that, in the analysis of outcome-dependent responses, adding the number of visits within a calendar year into models can help improve the model performance in estimating the effect of covariates and reduce the bias.

Despite aforementioned improvements, many of the estimates are still not satisfactory according to coverage (i.e., low C 95 ) and notably high bias, especially for time-varying covariates. More broadly, the consequences of misspecifying the random effects structure differ markedly between the dependent and independent scenarios. Under the independent scenario, even the clearly misspecified ‘ I n t_’ and ‘ S l o p e_’ models do not deviate substantially from the true model and can sometimes offer comparable performance in terms of bias and RMSE (for example, S l o p e_ n u l l in Tables 5 and S6). In contrast, under the dependent scenario, the impact of misspecification becomes much more severe. The worst observed scenario is the intercept under the simplest model, I n t_ n u l l, in which coverage drops to 0%, compared with 78% under the independent scenario.

The ‘_ J M’ models perform the best in estimating ‘intercept’, ‘Age’ and ‘Gender’ effects, compared to the ‘_ n u l l’ and ‘_ n v i s i t’ models within each random effect structure in all four simulation scenarios. This result is especially enhanced in the scenarios when the longitudinal responses depend on the visiting process. For the association parameter γ = − 2, ‘_ J M’ models provide accurate estimates, recovering the true value with small bias, which implies the joint models can detect the association between the longitudinal measurements and patient visits.

Comparison between the dependent and independent scenarios further highlights the advantages of the ‘_ J M’ models. Under the independent scenario, where the data generating process (DGP) follows the G P_ n u l l model, the G P_ n u l l model and the more complex G P_ n v i s i t and G P_ J M models yield similar inference in terms of bias and RMSE. Among the three GP-based models, ‘ G P_ J M’ is most frequently identified as the best performing model (Table S4). Under the dependent scenario, where DGP is based on the ‘ G P_ J M’ model, any model that excludes the joint modeling component results in zero coverage for the coefficient of the time-dependent covariate and poor coverage for the other two coefficients. Overall, while the ‘_ J M’ models perform well in both dependent and independent scenarios, model performance deteriorates rapidly when dependence is ignored. Although the ‘_ J M’ model is more complex, its additional flexibility allows it to better accommodate key features of the data.

This pattern also holds for the GP components of the random effects. Although our study does not include data replicates generated under a DGP without the GP component, Su et al. demonstrated that the GP model performs well for independent and identically distributed (iid) random effects, whereas fitting data that contain GP-structured random effects with an iid-only model leads to substantial loss of accuracy [27]. The consequences of such misspecification are similarly evident in our comparison of ‘ I n t_ n u l l’, ‘ S l o p e_ n u l l’ and ‘ G P_ n u l l’ in Table 5.

By comparing the models across the different frequencies of visits, we notice that both ‘_ n u l l’ and ‘_ n v i s i t’ models perform better in the sparse visits than the dense visits. Meanwhile, ‘_ J M’ models perform slightly better in the dense visits than in the sparse visits. By comparing across the random effect structures, when the longitudinal responses are dependent on the visits, we found a big improvement in parameter estimation of ‘Intercept’ and ‘Gender’ by including a random slope or GP term effects. The joint models provide a better estimation of ‘Age’ across the random effect structures by checking C 95 , bias and RMSE, but it is worth noting that running a single simulation replicate could take up to four hours for the most complex setting. This approximate run time is consistent with the computational intensity that we have observed in prior simulation studies of similar model structures [26]. When the responses are independent of visits, the improvement in estimating the ‘Intercept’ and ‘Gender’ is reduced. This evidence implies that, when there is an outcome-dependent response, the models with random slope or the GP term partially account for the visiting effect and improves the model estimation of the intercept and static covariates. Meanwhile, by estimating the visiting process, the proposed joint models appear to take care of the effect of dependency on the outcome.

In this paper, we have proposed a joint model for the continuous longitudinal measurement in the presence of an outcome-dependent visiting process. We examined the model performance among three different random effects structures with an application to CF data and simulation studies. Based on our simulation studies, we conclude that the joint models with a visiting process provide more accurate parameter inference in longitudinal analysis, compared to a single model structure. Furthermore, this is evident no matter how the longitudinal response pattern varies with respect to patients’ visiting process. In the studies with dense visits, this proposed joint model is preferred over the models with random intercept or random slope. We found that models adjusted by patient’s visiting behavior, such as the number of visits within the prior year, can help improve the model performance and reduce the bias in parameter estimation when the longitudinal measurements are dependent on patient’s visits; however, the coefficient estimate of a given time-varying covariate needs to be carefully examined in these settings. The random slope structure and a GP term can explain some part of patients’ visiting effect and improve the overall model fitting, but may not aid parameter estimation. We also observed that an improper assumption in the random effect structure impairs the model performance.

In the application to the CF study, the random intercept models do not work well compared with other models. Furthermore, the models with random intercept and GP over time have the best performance of those considered. This result is likely due to the large heterogeneity in CF data [30, 32, 1], which can be estimated by the nonlinear correlated structure from the GP term. The improvement from a single longitudinal model to the proposed joint model is large when only using a random intercept in modeling the CF FEV1 data. However, when adding the random slope or a GP term, though the joint model is the best fitting in terms of LPML, this improvement is reduced. From model estimation, the lung function measurement FEV1 and the patient’s visiting process have a negative association with estimation −1.99 (−3.47, −0.67), which shows that irregular visiting does affect the FEV1 trajectory estimation. This indicates patients who have lower FEV1 are likely to have more frequent visits. This result could be because patients at a certain period need more frequent visits due to a sudden lung function decline or more severe patients are recommended with more visits than usual. The average gap time between visits in this CF data is 0.18 years (2.16 months) with a range from 0.1 to 1.44 years; from the estimation of the joint model, the median gap time is estimated as 0.19 years (2.3 months).

Covariates were selected for each submodel in our real-data application based on information from prior CF studies. Our proposed joint modeling framework allows the same covariates to be included in both the longitudinal and event submodels, but additional care is needed when interpretating their effects. In the longitudinal submodel, a given covariate serves as a predictor of the lung function trajectory, describing systematic differences in the mean level and evolution of lung function. In the hazard submodel, the same covariate can exert a direct effect on the instantaneous risk of an observed visit that is not mediated by the longitudinal lung function. In our specification, the association structure is introduced through the longitudinal submodel, which captures how covariates affect lung function indirectly through their impact on the HR for having an observed visit. Consequently, when a covariate such as sex is included in both submodels, its overall effect on lung function decomposes into a direct effect (through the fixed effects in the longitudinal submodel) and an indirect effect mediated by the visiting process; that is, through f ( u i , β ∣ W ( i , j ) ). The association parameter γ summarizes the strength and direction of this indirect effect of the visiting process on lung function, conditional on the covariates and random effects. Adaptations to allow causal inference under the proposed model would require a different specification of the functional form and shared parameter to mimic the “reverse” shared parameter setup assumed in our proposed joint model [33].

Despite the insights gained on modeling properties and the application, our approach has some limitations. In the CF application, observations may be left truncated due to calendar time, since our visit records began in 2012. In larger analyses with CF registries, left truncation has been accounted for by using age at diagnosis or birth cohort [30, 32]; however, we did not consider these covariates, given the sample size of our center-level cohort. In addition, observations are truncated at the time of death or lung transplant, and we assume these dropouts are not informative. For pediatric applications, however, there are few deaths or lung transplants during visits (Table 1). Informative dropout would need to be considered if the CF analysis cohort includes older patients. Acute drops in a patient’s lung function may precede these events; thus, it might be helpful to consider them in estimating lung disease progression.

Extensions of our current work with application to a larger data source could be considered to account for the aforementioned types of censoring. For example, building the model as a multivariate joint model could account for a terminal event, such as death. We assumed in the CF application that population-level lung function declines linearly and that the GP captures nonlinear progression within an individual patient. This assumption is appropriate in the CF context, since our analysis cohort is pediatric. However, in longitudinal analysis, long-term sequences of visits, such as those over a life span, population-level linearity may not be a valid assumption. Potential extensions of our application are to use a nonlinear form for the mean structure, such as a semiparametric form in the continuous longitudinal model [30], incorporate non-Gaussian terms [1], or characterize survivorship via a multivariate joint model [22]. Furthermore, we assume the lung function measurements are associated with the hazard of having a visit, but different forms of association structures could be investigated in future applications. Those may include additional random effects to account for clustering in multi-center cohort studies [20]. Additional research is needed to extend the proposed model to settings with non-proportional hazards. Currently, we use the standard parametric Weibull distribution to model the gap times. If other parametric survival models are more suitable for the data, or if greater flexibility is needed, it would be worthwhile to explore alternative parametric distributions or even nonparametric approaches for modeling the gap times. Finally, another extension of our current work is jointly modeling two dependent outcomes via a latent process, wherein a general form of a joint model is linked through a latent stochastic process; this has been discussed by Henderson and colleagues [9].

The authors do not have permission to share the clinical data, but all code and the simulated data set have been provided as supplemental material for reproducibility purposes. For those interested in acquiring the clinical data, please contact the corresponding author for additional information.