Recurrent events modeling based on a reflected Brownian motion with application to hypoglycemia

Abstract

Patients with type 2 diabetes need to closely monitor blood sugar levels as their routine diabetes self-management. Although many treatment agents aim to tightly control blood sugar, hypoglycemia often stands as an adverse event. In practice, patients can observe hypoglycemic events more easily than hyperglycemic events due to the perception of neurogenic symptoms. We propose to model each patient’s observed hypoglycemic event as a lower boundary crossing event for a reflected Brownian motion with an upper reflection barrier. The lower boundary is set by clinical standards. To capture patient heterogeneity and within-patient dependence, covariates and a patient level frailty are incorporated into the volatility and the upper reflection barrier. This framework provides quantification for the underlying glucose level variability, patients heterogeneity, and risk factors’ impact on glucose. We make inferences based on a Bayesian framework using Markov chain Monte Carlo. Two model comparison criteria, the deviance information criterion and the logarithm of the pseudo-marginal likelihood, are used for model selection. The methodology is validated in simulation studies. In analyzing a dataset from the diabetic patients in the DURABLE trial, our model provides adequate fit, generates data similar to the observed data, and offers insights that could be missed by other models.

1. INTRODUCTION

Diabetes is a group of diseases characterized by elevated blood glucose. It is a major cause of kidney failure, non-traumatic lower-limb amputations, blindness, heart disease, and stroke. As a result, diabetes is one of the leading causes of death. Diabetes affects 37.3 million Americans which is 11.3% of the US population (Centers for Disease Control and Prevention 2022). The goal of treating diabetes patients is to lower their blood glucose in range. When blood glucose level is too high, patients will experience hyperglycemia. Lowering glucose too much, however, could result in an adverse event called hypoglycemia. Concern over hypoglycemia has a significant negative impact on diabetes management, making it a major factor to prevent patients’ glycemic control from reaching treatment target (Wild et al. 2007). It is desired to develop new anti-diabetes agents that lead to less hypoglycemic events while lowering the glucose toward the normal range.

Hypoglycemias are relatively easier observed than hyperglycemias in practice. Symptoms of hypoglycemia are the result of brain glucose deprivation or the perception of physiological change. Awareness of hypoglycemia is mainly the result of the perception of neurogenic symptoms (Towler et al. 1993; DeRosa and Cryer 2004). Times of hypoglycemic events are recorded through patients’ self-report diaries with easily observed symptoms, such as headache, shaking, sweating, hunger, and fast heartbeat (Cryer et al. 2003). In contrast, hyperglycemia symptoms are much less obvious, and may not be noticed unless a blood sugar level test is performed (Cryer et al. 2009). Therefore, only the times of recurrent hypoglycemic events are reliably available in self-reported data. These event times have been the target of modeling in diabetes clinical research (e.g. Fu et al. 2016; Ma et al. 2021; Doubleday et al. 2022).

A variety of recurrent event models have been available. Earlier works characterize the intensity function of the recurrent event process (Andersen and Gill 1982). The multiple events can be analyzed with the conditional intensity for the sequence of events given the history (Prentice et al. 1981; Lee et al. 1992) or with marginal intensities of the gaps between successive events (Wei et al. 1989; Huang and Chen 2003; Schaubel and Cai 2004). Alternative approaches to gap times have been proposed, such as the accelerated failure time model (Chang 2004), the additive hazards model (Sun et al. 2006), and the quantile regression (Luo et al. 2013). Heterogeneities beyond covariates can be incorporated into the models as a subject-level random effect, which also captures the dependency between events within an individual (Klein 1992; Duchateau et al. 2003; Box-Steffensmeier and De Boef 2006). Less restrictive model assumptions have been considered by marginal rate or marginal mean models (Lawless and Nadeau 1995; Lin et al. 2000), which do not fully specify the recurrent event process. More recently, a general scale-change model encompasses multiple existing event time models in a unified framework and allows a flexible form of informative censoring (Xu et al. 2020). The readers are referred to recent reviews for details (Cook and Lawless 2007; Charles-Nelson et al. 2019).

An alternative approach is the first hitting time (FHT) model. An FHT model has 3 components: a stochastic process, a boundary, and a time scale on which the process unfolds (Lee, 2019). An event occurs when the underlying stochastic process hits the preset boundary. It is well-known that the distribution of the FHT of a Brownian motion is the inverse Gaussian distribution (Schrödinger 1915; Folks and Chhikara 1978). Incorporating covariates into the model parameters leads to threshold regression (Lee and Whitmore 2006). Random effects can be incorporated to account for unmeasured covariates (Pennell et al. 2010). For recurrent events, Whitmore et al. (2012) modeled the event times as a sequence of independent and identically distributed hitting times of a Wiener process as it passes through successive equally-spaced levels. Economou et al. (2015) extended the inverse Gaussian threshold regression model with a random effect. Malefaki et al. (2015) further extended the model to allow censoring to occur at every intermediate stage during the recurrent event process. For recurrent hypoglycemic events of diabetic patients, it is natural to model them as hitting times when the glucose level crosses a lower barrier but the unobserved hyperglycemic event times need to be handled.

