Time-In-Range Analyses of Functional Data Subject to Missing with Applications to Inpatient Continuous Glucose Monitoring

Time in range (TIR) is defined as the percentage of time that CGM glucose values are within a target glycemic range, denoted by $G$, over a specified amount of time, represented by the time interval $[0,\tau]$. For example, to represent the TIR over 7 days with the target glucose range between 70 and 180 mg/dL, we set $G=[70,180]$ and $\tau=7$ (days). Let $\mathcal{T}\doteq\left\{0=t_{1}<t_{2}<\dots<t_{K}=\tau\right\}$ denote the set of time points when the CGM sensor reads interstitial glucose over the time interval $[0,\tau]$ (e.g., a time grid equally spaced by 5 minutes between time $0$ and $\tau$). We formulate an individual’s CGM glucose trajectory $Y(\cdot)$ as a right continuous piecewise constant function that only jumps at the time points in $\mathcal{T}$. With notation defined above, the TIR can be represented as a statistical functional of $Y(\cdot)$:

In a CGM study conducted in hospital, insufficient CGM sampling often occurs due to a short hospital stay. Let $C$ denote the follow-up duration of CGM (i.e., time from the beginning to the end of CGM). If $C<\tau$, then $Y(t)$ with $t>C$ is not observed, resulting in monotone missing of the CGM trajectory. Let $\delta^{*}(t)$ be a right continuous piecewise constant function. We use $\delta^{*}(t)$ to account for the intermittent missing of the CGM trajectory due to reasons such as random device errors. Specifically, for $t<C$, we have $\delta^{*}(t)=1$ when a CGM reading is available in the time interval $\left[t_{\iota(t)},t_{\iota(t)+1}\right)$, where $\iota(t)=\text{max}\{j:t_{j}\leq t<C,1\leq j\leq K\}$ and 0 otherwise. Define $\delta(t)=\delta^{*}(t)I(t\leq C)$.

In addition to the CGM glucose trajectory $Y(\cdot)$, other information may also be observed or collected in an inpatient CGM study, for example, individuals’ demographics and clinical characteristics. We capture such information by a $p$-dimensional vector of covariates, $\boldsymbol{Z}(\cdot)$, which is allowed to change over time. The observed data then include $n$ independent and identically distributed (i.i.d.) replicates of $\left\{\delta(t),C,\delta(t)Y(t),\boldsymbol{Z}(t):t\in(0,\tau]\right\}$, denoted by $\left\{\delta_{i}(t),C_{i},\delta_{i}(t)Y_{i}(t),\boldsymbol{Z_{i}}(t):t\in[0,\tau]\right\}_{i=1}^{n}$.

Mean TIR, namely $\mu_{W}\doteq E(W)$, is often used to summarize TIR at the population level and has been commonly used as a primary endpoint to evaluate group differences in CGM studies. When all CGM glucose trajectories (i.e., $Y_{i}(\cdot),\ i=1,\ldots,n$) are completely observed over the time interval $[0,\tau]$, the mean TIR, $\mu_{W}$, can be readily estimated by its empirical counterpart, namely, $n^{-1}\sum_{i=1}^{n}\left(\int_{0}^{\tau}I(Y_{i}(t)\in G)dt\right)/\tau$.

In practice, some $Y_{i}(\cdot)$’s are only partially observed, as often encountered in an inpatient CGM study. In the routine TIR analysis, one may estimate $\mu_{W}$ by simply averaging the within-individual proportion of the observed CGM readings belonging to $G$. This is equivalent to adopting a naive estimator of $\mu_{W}$, which is given by

where $\widecheck{W}_{i}=\{\int_{0}^{\tau}I(Y_{i}(t)\in G)\delta_{i}(t)dt\}/\{\int_{0}^{\tau}\delta_{i}(t)dt\}$. However, the glucose trajectory of a hospitalized patient is usually not stationary over time (see examples in Figure 2(A)). As a result, $\widecheck{W}_{i}$, which is obtained from a “snap shot” over a shorter time period with $\delta_{i}(t)=1$, may fail to adequately reflect $W_{i}$, the underlying TIR of individual $i$. Intuitively, this would lead to biased estimation of $\mu_{W}$. This is confirmed by our simulation studies (see Table 1).

