ushr captures dynamics of previously published data
We first demonstrate the general capability of ushr by using the package to analyze publicly available data from the ACTG 315 clinical trial. These data have been described previously [5, 15, 16], and are available with the package, or at https://sph.uth.edu/divisions/biostatistics/wu/datasets/ACTG315LongitudinalDataViralLoad.htm(accessed 15 September 2019). Briefly, longitudinal measurements of HIV viral load are documented for 46 chronically-infected adults undergoing RTI/PI-based ART. The assay detection threshold was 100 RNA copies ml −1, and data were recorded up to 28 weeks following treatment initiation.
Given such a dataset, one can use ushr to process the data and fit the applicable mathematical model (here, biphasic or single phase). First, individual trajectories that adhere to the specific inclusion criteria are identified i.e. those that form a consistent decreasing sequence towards the detection threshold, and include the minimum number of observations required for model fitting (six for the biphasic model and three for the single phase model). Our algorithm then fits the appropriate model to the data and computes the best-fit parameter estimates with 95% confidence intervals. Note that the single phase model is only applied to individuals ineligible for the biphasic model. Both the data processing and model fitting steps can be achieved in one line of code using the ushr() function (or ushr_triphasic() if fitting the triphasic model). The resulting model trajectories can be visualized using plot_model() (Fig 1 and Additional file 2).
For the ACTG 315 data, twelve individuals were fit with the biphasic model (Fig 2), and four with the single phase model (Fig 3). Although some single phase subjects had sufficient data to fit the biphasic model (i.e. at least six observations), the resulting 95% parameter confidence intervals were either unattainable or sufficiently wide to indicate an unreliable fit. This can occur, for example, when one of the decay phases is poorly documented (i.e. has few data points). As a result, the subjects were re-fit with the single phase model. This re-fitting step is automated in the package, although the user can control the confidence interval width above which a biphasic fit is deemed unreliable.
A summary of the model applied to each subject, and the associated infected cell lifespans, can be obtained using summarize_model() (Additional file 2). The biphasic fits from the ACTG 315 data had median lifespans of 2.1 days (SD = 0.8) and 25.6 days (SD = 11.9) for short and long-lived infected cells, respectively. These estimates are in line with previous analyses of these data; for example, 2.3 days and 31.3 days, respectively [5]. The single phase fits had a median lifespan of 9.1 days (SD = 5.8). To assess the robustness of parameter estimates for each subject, one can view their corresponding 95% confidence intervals (Additional file 2). Pairwise parameter correlation plots can also be viewed to assess dependencies at the population-level. We present further analyses comparing parameter estimates with their true values below.
Finally, one can calculate the TTS for all eligible subjects with the get_TTS() function. The user can specify the threshold viral load under which a subject is defined to have reached suppression. Here we set the suppression threshold equal to the detection threshold. The user can also specify whether the parametric or nonparametric method should be used. For the ACTG 315 data, the median TTS estimates were 65.7 (SD = 31) and 69.8 days (SD = 37.9) for the parametric and nonparametric methods, respectively (Additional file 2).
Parameter estimates are accurate for sufficiently high resolution data
Study sampling design can have a substantial impact on the reliability of parameter estimates obtained from model fitting. For example, short-lived cells have an average lifespan on the order of days, which may be difficult to estimate accurately in a clinical setting, given the challenge of obtaining frequent observations. Indeed, one may expect to overestimate this lifespan, as the most reliable fits will tend to derive from individuals with slower decay and thus more information regarding the first decline phase.
We investigated this possibility by comparing parameter estimates derived from simulated data with different sampling frequencies (Additional file 3). Briefly, we simulated data using the simulate_data() function in ushr, and specified the number of subjects with nsubjects = 200. This command generates individual trajectories from the biphasic model (Eq. 1), using underlying parameters that are sampled from lognormal distributions with user-specified mean and standard deviation. Gaussian noise (with zero mean and user-specified standard deviation) is added on the log10 scale, and the number of observed time points is randomly sampled for each subject from a pre-defined range. This approach assumes subjects are sampled at regular intervals until either the end of the study period or they are lost to follow-up.
To investigate the impact of sampling resolution, we varied the frequency of observation times using three different scenarios. In our high resolution dataset, individuals were sampled four times per month; in the medium resolution dataset this was reduced to twice per month; and in the low resolution dataset, the frequency was reduced further to once per month. In all cases, individuals were sampled for a minimum of three months and maximum of one year.
As above, we fit the high, medium, and low resolution datasets using ushr(), then calculated an average deviation score that summarized the difference between the true and estimated parameter values. For each parameter p, and study resolution r, we defined the average deviation as
$$ D_{p,r} = \frac{1}{n_{p,r}}\sum_{i = 1}^{n_{p,r}}\left(\theta_{i} - \hat{\theta}_{i}\right)/\theta_{i}, $$
(4)
where np,r is the number of subjects in the study that were fit with the biphasic model, θi is the true parameter value for subject i, and \(\hat {\theta }_{i}\) is the estimated value. Note that for the parameters A and B, we log10-transformed θi and \(\hat {\theta _{i}}\) prior to calculating the deviance. For TTS values, we used the parametric approach to calculate θi and \(\hat {\theta }_{i}\) from the corresponding true and estimated parameters, respectively. We repeated the process for 100 different randomly sampled underlying parameter distributions, resulting in 1200 average deviation scores (4 parameters, 3 study resolutions, and 100 repetitions).
In general, estimates of the short-lived lifespan from the low resolution studies had the largest deviation (Fig 4). The values were consistently negative, indicating overestimation; and the shortest lifespans exhibited the largest deviances. These results are consistent with the intuition that the fastest decay phases are most vulnerable to infrequent sampling. All other parameter values had relatively low deviation scores, although the long-lived lifespan estimates also exhibited a slight negative bias at low sampling frequencies.
Finally, we explored the extent to which estimates from the subset of individual trajectories fit using the biphasic model could be extrapolated across the whole study population. We focussed on the short and long-lived lifespans since these exhibited the greatest deviances. For each study and lifespan, we compared the true population median to the median estimate from the fitted subset. Unsurprisingly, estimates from low resolution studies were the poorest representation of the true population-wide values (Fig 5). Both lifespans were overestimated, reflecting problems with infrequent observations, as discussed above. In addition, fewer subjects in these low resolution studies had sufficient data for fitting, leading to smaller sample sizes and a greater departure from the population average.
In contrast, estimates from the intermediate and high resolution studies were in good agreement with the true population average, although the long-lived lifespans were slightly underestimated at higher values. At intermediate and high resolutions, the effects of sample size and infrequent observations are reduced. However, our algorithm for simulating data is such that individuals with slow second decline phases are more likely to be lost to follow-up before reaching suppression. Thus our fitting may be biased towards subjects with shorter long-lived lifespans than the population average. The effect will be most pronounced when the average long-lived lifespan is high and restrictions from infrequent sampling are minimal.
Overall, these results demonstrate that our fitting procedure can reliably estimate individual parameters at intermediate and high sampling resolutions, and that these estimates represent an unbiased sample of parameters across the population. Conversely, caution must be taken when interpreting parameters at low resolutions, and in extrapolating these estimates to the wider population.
ushr performs comparably to nonlinear mixed effects modeling
In ushr we take an individual-based approach by independently fitting data from each subject that meets our inclusion criteria and our minimum number of observations. Since the model equations are pre-defined (i.e. do not need to be specified by the user), our package has the advantage of straightforward implementation. Conversely, another common method for fitting longitudinal clinical data is to use a population-based approach such as nonlinear mixed effects (NLME) modeling [5, 17]. Briefly, the NLME framework aims to describe inter-subject variation in the population by assuming a form for the distribution of parameters at the population-level. This approach has the advantage that information from all individuals who meet the inclusion criteria can be leveraged to inform the overall fits, regardless of the number of available measurements. However, existing NLME software packages are typically general purpose, and require users to manually define model structure according to their specific need. In addition, imposing population-level distributions may cause NLME to perform poorly on individuals with outlying dynamics.
In order to compare estimates obtained from NLME modeling with those obtained above from ushr, we re-fit the low, medium, and high resolution datasets using the NLME package saemix in R (Additional files 4, 5, 6, 7, and 8) [18]. In brief, we defined the parameters A,B,δ, and γ as log-normal distributions with individual random effects, and used saemix to fit Eq. 1 and estimate the corresponding population-level distributions.
In general, we found lifespan estimates from the NLME approach were similar to, or marginally worse than, those obtained from ushr for the intermediate and high resolution data (Fig 6). Conversely, for the low resolution data, the NLME lifespans were closer to the true population average, although the short-lived lifespans were slightly underestimated. In summary, our individual-based algorithm performs comparably to the population-based NLME approach, with clear discrepancies only becoming apparent at low sampling resolutions.