We propose to model hypoglycemic event times as the FHTs of a Brownian motion with an upper reflection boundary hitting a lower barrier. The upper reflection barrier represents that patients’ glucose levels will eventually decrease, potentially due to the use of glucose-lowering medications that are routinely prescribed for this patient population. Consequently, hyperglycemic events at which the patient’s glucose levels hit the upper reflection barrier are typically not observed. Conversely, the lower barrier represents a boundary beyond which patients experience uncomfortable symptoms associated with hypoglycemia. As a result, hitting times for the lower barrier are typically observed. The distribution of the FHT of the lower barrier of a Brownian motion with an upper reflection barrier has been studied by Hu et al. (2012). Taking advantage of this distribution, we model the recurrent hypoglycemic events based on a sequence of reflected Brownian motions hitting a lower barrier, after which the process restarts at a preset point between the 2 barriers. The gap times between the successive events of the same patient are assumed to be independent conditional on a subject-level random effect. Covariates and a subject-level random effect are linked to parameters of the FHT model. The model provides unique opportunities to characterize the variability and heterogeneity of hypoglycemic events with an intuitive interpretation.

Replacing the Brownian motion with a reflected Brownian motion in an FHT model introduces a significant level of complexity. The density and the distribution functions of the FHT derived in Hu et al. (2012) have not been used in the statistical literature. Both functions involve infinite series, which present challenges to accurate evaluation of the loglikelihood and efficient design of random number generation. To address the challenges, we first demonstrated that the right tail of the FHT is bounded by an exponential rate. Drawing on this result, we implemented the density and distribution functions and formulated an efficient rejection sampling algorithm with a 3-piece proposal density. Inferences are conducted within a Bayesian framework using Markov chain Monte Carlo (MCMC) with our implementation of the FHT distribution incorporated into the generic algorithm in R package NIMBLE (de Valpine et al. 2017). Different models, especially those with different types of random effects, are compared with 2 Bayesian model comparison criteria, deviance information criterion (DIC) (Spiegelhalter et al. 2002) and logarithm of the pseudo-marginal likelihood (LPML) (Geisser and Eddy 1979; Gelfand and Dey 1994). The whole methodology was validated in numerical studies before being applied to analyze the hypoglycemic event data in our motivating application. The R code is publicly available in https://github.com/YingfaX/reflbrown.

The remainder of the article is organized as follows. In Section 2, an FHT model of a reflected Brownian motion is set up and an efficient rejection sampling algorithm is proposed for simulation from the model. Bayesian parameter estimation and model selection methods are presented in Section 3. Simulation studies that validate the Bayesian inferences are reported in Section 4. The methodologies are applied to the hypoglycemic event times from the DURABLE trial (ClinicalTrials.gov Identifier: NCT00279201) in Section 5. Section 6 concludes with a discussion.

2. MODEL

The fundamental component of our model is the FHT model for the first time when a Brownian motion with a reflecting upper barrier hits a lower barrier. A sequence of such models is assembled for recurrent events. Subsequently, we incorporate covariates and frailty terms in the model parameters to further characterize variability and heterogeneity.

2.1. FHT of a reflected Brownian motion

Consider a no-drift Brownian motion |$X(t)$| with volatility |$\sigma$|⁠. Let |$\kappa$| be the upper reflection barrier and |$\nu$| be the lower hitting barrier, where |$\kappa \gt\nu$|⁠. Suppose that |$X(0)=x_{0}\in(\nu,\kappa]$| is the starting point. The first time when |$X(t)$| hits |$\nu$| is |$\tau:=\inf\{t\,\gt\,0; X(t)=\nu\}$|⁠. For any |$\nu\in[0, x_{0})$|⁠, the density and distribution function of |$\tau$| are, respectively (Hu et al. 2012),

where for |$n\,=\,1,2,\ldots$|⁠,

Note that |$0 \lt\lambda_{1} \lt\lambda_{2} \lt\cdots$|⁠, |$\lambda_{n}\to\infty$|⁠, and |$\sum_{n\,=\,1}^{\infty}c_{n}=1$|⁠. It is tempting to think of |$f(t)$| as a mixture of exponential densities, but this is not true because |$c_{n}$|’s can be negative.

Random number generation from density (2.1) has not been investigated in the literature, we propose a rejection sampling algorithm. By construction, the proposal density |$g(t)$| is required such that |$f(t)/g(t)$| is bounded over the support |$t\,\gt\,0$|⁠. It is known that |$f(t)/t\to 0\, \text{as}\, t\downarrow 0$| (Hu et al. 2012: 13). In  Appendix A, we show that the right tail of |$f(t)$| is bounded by that of an exponential density with rate |$\lambda_{1}$|⁠:

Given these properties, we propose a 3-piece proposal density which consists a left tail, a body, and a right tail. The left tail is a triangular density that connects the origin to the peak of the target density peak |$(t_{m},f(t_{m}))$|⁠, where |$t_{m}$| is the mode of the target density |$f$|⁠. The body is a trapezoid density that connects the peak |$(t_{m},f(t_{m}))$| to a user-defined |$q^{th}$| quantile point |$(t_{q},f(t_{q}))$|⁠, where |$t_{q}$| is the upper |$q$|th quantile of |$f$|⁠, and |$q$|⁠, for example, can be set to be 0.95. The right tail beyond |$t_{q}$| is the tail of the exponential density with rate |$\lambda_{1}$|⁠. The details of the implementation of the algorithm are in Appendix B.