To address the issue caused by the incompletely observed $Y_{i}(\cdot)$’s, we propose an estimator of $\mu_{W}$ by discovering the following probabilistic representation of the mean TIR:

where $p_{G}(t)=\Pr\left(Y(t)\in G\right)$, representing the probability of CGM glucose falling into the target range at time $t$. This representation of $\mu_{W}$ enlightens a viable direction to estimate $\mu_{W}$ based on an estimate for $p_{G}(t)$. Adopting this strategy can lead to multiple benefits: (i) $p_{G}(t)$ is a population quantity depending on $Y(t)$ cross-sectionally; thus tackling $p_{G}(t)$ (instead of $W_{i}$) allows for information sharing across subjects and circumventing the difficulty in directly imputing the missing data in within-subject glucose trajectories; (ii) with a proper estimate for $p_{G}(t)$ bounded between $0$ and $1$, the corresponding plug-in estimate for $\mu_{W}$ automatically satisfies the bounded constraint for $\mu_{W}$.

To estimate $p_{G}(t)$, we need to properly account for the possible missingness of $Y(t)$’s at time $t$. To tackle this problem, we adopt the following assumptions:

(M1): $\delta^{*}(t)$ is independent of $C$ and $Y(t)$;

(M2): $C$ is independent of $Y(t)$ conditional on ${\mathcal{H}}(t)$, where ${\mathcal{H}}(t)$ denotes the subject-specific history data before time $t$.

Assumption (M1) is reasonable in CGM studies because the intermittent missing is mostly caused by some external technical issues with the CGM device and is less likely related to the underlying glycemic control. Assumption (M2) admits a realistic missing-at-random (MAR) mechanism for the monotone missingness related to $I(t\leq C)$. By this assumption, $I(t\leq C)$ is independent of the current glucose value $Y(t)$ conditional on ${\mathcal{H}}(t)$, while $I(t\leq C)$ is allowed to depend on the glucose history captured by ${\mathcal{H}}(t)$, thereby accommodating realistic scenarios, such as better glycemic control leading to a shorter hospital stay.

We can derive an estimator of $p_{G}(t)$ by applying the idea of inverse probability weighting. Specifically, under assumptions (M1) and (M2), we can show that

Similarly, we can get $E\left\{\frac{\delta^{*}(t)I(t\leq C)}{\Pr(t\leq C|\mathcal{H}(t))}\right\}=E\{\delta^{*}(t)\}.$ By these results, we propose the following estimator of $p_{G}(t)$:

where $\widehat{p}_{C,i}(t)$ is a reasonable estimator of $p_{C,i}(t)\doteq\Pr(t\leq C_{i}|{\mathcal{H}}_{i}(t))$ and ${\mathcal{H}}_{i}(t)$ is the sample analogue of ${\mathcal{H}}(t)$ (i.e., the history data of subject $i$ at time $t$).

Remark 1: Consider a special case where $C$ is independent of both $Y(t)$ and ${\mathcal{H}}(t)$ for all $t\in[0,\tau]$. This points to a practical scenario where the termination of CGM is completely independent of the underlying CGM glucose trajectory. In this special case, the proposed estimator reduces to a simplified estimator,

The remaining task is to obtain $\widehat{p}_{C,i}(t)$. To this end, we assume that the distribution of $C$ is influenced by ${\mathcal{H}}(t)$ through the covariates in $\boldsymbol{Z}(t)$, which may include summaries of glucose history up to time $t$ (e.g., average or maximum CGM glucose before time $t$). Following this direction, we further model the relationship between $C$ and $\boldsymbol{Z}(t)$. A popular choice of such a model is the semiparametric Cox proportional hazard model, which assumes

where $\lambda_{C}(t|\mathcal{H}(t))$ denotes the conditional hazard function of $C$ given ${\mathcal{H}}(t)$, $\lambda_{0}(t)$ denotes an unspecified baseline hazard function, and $\boldsymbol{\beta}_{0}$ is a vector of unknown coefficients. Let $\Lambda_{0}(t)=\int_{0}^{t}\lambda(u)du$. Using the observed data on $C$ and $\boldsymbol{Z}(t)$, we can obtain the partial likelihood estimator of $\boldsymbol{\beta}_{0}$, denoted by $\widehat{\boldsymbol{\beta}}$, and the Breslow estimator of $\Lambda_{0}(t)$, denoted by $\widehat{\Lambda}(t)$ (Andersen and Gill, 1982; Lin, 2007). The corresponding estimator of $p_{C,i}(t)$ is then given by