2.2. Modeling gap times of hypoglycemic events

We model the recurrent hypoglycemic events based on a sequence of reflected Brownian motions hitting a lower barrier |$\nu$|⁠. The Brownian motion has volatility |$\sigma$| and an unobserved reflection upper barrier |$\kappa$|⁠, beyond which hyperglycemia occurs but may not be recorded. The reflected Brownian motion restarts at a known starting point |$x_{0}$| after hitting the lower barrier |$\nu$|⁠, where |$x_{0}$| can be set as the average level from the patients group. This is motivated by that patients’ glucose level arises back, possibly after eating to relieve quickly the symptoms of hypoglycemia. The gap times between 2 successive hypoglycemic events thus follow an FHT distribution. Figure 1 shows a sample path for a subject who experiences 2 events at time |$T_{1}=19$| and |$T_{2}=179$|⁠, and finally is censored at time 220.

To incorporate heterogeneity, we allow volatility |$\sigma$| and upper barrier |$\kappa$| to depend on subject-level covariates. Appropriate link functions are necessary to ensure that |$\sigma \gt 0$| and |$\kappa \gt x_{0}$|⁠. Subject-level random effects can be further added to the regression. In particular, the models for |$\sigma_{i}$| and |$\kappa_{i}$| with the log link are

respectively, where |$\boldsymbol{X}_{i}$| is a |$p$|-dimensional covariate vector for subject |$i$|⁠, |$\boldsymbol{\beta}$| and |$\boldsymbol{\alpha}$| are |$p$|-dimensional regression coefficient vectors, and |$(Z_{i1},Z_{i2})$| is a bivariate normal random effect with mean |${\bf 0}$|⁠, variance |$\boldsymbol{\theta}=(\theta_{1},\theta_{2})$|⁠, and correlation |$\rho$|⁠. Also note that |$\boldsymbol{X}_{i}$| is considered a time-independent covariate vector.

For ease of implementation, the bivariate normal vector |$(Z_{i1},Z_{i2})$| can be reparametrized as |$(Z_{i1},\gamma Z_{i1}+Z_{i2}^{\prime})$|⁠, where |$Z_{i1}$| and |$Z_{i2}^{\prime}$| are independent normal variables with mean zero and variance |$\theta_{1}$| and |$\theta_{2}^{\prime}$|⁠, respectively. Then the variance of |$Z_{i2}$|⁠, |$\theta_{2}$|⁠, is reparametrized as |$\gamma^{2}\theta_{1}+\theta_{2}^{\prime}$|⁠, and the correlation |$\rho$| between |$Z_{i1}$| and |$Z_{i2}$| is reparametrized as |$\gamma\sqrt{\theta_{1}/(\gamma^{2}\theta_{1}+\theta_{2}^{\prime})}$|⁠. The reparametrized models for |$\sigma_{i}$| and |$\kappa_{i}$| are, respectively,

The starting point |$x_{0}$| and lower barrier |$\nu$| are set to be fixed depending on the real data and clinical trial standard, which are shared for all subjects. In general, the subject with higher volatility and smaller reflected barrier is associated with higher risk of hypoglycemic event.

3. BAYESIAN INFERENCE

Suppose that the event times of hypoglycemia are observed for a group of |$n$| patients during their follow-up times. For |$i\,=\,1,\ldots, n$|⁠, let |$t_{ij}$|⁠, |$j\,=\,1,\cdots, n_{i}$| be the |$j$|th of the |$n_{i}$| observed gap times of subject |$i$|⁠, with the final gap time being censored. To ensure notational generality, define |$\delta_{ij}=1$| if |$t_{ij}$| represents an event and |$\delta_{ij}=0$| if it is censored. In addition to a |$p$|-dimensional covariate vector |$\boldsymbol{X}_{i}$|⁠, the observed data include |$\boldsymbol{t}_{i}=(t_{i1},\cdots, t_{in_{i}})^{\top}$| and |$\boldsymbol{\delta}_{i}=(\delta_{i1},\cdots,\delta_{in_{i}})^{\top}$|⁠, |$i\,=\,1,\ldots, n$|⁠. The parameters to be estimated are |$\mathbf{\Omega}=(\boldsymbol{\alpha}^{\top},\boldsymbol{\beta}^{\top},\theta_{1},\theta_{2}^{\prime},\gamma)^{\top}$|⁠.

3.1. Likelihood, prior, and posterior

Given the covariate vector |$\boldsymbol{X}_{i}$| and subject-level frailties |$\boldsymbol{z}_{i}=(z_{i1},z_{i2}^{\prime})$|⁠, the gap times of subject |$i$| are assumed to be conditionally independent. The conditional likelihood contribution of subject |$i$| given the unobserved frailties is

where the density and distribution functions are given by (2.1) and (2.2), respectively, and |$\sigma_{i}$| and |$\kappa_{i}$| are defined in (2.3).

In our case, the event times were recorded in days. As the result, |$t_{ij}$| can be treated as interval-censored. Consequently, the contribution to the likelihood of subject |$i$| is then rewritten as

where

Prior distributions of the parameters need to be specified to complete the Bayesian model. For regression coefficients |$\boldsymbol{\alpha}$| and |$\boldsymbol{\beta}$|⁠, normal priors with zero mean and a large variance are specified. Inverse gamma priors are specified for the normal variances of the frailties |$\theta_{1}$| and |$\theta_{2}^{\prime}$|⁠. For the reparametrized coefficient of the frailty |$\gamma$|⁠, a vague normal prior with mean zero is imposed. The priors are summarized as follows:

where |$N(0,\nu^{2})$| is the normal distribution with mean zero and variance |$\nu^{2} \gt 0$|⁠, and |${\rm IG}(a, b)$| is the inverse gamma distribution with shape |$a\,\gt\,0$| and scale |$b\,\gt\,0$|⁠. In our numerical studies, the hyper-parameters were set to be |$\sigma^{2}_{\alpha}=\sigma^{2}_{\beta}=\sigma^{2}_{\gamma}=10^{2}$|⁠, and |$a\,=\,b\,=\,1$|⁠.

With |$\boldsymbol{z}=(\boldsymbol{z}_{1},\cdots,\boldsymbol{z}_{n})^{\top}$|⁠, combining the likelihood function and prior distributions leads to the joint posterior density

where |$q(\cdot)$| denotes a generic density function of its argument, |$L_{i}(\mathbf{\Omega}|\boldsymbol{t}_{i},\boldsymbol{\delta}_{i},\boldsymbol{X}_{i},\boldsymbol{z}_{i})$| is the conditional likelihood function for subject |$i$| given in (3.4), |$g(\boldsymbol{z}_{i}|\theta_{1},\theta_{2}^{\prime})$| is the independent bivariate normal density of the reparametrized frailties of subject |$i$|⁠, and all the |$q(\cdot)$|’s are priors given in (3.5). Since all the priors are proper, the posterior is proper.

To make inferences about the parameters, we use MCMC. The full conditional distributions of all the parameters are sampled with the Metropolis–Hasting algorithm because the FHT density is not of any existing standard form. We incorporated our customized FHT distribution into a generic implementation through R package NIMBLE (de Valpine et al. 2017). A good number of draws can be obtained after thinning those long chains, which also helps reduce the computational burden in the calculation of model comparison criteria discussed in a later section.

3.2. Model selection

With Dev|$(\mathbf{\Omega})$| denoting the deviance calculated by negated observed-data loglikelihood evaluated at |$\mathbf{\Omega}$|⁠, DIC is given as

where |$\bar{\mathbf{\Omega}}$| is the posterior mean of |$\mathbf{\Omega}$|⁠, |$p_{D}=\overline{{\rm Dev}}(\mathbf{\Omega})-{\rm Dev}(\bar{\mathbf{\Omega}})$| is the effective number of parameters, and |$\overline{{\rm Dev}}(\mathbf{\Omega})$| is the posterior mean of |${\rm Dev}(\mathbf{\Omega})$|⁠. This calculation based on the observed-data likelihood is the |${\rm DIC}_{3}$| in Celeux et al. (2006). Since there is no closed form of the observed-data likelihood, we used Monte Carlo integration to approximate it. The observed-data likelihood of subject |$i$| has the form

where |$L_{i}(\mathbf{\Omega}|\boldsymbol{t}_{i},\boldsymbol{\delta}_{i},\boldsymbol{X}_{i},\boldsymbol{z}_{i})$| is given in (3.4). The Monte Carlo approximation for the observed-data likelihood for subject |$i$| is

where |$\boldsymbol{z}_{i}^{(1)},\cdots,\boldsymbol{z}_{i}^{(M)}$| are Monte Carlo samples that each can be easily generated by the independent normal distributions, and |$M$| is the Monte Carlo sample size. The deviance is given by

Models with smaller DIC values are preferred.

The other criterion LPML is calculated based on conditional predictive ordinate (CPO). Let |$D_{-i}=\{(\boldsymbol{t}_{j},\boldsymbol{\delta}_{j},X_{j}):j\,=\,1,\ldots, n; j\neq i\}$|⁠, denote the observed data with the |$i$|th subject excluded. The CPO for the |$i$|th subject is the leave-one-out predictive likelihood

where |$q(\mathbf{\Omega}|D_{-i})$| is the marginal posterior distribution of |$\mathbf{\Omega}$| with the |$i$|th subject excluded. The Monte Carlo estimate of |${\rm CPO}_{i}$| (Dey et al. 1997) is

where |$\mathbf{\Omega}^{(k)}=(\boldsymbol{\alpha}^{(k)},\boldsymbol{\beta}^{(k)},\theta_{1}^{(k)},\theta_{2}^{^\prime(k)},\gamma^{(k)})^{\top}$|⁠, |$k\,=\,1,\cdots, K$|⁠, is the |$k$|th iteration of posterior draws. Each term |$L_{i}$| in Equation~(3.7) can be approximated the same way as in (3.6). Then the LPML can be calculated by

Models with higher LPML values are preferred.