Plugging the $\widehat{p}_{G}(t)$ obtained based on (3) and (5) into the representation of $\mu_{W}$ in (2) leads to the proposed estimator of $\mu_{W}$,

Note that $\widehat{p}_{G}(t)$ is a piecewise constant function. Thus, the integration involved in $\widehat{\mu}_{W}$ can be evaluated via finite summation.

The proposed estimation of $\mu_{W}$ can readily accommodate other estimators of $p_{C,i}(t)$ obtained under different modeling of $C_{i}$ given ${\mathcal{H}}_{i}(t)$, such as the additive hazards model (Aalen, 1980). More specifically, we may consider a general setting where

Here $\boldsymbol{\eta}(\boldsymbol{Z},t)=\{\boldsymbol{Z}(u):\ 0\leq u\leq t\}$ and $\boldsymbol{\theta}_{0}$ captures either real valued or function valued model parameters. It is easy to see that the Cox model (4) is a special case with $\Psi(\boldsymbol{\eta}(\boldsymbol{Z},t),t|\boldsymbol{\theta}_{0})=\exp\left\{-\int_{0}^{t}\exp\{\boldsymbol{Z}(u)^{\top}\boldsymbol{\beta}_{0}\}d\Lambda_{0}(u)\right\}$. The estimator of $p_{C,i}(t)$ can be generally formulated as $\Psi(\bar{\boldsymbol{Z}}_{i}(t),t|\widehat{\boldsymbol{\theta}})$, where $\widehat{\boldsymbol{\theta}}$ is an existing estimator of $\boldsymbol{\theta}_{0}$,

Our data generation scheme closely mimic real data scenarios in inpatient CGM studies. Specifically, we generate CGM glucose trajectories as Gaussian processes, $\mathcal{GP}(\mu(t),k(t,t^{\prime}))$, where $\mu(t)$ denotes the mean function and $k(t,t^{\prime})$ denotes the covariance function characterized by a periodic kernel that takes the form, $k(t,t^{\prime})=\sigma^{2}\text{exp}\left(-\frac{2}{l^{2}}\text{sin}^{2}\left(\pi\frac{|t-t^{\prime}|}{p}\right)\right)$. The kernel parameters, amplitude ($\sigma$), length-scale ($l$), and period ($p$), govern the mean deviation, trajectory smoothness, and repetition frequency, respectively. To reflect the natural 1-day (i.e., 1440 minutes) periodicity in glucose, we set the the period $p=1440$, length-scale $l=1$, and amplitude $\sigma=62$. The resulting covariance matrix is shown in Figure 1(B).

In each simulation setting, three groups of CGM glucose trajectories are generated and are referred to as Group 1, Group 2, and Group 3. For Groups 1 and 2, the corresponding mean functions, denoted by $\mu_{1}(t)$ and $\mu_{2}(t)$, are specified as the empirical mean functions of the CGM glucose trajectories observed for the two study groups considered in the real data application presented in Section 5. For Group 3, we use the same mean function specified for Group 1 so that comparing Group 1 versus Group 3 can serve as the null case with respective to the hypothesis testing on the equivalence in mean TIR. In addition, we set ${\mathcal{T}}$ as an equally spaced time grid with the step size of 5 minutes, mimicking the data capturing scheme of a Dexcom G6 CGM system.

To induce missing data for the generated CGM glucose trajectories, we first impose intermittent missing. Specifically, we draw the starting time of the intermittent missing from an exponential distribution with the scale parameter of 3424, which is derived from the real data application discussed in Section 5, and then set the missing duration as a random number generated from a uniform distribution between 10 and 70 (minutes). By doing so, we determine the time interval when the intermittent missing occurs, and subsequently $\delta_{i}^{*}(t)$.