4. SIMULATION

Simulation studies were conducted to evaluate the performance of the Bayesian estimator and the model selection criteria. To mimic the real application analyzed in the next section, we considered a setting of 2 covariates for each subject |$i$|⁠: |$X_{i1}$|⁠, baseline insulin level and |$X_{i2}$|⁠, baseline body mass index (BMI). A bivariate gamma distribution with a normal copula was fitted to the real data with R package copula (Hofert et al. 2018). The covariates were then generated from the fitted bivariate gamma distribution.

Three frailty models were used to generate recurrent events with covariate vector |$\boldsymbol{X}_{i}=(1, X_{i1},X_{i2})^{\top}$|⁠. Model 1 is the full model (2.3) with |$\boldsymbol{\alpha}=(2.9,0.2,-0.1)^{\top}$|⁠, |$\boldsymbol{\beta}=(0.9,-0.2,-0.1)^{\top}$|⁠, and |$\gamma=-0.55$|⁠; the frailties |$z_{i1}$| and |$z_{i2}^{\prime}$| were generated from the normal distribution with mean zero and variance |$\theta_{1}=0.2$| and |$\theta_{2}^{\prime}=0.3$|⁠, respectively. The starting point |$x_{0}$| and the lower barrier |$\nu$| were set to be |$x_{0}=10$| and |$\nu\,=\,3.9$|⁠, respectively, which will be discussed in Section 5. Model 2 is the same as Model 1 except that |$\gamma\,=\,0$|⁠, which means that the frailties in the upper reflection barrier and the volatility are independent. Model 3 is also a reduced model of Model 1 with the restrictions |$\gamma=-1$| and |$\theta_{2}^{\prime}=0$|⁠, which implies that the upper reflection barrier and the volatility share the same frailty. The negative sign of gamma was chosen as suggested by the real data analysis. For ease of reference, the 3 models are referred to as, respectively, correlated-frailty (CF), independent-frailty (IF), and shared-frailty (SF).

For each model, the gap times were generated with the rejection sampling algorithm in Section2.1 until the sum of generated gap times is equal to or larger than the follow-up time of each subject |$i$|⁠, where the follow-up times were independently generated from the empirical distribution of the follow-up times of the real data. The generated gap times were then rounded into integers as days. Several sample size levels, |$n\in\{200,400,800,1,600\}$| were considered. For each simulation setting, 200 datasets were generated.

The Bayesian analysis of each dataset used prior distributions specified in Equation (3.5). For the coefficient parameters, we set |$\sigma^{2}_{\alpha}=\sigma^{2}_{\beta}=\sigma^{2}_{\gamma}=10^{2}$|⁠, which leads to vague but proper prior distributions. The hyper-parameters |$(a, b)$| of the inverse gamma prior for the frailty variances are often set to be |$a\,=\,b\,=\,1$| (Manda and Meyer 2005; Gelman 2006). Given a dataset, an MCMC was run with 90,000 iterations. The first 15,000 iterations were discarded as burn-in period and the rest were thinned by 10. The convergence of MCMC was checked with Heidelberger and Welch’s diagnostic (Heidelberger and Welch 1983), which is available in R package CODA (Plummer et al. 2006). The posterior inferences were obtained from the resulting 7,500 MCMC samples. However, the calculation of model comparison criteria, which includes integrating random effects with Monte Carlo approximation, is computational intensive. To reduce the computational cost, we further thinned the remaining posterior samples by 15, resulting in 500 posterior samples for the calculation of the model comparison criteria.

Table 1 summarizes the results of the Bayesian estimation for all 3 models when fitted with correct specifications. The posterior mean of each parameter was considered as the point estimator. The empirical bias of the estimates for all parameters is close to zero except that of |$\theta_{2}^{\prime}$| when |$n\,=\,200$| for the CF and the IF model; in both cases the bias reduces as |$n$| increases. The mean of the posterior standard deviation of the estimates agrees closely with the empirical standard deviation of the point estimates for most parameters even for sample size |$n\,=\,200$|⁠. Consequently, the empirical coverage rates of the |$95\%$| highest posterior density (HPD) credible intervals are close to the nominal level, and the agreement improves in general as |$n$| increases from 200 to 400. Results for larger sample sizes (800 and 1,600) are provided in Section S1.

To evaluate the performance of model comparison criteria, we fitted 3 models to each dataset. DIC and LPML were used to select among the 3 fitted models. The Monte Carlo sample size used for approximating the integrals was set to be |$M\,=\,500$|⁠. Table 2 presents the selection frequencies for candidate models under each scenario, based on either DIC or LPML, for sample size |$n\in\{200,400,800,1,600\}$|⁠. On average, the correctly specified models have the smallest DIC and highest LPML. The results suggest good performance of both DIC and LPML in selecting the right models. For example, when the CF model is the data-generating model, the proportion of correctly identifying the true model increases for both DIC and LPML as the sample size increases, reaching |$95\%$| and |$90\%$|⁠, respectively, with DIC slightly outperforming LPML. When the 2 reduced models are the data generating models, both DIC and LPML still select them with the highest frequency, but the tendency to choose the full model also increases when the sample size increases. Such observations echo the limitations of DIC and LPML in distinguishing the true model and an overfitted model (Maity et al. 2021). In our application, both criteria select either the correct model or the full model, which still provides valuable information for practitioners.

5. HYPOGLYCEMIC EVENT TIME ANALYSIS

The proposed models were applied to analyze the hypoglycemic event times from the DURABLE trial (Buse et al. 2009). Between 2005 and 2007, 2,187 patients with type 2 diabetes from 11 countries were enrolled in the study. The dataset contains the possibly censored times of hypoglycemic events of the patients during their follow-up periods. Also, available are a collection of baseline covariates, which allows assessments of risk factors for hypoglycemia among the patients. The median follow-up time of the patients is 168\,days. Continuous baseline covariates include fasting blood glucose, fasting insulin, adiponectin, weight, height, BMI, systolic blood pressure, diastolic blood pressure, heart rate, and duration of diabetes. Summaries of the continuous covariates are presented in Table 3. Three important variables, fasting glucose level, adiponectin level, and fasting insulin level, have extremely high values, which call for prudence in data analysis. Two categorical variables are available. The first one is starter insulin regimens with 2 levels, twice-daily lispro mix (LM) |$75/25$|⁠, |$75\%$| lispro protamine suspension and |$25\%$| lispro (referred to as LM 75/25 hereafter), versus once-daily insulin glargine. The second one is the usage of oral antihyperglycemic drugs with 3 levels: thiazolidinedione, sulfonylurea, and both. All the available covariates are subject-level and time-independent.

After excluding the subjects with missingness in covariates or outside reference range, the dataset contains |$n\,=\,1,943$| patients. Prior to model fitting, all the continuous covariates were standardized. Log transformation was applied to 2 right-skewed covariates, baseline adiponectin and baseline fasting insulin, before standardization. Among the 1,943 patients, 570 (29%) received both oral antihyperglycemic drugs, 1,207 (62%) only received sulfonylurea, and 166 (9%) only received thiazolidinedione. For ease of discussion, the group that received both drugs was used as the reference group; 2 dummy variables, |$\mathsf{sulf{-}Only}$|⁠, which is 1 if only received sulfonylurea, and |$\mathsf{tzd{-}only}$|⁠, which is 1 if only received thiazolidinedione, were included. Define an indicator variable for the insulin regime |$\mathsf{LM}$|⁠, which equals 1 for the 959 (49%) patients who received LM |$75/25$| and 0 for the 984 (51%) patients who received glargine. Some patients had multiple hypoglycemic events within a single calendar date. In this case, the gap times between successive hypoglycemic events were recorded as zero. This was handled by treating the gap times in days as interval-censored, using the likelihood constructed with (3.4) in Section 3.1. The daily hypoglycemic event rates of the patients have a wide range from 0 to 0.77 with mean 0.07. These descriptive statistics indicate the existence of severe heterogeneity risk of hypoglycemia among subjects.

The 3 models along with their priors investigated in Section 4 were fitted to the DURABLE data. The lower boundary was set to 3.9\,mmol/l (70\,mg/dl), which is the clinical standard for hypoglycemic events (Seaquist et al. 2013). The starting point |$x_{0}$| of the Brownian motion after each hypoglycemic event was set to 10, which is the rounded integer of the median baseline fasting glucose level of all patients. Sensitivity analyses of alternative choices for the starting point |$x_{0}$| and priors were conducted and will be discussed later in this section. For each model, an MCMC was run for 110,000 iterations and thinned by 10 after discarding the first 10,000 iterations as burn-in. The choice of burn-in and convergence of the MCMC chains were monitored by traceplots that are provided in Section S2. The resulting 10,000 posterior samples were used for inference, and we further thinned the posterior samples for model comparison criteria calculation to reduce computational cost. The results of DIC and LPML for the 3 models are presented in Table 4. Both criteria suggest that the CF model and the IF model are similar, both of which are preferred to the SF model. Given that the CF and the IF have close model fit, we chose the IF model as it is more parsimonious. That is, the 2 frailties in the upper reflection barrier and the volatility can be treated as independent.

Table 5 summarizes the estimated model parameters, their standard errors, and 95% HPD credible confidence intervals from the fitted IF model. The results from the volatility model suggest that patients with higher baseline fasting blood glucose level, lower adiponectin, higher fasting insulin, higher height, higher BMI, higher diastolic blood pressure, and lower duration of diabetes are significantly associated with lower volatility and, hence, lower risk of hypoglycemia. Patients who received LM |$75/25$| appear to have higher volatility or higher risk of hypoglycemia compared to those who received glargine. For the oral antihyperglycemic drugs, patients who received only thiazolidinedione appear to have lower volatility or lower risk of hypoglycemia compared to those who received both thiazolidinedione and sulfonylurea; patients who received only sulfonylurea are not significantly different from those who received both.