Next, we generate $C_{i}$, which represents the ending time of CGM, in the following two ways:
(I) $C_{i}$ is completely non-informative and is sampled from a distribution not related to the generation of glucose trajectories. Specifically, for Groups 1 and 2, we sample $C_{i}$’s with replacement from the empirical distribution of the observed CGM follow-up durations in the real data application discussed in Section 5. For Group 3, $C_{i}$’s are generated from a mixture distribution, $0.8F_{1}(c)+0.2F_{2}(c)$, where $F_{1}(\cdot)$ is the uniform distribution function between 0 to 2 (days) and $F_{2}(\cdot)$ is the uniform distribution function between 2 and 9 (days).
(II) $C_{i}$ is informative, depending on the history of glucose trajectories, and is generated from the Cox’s proportional hazards model (4) with $\boldsymbol{Z}_{i}(t)=(Z_{i}^{(1)}(t),\ Z^{(2)}_{i})^{\top}$, where $Z^{(1)}_{i}(t)$ captures the average glucose during the previous day as of time $t$ (which equals $0$ when $t$ indicates a time point during the first day), and $Z^{(2)}_{i}$ is drawn from a uniform distribution $\text{Unif}(-0.5,0.5)$. To make $Z^{(1)}_{i}(t)$ and $Z^{(2)}_{i}$ to range within a similar scale, we normalize $Z^{(1)}_{i}(t)$ by dividing it by $100$. We let $\lambda_{0}(t)=\frac{1}{16}t^{-15/16}$ and $\boldsymbol{\beta}_{0}=(-2,-2)^{\top}$ for Group 1, $\lambda_{0}(t)=\frac{1}{4}t^{-3/4}$ and $\boldsymbol{\beta}_{0}=(2,2)^{\top}$ for Group 2, and $\lambda_{0}(t)=\frac{1}{2}t^{-1/2}$ and $\boldsymbol{\beta}_{0}=(2,2)^{\top}$ for Group 3.

With the $\delta^{*}(t)$ and $C$ obtained as above, we compute the availability indicator $\delta(t)=\delta^{*}(t)I(C\leq t)$, and $Y(t)$ is observed if and only if $\delta(t)=1$. The empirical distribution of $C$ and the empirical mean of $\delta(t)=1$, which reflects the proportion of the available glucose readings at time $t$, in case (I) and case (II) are plotted in Figure 1(C)–Figure 1(F).

In each simulation setting, we consider sample sizes $n=50,100$ and $200$, and generate 1000 replicated datasets. The bootstrap based inferences are implemented based on 200 bootstrapping samples.

In each simulation setting, we analyze the mean TIR over 7 days with the target glucose ranges, $<70$ mg/dL, between 70 and 180 mg/dL, and $>180$ mg/dL, using three different methods: (1) the oracle approach that computes $n^{-1}\sum_{i=1}^{n}\left(\int_{0}^{\tau}I(Y_{i}(t)\in G)dt\right)/\tau$ based on the fully observed $Y_{i}(\cdot)$ without any missing data; (2) the proposed approach that uses $\widehat{\mu}_{W}$; (3) the naive approach that completely ignores the missing data and estimates $\mu_{W}$ by $\widecheck{\mu}_{W}$.

Table 1 summarizes the simulation results in the setting with informative CGM follow-up durations (i.e., $C_{i}$’s). The presented results include the empirical averages, relative biases, empirical standard deviations, and estimated standard errors (based on bootstrapping) of the three different estimators of mean TIR. The relative bias is defined as the absolute difference between the average estimate and the average oracle estimate divided by the average oracle estimate. We also report the empirical rejection rates from testing the mean TIR equivalence between Groups 1 and 2 and between Groups 1 and 3. Since the mean TIRs are the same between Groups 1 and 3 and are different between Groups 1 and 2, the empirical rejection rates in these cases respectively reflect empirical size and empirical power.

From Table 1, we observe that the proposed estimator of mean TIR is virtually unbiased with the relative biases mostly below 5% and exhibiting a decreasing trend with the sample size. The empirical standard deviations and the estimated standard errors closely agree with each other. The naive estimator demonstrates apparent biases in many cases. For example, even when the sample size is as large as 200, the relative biases of the naive estimates for the TIR $<70$ mg/dL are still pretty substantial, equal to 23.7%, 18.7%, and 33% for Groups 1, 2, and 3 respectively. In contrast, the corresponding relative biases of the proposed estimator are much smaller, equal to 0.4%, 0.1%, and 6.6%. It is noted that the proposed estimator generally produces larger relative biases for Group 3 as compared to that for Group 1, though Group 1 and Group 3 share the same mean TIR. This can be explained by the different distributions of $C$ between these two groups (see Figure 1(D)), which entail more frequent early termination of CGM and thus more intensive monotone missing in the data for Group 3 (as confirmed by Figure 1(F)). As expected, the lower proportions of available CGM glucose readings in Group 3 (as compared to Group 1) also result in larger estimation variability captured by either empirical standard deviations or estimated standard errors.