In the upper reflection barrier model, fewer covariates are significant and they are a subset of those that are significant in the volatility model. Patients with higher baseline fasting blood glucose level, higher fasting insulin, and lower BMI are associated higher reflection barrier and, hence, lower risk of hypoglycemia. For the oral antihyperglycemic drugs, patients who received only thiazolidinedione appear to have higher reflection barrier and, hence, lower risk of hypoglycemia compared to those who received both thiazolidinedione and sulfonylurea. Interestingly, the effect of baseline fasting blood glucose level, fasting insulin, and received only thiazolidinedione, the coefficients have the same direction on the risk of hypoglycemia in the models for volatility and upper reflection barrier (with opposite coefficient signs). In contrast, BMI is significant in affecting both the volatility and the upper reflecting barrier, but with opposite directions (with the same coefficient signs). That is, the overall effect of BMI on the risk of hypoglycemia is complicated by lowering the volatility (or lowering the risk) while decreasing the upper reflection barrier (or increasing the risk). This discovery has not been reported in the quantile regression analysis of Ma et al. (2021). But the overall effects of baseline BMI are worth further investigating. Figure 2 gives a set of the FHT distributions with the fitted parameters of the IF model using different levels of baseline BMI and frailties. From Fig. 2, while individual differences due to large frailty have a greater impact than BMI, smaller BMI is consistently associated with a higher risk of hypoglycemic events within the same frailty level.

The proportional hazards model of gap time between recurrent events (Huang and Chen 2003) is considered as a comparison method. Before applying the model to the DURABLE data, the censored gap for patients with more than one observation has been excluded. The estimates of the comparative model are also given in Table 5. With the exception of the covariate Adiponectin, the other covariates exhibit similar levels of significance between proportional hazards model and volatility in the proposed model. For Adiponectin, the estimates in the volatility and upper reflection barrier model are |$0.046[0.011,0.085]$| and |$0.030[-0.035,0.088]$|⁠, respectively, in the proposed model. This suggests that Adiponectin contributes in different directions on volatility and upper reflection barrier model for the recurrent risk. Given that the volatility component exerts a stronger influence on recurrence, our estimate indicates that the higher value of Adiponectin is associated with an increased risk of recurrence. In proportional hazards model, the estimate of Adiponectin is |$0.042[-0.008,0.093]$|⁠, aligning with the trend observed in the proposed model. For the other significant covariates, the signs of estimates between the volatility model and proportional hazards model are consistent. This observation suggests that, for this dataset, the proposed model offers a more detailed characterization of the concealed glucose levels linked to recurrent events.

The 2 IFs in the volatility and the upper reflection barrier capture much of the heterogeneity among the patients beyond the covariates. The variances of the 2 frailties are estimated to be far away from zero. From the fitted model, the range of volatility spans from 0.65 to 9.52 with median 2.75 and mean 3.03; the range of upper reflection barrier spans from 12.69 to 83.47 with median 29.57 and mean 29.88. The model without the frailties fits much poorer in terms of DIC and LPML (not reported).

The sensitivity analyses on the DURABLE data, using alternative starting points |$x_{0}$| and priors, including both noninformative and informative priors, were conducted. The results indicate that the significance of the covariate effects is stable when using different noninformative priors and starting points. However, the significance of certain covariate effects is sensitive to certain informative priors, as detailed in Section S3.

6. DISCUSSION

The risk of hypoglycemia is an important concern in diabetes management. It is natural to model the underlying blood glucose level as hypoglycemia occurs when it hits a lower boundary and hyperglycemia occurs when it hits an upper boundary. Because of the unique setting where hyperglycemia cannot be reliably observed in self-reported data, it is challenging to model the blood glucose level as a stochastic process. The proposed Brownian motion model with an upper reflection barrier allows bypassing the need for observing hyperglycemic event times. Only hypoglycemic times are needed for the model fitting. This model fitting is made possible by the FHT density and distribution (Hu et al. 2012). The recurrence of the hypoglycemic events is captured by a sequence of stochastic processes reaching the lower boundary. The upper reflection barrier and volatility of the reflected Brownian motion are linked to patient-level covariates and frailties. Due to the unobserved frailties, we resorted to Bayesian inference for the parameters with MCMC implemented with NIMBLE (de Valpine et al. 2017). The computation of our work relies on an accurate implementation of the FHT density/distribution functions as well as the rejection sampling algorithm with the 3-piece proposal kernel. Another computation challenge is the complexity brought by unobserved frailties in calculating model selection criteria DIC and LPML. We applied Monte Carlo integration for an approximation, with which the 2 criteria were shown to be reasonably effective in selecting the correct frailty model by simulation studies.

It is worthwhile to revisit the key model assumptions. In our model, the imposed upper reflection barrier reflects the expectation that the blood sugar level is bounded and will not reach infinity. This modeling choice aids in illustrating that glucose levels will ultimately be controlled (dropping down) in the range. We believe this assumption is reasonable, especially for participants who are consistently under regular blood sugar maintenance. On the other hand, we acknowledge that the reflecting nature of the barrier may be debatable, as blood sugar levels might remain elevated for a period rather than immediately reflecting back. Nonetheless, this reflecting approach offers a straightforward modeling strategy for this analysis. Another key assumption is that the glucose level returns to normal immediately following the occurrence of the hypoglycemia event. In reality, it certainly takes time to consume food and bring back the glucose level, but the recovery is usually in minutes and, thus, quick enough so we neglect the time used. We also assumed that, after a hypoglycemia event, the glucose level restarts at a fixed level |$x_{0}$|⁠. A sensitivity analysis with different values of |$x_{0}$| showed little difference in the resulting regression coefficient estimates.

Under the model framework, a subject with larger volatility and lower upper reflection barrier is associated with higher risk of the hypoglycemia event. In our experience, the covariate with the same direction on the risk of hypoglycemia (with opposite coefficient signs in volatility and upper reflection barrier) is usually consistent with that of the classic gap time models. On the other hand, if a covariate contributes a positive (negative) effect to both the volatility and upper reflection barrier, its overall impact is less straightforward. In this case, it is possible that the covariate shows no significant impact from the classic gap time models, while play an import role for volatility or upper reflection barrier. An example of this is the covariate Adiponectin in our data application. Therefore, despite its complexity, the proposed model provides additional insights for a better understanding of the data. To further investigate the covariate overall effect of proposed model, we recommend plotting the FHT distribution to provide an overall characterization of the impact of this covariate.

Several directions are worth further investigation. In the broad sense, the proposed model can be applied to scenarios where an event occurs when an underlying health level process hits a boundary on one side. Therefore, many of the examples based on Wiener process reviewed in Lee and Whitmore (2006) are potentially applicable. The Wiener process approach ensures finiteness of the FHT by a nonzero drift. Our process does so with a reflecting boundary on the other side for a driftless Wiener process. When the reflecting boundary is removed, the FHT has a positive probability of being infinity, which makes it applicable when a cure rate is needed (Lee and Whitmore 2006, Section 5). Furthermore, when a large number of unobserved frailties are included in the model, the MCMC chains for some parameters exhibit high autocorrelation, resulting in a relatively small effective sample size. This suggests that there is potential for improving the algorithm. Finally, the DURABLE dataset has additional longitudinally observed blood glucose levels. To combine these longitudinal observations with the hypoglycemic events into a joint modeling framework, the transition density of a reflected Brownian motion would be needed. Incorporating this density into our framework would be interesting but not trivial.

SUPPLEMENTARY MATERIAL

Supplementary material is available at Biostatistics Journal online.

FUNDING

None declared.

CONFLICT OF INTEREST

None declared.

APPENDIX

A. TAIL OF FHT DENSITY

Here we show that the right tail of the FHT is bounded by an exponential rate. Note that rates |$\lambda_{n}$| monotonically increase to |$\infty$| as |$n\to\infty$| (so, |$\lambda_{1}$| is the slowest rate), |$c_{1} \gt 0$|⁠, and |$|c_{n}| \lt 1$| for |$n\geq 2$|⁠. It can be shown that density |$f(t)$| is asymptotically equivalent to |$c_{1}\lambda_{1}e^{-\lambda_{1}t}$| as |$t\to\infty$|⁠.

First, let us rewrite the FHT density as follows:

where |$\lambda_{n}=b(2n-1)^{2}$| and |$b=\frac{\sigma^{2}\pi^{2}}{8(\kappa-\nu)^{2}}$|⁠.

Now, consider |$h(t)=\sum_{n\,=\,2}^{\infty}c_{n}\lambda_{n}e^{-(\lambda_{n}-\lambda_{1})t}$| and observe that

Therefore, if |$\sum_{n\,=\,2}^{\infty}\lambda_{n}e^{-(\lambda_{n}-\lambda_{1}-b)t}$| is bounded for all sufficiently large |$t$|⁠, then |$|h(t)|\to 0$| as |$t\to\infty$|⁠. Indeed, since |$\lambda_{n}-\lambda_{1}-b\,=\,b(2n-1)^{2}-2b \gt 0$| for |$n\geq 2$|⁠, for |$t\,\gt\,1$| we have

because |$\sum_{n\,=\,2}^{\infty}\lambda_{n}e^{-\lambda_{n}}$| is obviously a convergent series. Thus, hitting time density |$f(t)\sim c_{1}\lambda_{1}e^{-\lambda_{1}t}$| as |$t\to\infty$|⁠.

This result allows the use of an exponential distribution, with proper scaling, on the right tail to bound the FHT density as detailed next.

B. REJECTION SAMPLING ALGORITHM

To sample from the FHT density |$f$|⁠, we handle the left tail, body, and right tail separately. Define |$g_{1}(t)$|⁠, |$g_{2}(t)$|⁠, and |$g_{3}(t)$| as the proposal density for the left, body, and right components, respectively,

where |$k_{1}=f(t_{m})/t_{m}$| and |$k_{2}=[f(t_{m})-f(t_{q})]/(t_{m}-t_{q})$|⁠. For illustration, Fig. B1 shows the actual (target) density of first hitting distribution in red, the kennels of the 3 proposal densities in grey, and the envelopes derived by multiplying corresponding constants to the kernels in black. Given the shape properties of |$f$|⁠, the |$i$|th component of |$f$| can be bounded by |$M_{i}g_{i}(t)$|⁠, where |$M_{i}$| can be identified by maximizing |$f(t)/g_{i}(t)$| over the domain of |$g_{i}$|⁠, |$i\,=\,1,2,3$|⁠.

The rejection sampling algorithm for generating one observation from |$f$| is given in Algorithm 1.

References