Regarding the hypothesis testing that compares the mean TIR between Groups 1 and 3, the proposed method yields empirical sizes quite close to the nominal significance level of 5%, while the naive approach may be prone to inflated type I errors (e.g., the empirical size of 10% when comparing the mean TIR $>180$ mg/dL with $n=200$). By examining the empirical rejection rates for comparing Groups 1 vs. 2, we find that the test based on the proposed method is much more powerful as compared to the test based on the naive approach, as evidenced by more than two fold larger empirical rejection rates. For example, with the larger sample size 200, the naive method only has around 11% power to detect the difference in mean TIR between Group 1 and Group 2, while the proposed method has 74% power to do so.

The simulation results in the setting with non-informative $C_{i}$’s are presented in Table S1 of Web Appendix B.1. The finding from Table S1 is consistent with that from Table 1. Both suggest that failing to properly account for the special missing features of the inpatient CGM glucose trajectories can lead to considerably biased estimation of mean TIR as well as problematic inferences that are reflected by the inflated type I errors and substantial power loss associated with the naive method.

We conduct simulation studies under other configurations that represent different realistic characteristics of inpatient CGM data. For example, we consider settings with different sizes of group difference by modulating the mean functions used to generate CGM glucose trajectories. We also simulate data with different missing data patterns by varying the distributions of $C$. Details on these additional simulation studies are provided in Web Appendix B.2. The results consistently demonstrate the advantages of the proposed method over the naive method.

We further conduct simulation studies to evaluate the robustness of the proposed method to the misspecification of the assumed Cox model for $C$ given ${\mathcal{H}}(t)$. Specifically, we generate the CGM glucose trajectory $Y_{i}(\cdot)$ as a Gaussian process with the mean function $\mu^{[i]}(t)$ and the covariance function $k(t,t^{\prime})$, which is the same as that plotted in Figure 1(B). We set the mean function $\mu^{[i]}(t)=\mu_{1}(t)-3\zeta_{i}$ for Group 1 or 3 and $\mu^{[i]}(t)=\mu_{2}(t)+3\zeta_{i}$ for Group 2, where $\mu_{1}(t)$ and $\mu_{2}(t)$ are the same as those plotted in Figure 1(A) and $\zeta_{i}$ is a $Bernoulli(0.5)$ random variable. We generate the ending time of CGM, $C_{i}$, based on a transformation model (Cheng et al., 1995), $C_{i}=s\cdot\exp(-\zeta_{i}+\epsilon_{i})$, where $s=8$ for Group 1 and Group 2 and $s=3$ for Group 3, and the error term $\epsilon_{i}$ has the distribution function, $pF_{3}(x)+(1-p)F_{4}(x)$. Here $F_{3}(\cdot)$ denotes the standard logistic distribution function and $F_{4}(\cdot)$ denotes the extreme value distribution function. It is clear that $C_{i}$ follows the Cox proportional hazards model (4) with $\boldsymbol{Z}_{i}(t)=\zeta_{i}$ when $p=0$ and a proportional odds model when $p=1$. As such, the magnitude of the parameter $p$ reflects the degree of departure from a Cox model.

Table 2 presents the simulation results for the settings described above with $p=0.5$ and $p=1$. It is shown that the proposed method yields quite small relative biases, ranging from 1.2% to 2.0%, when $p=0.5$ or $1$. The empirical sizes for testing the mean TIR equivalence between Group 1 and Group 3 are still quite close to the nominal significance level of $5\%$. At the same time, the relative biases of the naive estimator are several fold larger than those of the proposed estimator. The naive method also results in tests which yield inflated type-I errors or suffer from substantial power loss. These observations suggest that the proposed method can offer robust estimation and inference on mean TIR and consistently outperform the naive method even when the Cox model assumption for the CGM follow up duration (i.e., $C_{i}$) does not hold.