Repeated Time-to-event Analysis of Consecutive Analgesic Events in Postoperative Pain

PATIENT-controlled analgesia is generally accepted as a safe and efficient way to individualize dosing with opioids in clinical postoperative pain management.^1,2  However, it has frequently been suggested that combining opioids with adjuvant analgesics such as gabapentin or that dosing before surgery (preemptive treatment) could lead to pain relief with less opioid consumption and opioid-related side effects.^3–5  Pain intensity is the primary measure that defines whether a drug has analgesic efficacy, but the measure of opioid consumption (or rescue medication) itself can be the primary endpoint when opioid-sparing interventions are investigated.⁶ 

Analgesic consumption has traditionally been reported as a mean with SD in clinical trials.^5,7  However, analgesics are given in discrete doses, and analgesic consumption data are therefore not continuous and have been shown not to be normally distributed among individuals.^8,9  A number of parametric and nonparametric statistical tests have been used to compare analgesic consumption between intervention groups, but the discrete nature of analgesic consumption data violates assumptions of most traditional statistical tests.⁸  Meta-analyses frequently analyze weighted mean reduction in 24-h analgesic consumption, and significant reduction is often considered indicative of opioid-sparing efficacy.^5,7 

The timing of first analgesic event after surgery has also been studied as a measure of analgesic efficacy, for example, by using time-to-event (TTE) analysis.^5,10  However, for postoperative conditions where multiple analgesics are needed, TTE of a single event involve censoring of data and loss of information. Furthermore, it remains debated whether an effect on the timing of first analgesic event is relevant to clinical pain management.^5,11 

Concerns have been raised that the traditional ways of analyzing analgesic consumption do not sufficiently handle time-varying factors such as nausea and vomiting and that a direct causal relation between pain intensity and analgesic consumption remains to be demonstrated.^12,13 

Despite the abundance of clinical studies investigating the value of opioid-sparing regimens with adjuvant analgesics, large discrepancies are found between studies.^7,14  In addition, the efficacy and dose-exposure–response relations of opioid-sparing effects of adjuvant analgesics in combination therapy generally remain unclear. Inappropriate analysis methods of analgesic consumption may be one important factor in relation to the difficulties in demonstrating efficacy and dose-exposure–response relations.⁹ 

We propose that some or all of these issues with the discrete nature of analgesic consumption, the time-varying factors, and the difficulties in demonstrating efficacy and dose-exposure–response relations may be solved by applying model-based analysis of the timing of consecutive analgesic events after surgery. Repeated time-to-event (RTTE) modeling is an extension to parametric TTE survival analysis using nonlinear mixed-effects modeling^15–18  and builds on methods that are well described in the statistical literature.^19,20  Methods to study RTTE have been implemented in other fields of health science and have proven superior to TTE for investigation of recurrent functional disability.²¹  To our knowledge, RTTE modeling has not previously been applied to study analgesic events in postoperative pain. It is the aim of this study to demonstrate how RTTE modeling can be applied to analyze consecutive analgesic events in postoperative pain.

The study was approved by the Danish Data Protection Agency, Copenhagen, Denmark (J. nr.: 2008-58-0028). Retrospective data were collected from medical records of patients admitted for surgery after hip fracture at Orthopaedic Ward, Aalborg University Hospital, Aalborg, Denmark, in the period from May to December 2012.

Patients were included from a population of 109 randomly selected patients and included if they received hospital standard analgesic treatment consisting of 1 g acetaminophen four times daily plus morphine on request. Patients who had either none or multiple surgeries registered were excluded. The sample size was chosen for illustrative purposes and was not estimated from power calculations.

Times of morphine dosing (analgesic events) were extracted from the medical records. Data were collected until 96 h after surgery or until censoring due to discharge from hospital or initiation of nonstandard analgesic treatment.

An RTTE model was developed to describe the probability of analgesic events over time using the Laplace estimation method in NONMEM^®, version 7.2 (ICON Development Solutions, USA).^15,22  RTTE modeling was performed as described by Karlsson et al.¹⁵  and Plan et al.¹⁷ 

The probability density for an event to occur was formulated as follows¹⁵ :

where t is the time after surgery, S(t) is the survival function describing the probability of not having an event in the time interval between t_j−1 and t_j (t_j−1 being the time of surgery or the time of last event), and h(t) is the hazard function. Censoring, defined as the last time point a patient was in the study, was set at time of dropout or 96 h after surgery.

Hazard (h) was modeled over time (t) as follows:

where h₀(t) is the baseline hazard and β_i is the coefficient describing the effect of covariate X_i, and η is a random effect (frailty) describing a log-normal distribution of the hazard in the population. β_i was modeled as (β_i · X_i) for dichotomous risk factors and as a change from the median risk factor (β_i · [X_i − X_{i, median}]) for continuous risk factors.

For an exponential distribution, h₀(t) was parameterized as follows:

where λ is the scale parameter of the exponential distribution.

For a Weibull distribution, h₀(t) was parameterized as follows:

where λ and γ are the scale and shape factor of the Weibull distribution, respectively.

For a Gompertz distribution, h₀(t) was parameterized as follows:

where λ and γ are the scale and shape factor of the Gompertz distribution, respectively.

The hazard model with the most appropriate hazard distribution was first identified. Subsequently, the final model was found by implementing covariates as proportional hazards. Night, defined as the period between 11 pm and 7 amversus day (between 7 am and 11 pm), was implemented as a time-varying discrete covariate. Sex was implemented as a discrete covariate on hazard. Age was implemented as a continuous covariate on hazard and both a linear and a power function scaled to the median age in the population were tested.

Model selection was based on the comparison of the objective function value (OFV) between nested models, precision in parameter estimates, and scientific plausibility. A difference in OFV of −3.84 or less is significant at the P value of 0.05 or less level, for 1 added degree of freedom. Performance of the RTTE models was evaluated by a visual predictive check (VPC), which is the Kaplan–Meier (KM) curve of the observed data plus a 95% prediction interval of simulated data from the RTTE model with 1,000 replicates. The observed censoring in the original data was used in the VPC. CIs for parameter estimates were derived using 1,000 bootstraps performed in PsN (version 3.5.3; Department of Pharmaceutical Biosciences, Uppsala University, Sweden).²³ 

The probability of requesting morphine during a 24-h interval for each consecutive event was visualized in a typical 80-yr-old female patient having surgery at 8 am with no dropout before 96 h postsurgery. The time of analgesic events for the typical subject was simulated 1,000 times using the final RTTE parameter estimates with a fixed set of covariates corresponding to this subject. Probabilities for each event to occur at a given time were calculated as follows²⁰ :

where F(t) is the cumulative distribution function giving the probability of having an event before time t and S(t) is the KM survivor function of the 1,000 simulations.

The cumulative hazard ratios for days 1 to 4 were calculated for the typical patient with average frailty as the integral (area under the curve) of hazard in time periods 0 to 24 h, 24 to 48 h, 48 to 72 h, and 72 to 96 h by using R (version 2.15.3; The R Foundation for Statistical Computing, Austria).²⁴ 

Hazard ratios were calculated to compare the effect of the log-linear covariates on the probability of an analgesic event. For a dichotomous covariate (day vs. night and sex), the hazard ratios represent the change in probability of event for a patient during the night compared with during the day. For continuous covariate (age), the hazard ratio quantified the change in probability of event for every year.

Data set preparation, data exploration, graphical analyses, modeling, and simulations were performed by using the software packages R (version 2.15.3; The R Foundation for Statistical Computing), NONMEM^® (version 7.2; ICON Development Solutions), Xpose (version 4.5.0; Uppsala University, Sweden), and PsN (version 3.5.3).

Of 109 eligible patients, 41 patients were excluded for receiving other analgesics than the hospital standard (acetaminophen plus morphine on request). Five patients were excluded because of either none or multiple surgeries. Descriptive statistics for the 63 included patients is available in table 1. The distribution of times of analgesic events as well as censoring time is seen in figure 1. Eight patients were censored in the interval 0 to 24 h; another 11 patients in the interval 24 to 48 h; 8 patients in the interval 48 to 72 h; 9 patients in the interval 72 to 96 h; and the remaining 27 patients at 96 h. A total of 302 analgesic events followed by a morphine dose were included in the analysis, ranging from 0 to 22 analgesic events among patients. Morphine was given as IV doses ranging from 2.5 to 10 mg, oral immediate-release doses ranging from 2.5 to 30 mg, and oral controlled-release doses ranging from 5 to 30 mg. No information was available in records on the basis for dose selection, and all analgesic events were treated identical for this illustration.

The fitted baseline hazard function of the three tested distribution models is shown in the left panels of figure 2. On the right are the VPCs that describe the adequacy of the three models’ ability to predict the observed KM survival for the first, second, third, and fourth analgesic events over time. The first plot thus shows that 89% of the patients had at least one analgesic event, the second plot shows that 80% had at least two analgesic events during the 96 h postsurgery, and so on. The interval around each observed KM curve is the 95% CI, which shows how well each model distribution described the observed data.

A Gompertz baseline model with decreasing hazard of analgesic events over time described the data significantly better than an exponential model with constant hazard (ΔOFV = −26.8). A Weibull distribution model had slightly lower OFV than the Gompertz model (ΔOFV = −3.7), but the Gompertz model better predicted the event numbers 3 and 4 (fig. 2) and was therefore chosen as the base model.

The covariate night (11 pm to 7 am) was statistically significant (P < 0.001) with a predicted hazard ratio of 32% relative to day (7 am to 11 pm), that is, the probability of analgesic events during the night was 32% of the probability during the day. Neither age (P = 0.057) nor sex (P = 0.056) was found to be significant covariates at the classical criterion for statistical significance of P value less than 0.05. Table 2 shows the parameter estimates of the final model. Figure 3 shows the VPCs of the final model with covariates for 1st to the 12th event. No tendency was found to indicate that frailty predicted censoring.

The probability of having one or more analgesic event(s) in a typical female patient of 80 yr age with surgery at 8 am is seen in figure 4. The probability of having one analgesic event was 80% within 24 h, 91% within 48 h, 95% within 72 h, and 97% within 96 h. The probability of having at least two analgesic events within 24 h was 55%, 31% for three analgesic events, and 18% for four analgesic events. The cumulative hazard given by the integral of hazard was 1.8 within 0 to 24 h, 1.4 within 24 to 48 h, 1.2 within 48 to 72 h, and 0.95 within 72 to 96 h for a typical 80-yr-old female patient having surgery at 8 am.

In the final model, the frailty was estimated to be 80% (expressed as coefficient of variation) in the hazard of analgesic events between patients. The frailty and dynamical changes in the probability over time and at night, as predicted by the RTTE model, are illustrated in figure 5. The figure shows the hazard versus time profile for five patients who were censored at 96 h and had estimated minimum, median, maximum, and interquartile (25th and 75th percentiles) estimates of individual frailty.

In this work, we have demonstrated the application of RTTE modeling for the study of timing of consecutive analgesic events in postoperative pain. It was seen that an RTTE model overall appropriately described the probabilities of consecutive analgesic events in time and thus gives a correct description of data. By using RTTE modeling, we described time-related changes in probability of analgesic events and found a general reduction in probability over time and a significant reduction in probability during night compared with day.

A main point of this demonstration was to present the possibility of studying consecutive analgesic events in postoperative pain. To study this, registration of the accurate timing of each consecutive analgesic event is needed on a subject-level basis. As demonstrated, this could be analyzed directly in the clinical setting but could equally well be applied to data from clinical trials as an alternative or a supplement to traditional analysis methods of analgesic consumption.

Repeated time to event analysis is related to TTE analysis as both build on mathematical methods that can appropriately analyze timing of discrete events, as analgesic events (requests) inherently are. Furthermore, RTTE and TTE do not rely on normal distribution or continuity of data and can handle censoring.^18,25  TTE modeling of first analgesic event has been used in analgesic efficacy trials^5,10  and has proven useful for studying covariates, accounting for informative dropouts, and studying pharmacokinetic–pharmacodynamic relations using nonlinear mixed-effects modeling.^18,25  The RTTE modeling shares the qualities needed for these types of analysis. However, RTTE analysis is not limited to study only a single analgesic event and thus does not share the same problem with censoring and loss of information as TTE of first analgesic event. This makes RTTE modeling better suited for longer studies and studies where multiple rescue medications may be given per individual. For studying health conditions, RTTE modeling has previously proven superior to a TTE model for detecting sex differences in functional disability.²¹ 

The RTTE model mathematically defines the probability of having analgesic events in time and thus allows for a continuous description of analgesic consumption after surgery. In this study, two commonly used parametric models (Gompertz and Weibull) with different shapes of time-changing probability were compared to a model with a constant probability (exponential distribution). The adequacy of the models was evaluated by using both simulation-based visual diagnostic (VPCs) and statistical fit. The significance of time dependency was tested by using RTTE modeling by comparing these models versus the exponential model with constant probability in time.¹⁸  The Gompertz model with a declining probability of analgesic events in time was found to describe data significantly better, signifying that the probability of analgesic events was not constant. The fits of the Weibull and Gompertz models were very similar when using the OFV criteria. However, using the VPC, a deviation to the trend was seen in event numbers 3 and 4 that was more pronounced with the Weibull than with the Gompertz model. Other more advanced models such as a surge function model might have given a better fit,¹⁷  but with an overall satisfactory fit of the Gompertz model, this was considered out of scope for this analysis.

The time dependency of observed consecutive analgesic events was well described by the RTTE model as evaluated on VPC (fig. 3). In this study, the probability of analgesic events was highest immediately after surgery and was declining with a rate corresponding to a 50% reduction on 3.3 days after surgery. Thus, RTTE modeling is a tool that allows for identification and quantification of time periods when the treatment need and probability of analgesic events is high. Also, an RTTE model may be used to predict the probability of analgesic events in time given a set of covariates. This was illustrated in figure 4 where the probabilities of analgesic events in time (F[t]) were predicted for a typical female subject 80 yr of age by using simulations. In the clinic, estimates of frailty and hazard profiles in the individual patients based on analgesic consumption history as illustrated in figure 5 could further inform the probability of consecutive analgesic events and identify patients who are at high risk or low risk of needing analgesics. Such predictions could be used to develop quantitative clinical rationales for adjusting doses, changing analgesia, or initiate adjuvant therapy based on covariates or estimates of individual frailty, where clinical decisions would otherwise be done empirically.

The hazard describes the instantaneous rate of events with a unit of events per time. In conditions where events are repeatable and incidental (i.e., events can be assumed to be independent of each other), the hazard and parameters of the hazard function can be interpreted as the number of events that occurs in time (count interpretation).²⁰  The cumulative hazard thus gives the probable number of events in a time interval. In this study, the baseline parameter of the Gompertz function described a hazard of 2.5 day⁻¹ for a typical female patient of 80 yr age. Such a patient would thus be expected to request 2.5 doses of morphine per day if she remained at the same hazard as immediately after surgery. However, the model also described individual variability and time-dependent hazard with a covariate describing reduction at night such as visualized in figure 4. The cumulative hazard over a time interval accounts for time dependency in hazard and gives the predicted number of analgesic events in the count interpretation. Thus, the predicted numbers of analgesic events were 1.8 on day 1, 1.4 on day 2, 1.2 on day 3, and 0.95 on day 4 for the typical 80-yr-old female subject.

Repeatable and incidental events are essential to the count interpretation of hazard. For example, it would not be meaningful to interpret a hazard of 2.5 day⁻¹ as the probable events per day in a TTE analysis, where only event can occur per subject. Censoring should also be considered, as informative dropout may affect the apparent hazard. Informative dropout occurs if the hazard and censoring time is related, for example, when rescue medication leads to dropout because of high pain intensity in a control group.²⁵  If unaccounted for the apparent hazard may drop over the course of a study as a result of high-risk patients dropping out. Tendencies between individual frailty or other predictors and censoring time should therefore be explored. In case of informative dropout, a dropout model can be implemented to obtain unbiased prediction of events in time.²⁵  In this study, we did not find tendencies to informative dropout when evaluating individual frailty versus censoring time.

When investigating new approaches to pain management, effects on different parameters of the hazard function may give different information about the type of intervention. Thus some interventions may reduce the baseline hazard proportionally (e.g., concentration-related symptomatic treatment with adjuvant analgesia), whereas others may affect the shape parameter and thus accelerate the decline in probability of events (e.g., interventions that improve recovery).

The coefficient of variation (80%) for the random effect on hazard described how much the probable number of events varied between individuals.¹⁷  This variability is termed frailty in statistical literature to denote that some subjects are more susceptible to events than others²⁰  and is sometimes referred to as interindividual variability in RTTE modeling.¹⁵  The possibility to describe discrete data as probabilities with a continuous distribution of hazards between individuals is the key to appropriate description of analgesic consumption data and is not readily done with traditional summary statistics.^8,9  We used a log-normal distribution for frailty, which is commonly used for RTTE modeling with nonlinear mixed effects such as drug concentration–effect relations.^15,17  Other distributions such as the gamma distribution are also frequently used to describe frailty in RTTE modeling²⁰ ; however, this random distribution is not available in NONMEM^®, version 7.2.²² 

One major challenge in analgesic consumption analysis has been to correlate and account for time-varying factors, most notably pain intensity and side effects such as nausea and vomiting.^12,13  This is particularly relevant when comparing analgesics, as a reduction in analgesic consumption could have several causes that may fluctuate in time. It is thus of essence that the analysis method allows these factors to be analyzed as time-varying covariates.

In this study, a time-varying covariate was implemented in the RTTE model as the probability of analgesic events was found to be significantly reduced during the hours of the night (11 pm to 7 am) with a reduction in probability of 32% compared with day. The exact cause of this reduction is unknown, but one causal factor could be sleeping patterns.

In any case, failure to take into account the reduction in probability at night would have underestimated the probability of having analgesic events during the hours of the day. Age and sex were also implemented as covariates but were not found significant at a P value less than 0.05 level. However, the ability to study multiple time-dependent factors is an important feature of RTTE modeling that could be used to study the correlation between pain intensity, nausea and vomiting, and exposure–response relations of opioid-sparing drugs in future trials.

High-quality clinical research on opioid-sparing interventions in postoperative pain management has been conducted.^3,7  However, the overall value of opioid-sparing interventions remains unclear. Opioid-sparing efficacy is often claimed for drugs based on a reduction in mean 24-h analgesic consumption from meta-analyses of small clinical trials.^7,14  This is despite knowledge that analgesic consumption is not normally distributed, as demonstrated by Moore et al.⁹  who studied the 0- to 24-h, 0- to 48-h, and 24- to 48-h analgesic consumption of fentanyl in a phase 2 and phase 3 drug development program with 917 patients undergoing lower-segment cesarean section, lower abdominal surgery, or hip arthroscopy.

The heterogeneity of trials and inconsistency in use of statistical tests also constitute a major problem, as large discrepancies in apparent effects have been found between trials.^5,7  Traditional statistical tests such as ANOVA, t test, Mann–Whitney test, and Kruskal–Wallis test all involve violation of statistical assumptions when used to study analgesic consumption data.⁸  Furthermore, the yes/no type of information that is derived by such tests may be difficult to translate into other clinical settings. We hope that this illustration of RTTE modeling as an analysis tool will facilitate better knowledge from postoperative pain trials and in particular that it can be used to bridge results from clinical trials to make evidence-based rationales for clinical pain management.

This work was done using retrospective data from an actual clinical setting, and some assumptions were made in this analysis. First, it was assumed that all opioid doses registered in the medical journal were given following a request by the patient as per hospital guideline and that timing was accurately registered in the medical records. Second, analgesic events were treated identically and it was assumed that the choice of dosage form and dose was uninformative. However, none of these assumptions are regarded as influencing the possibility to show how RTTE can be applied using these types of data. Finally, the sample size used in this analysis was intended to illustrate the RTTE approach and should not be considered a model for future studies.

In conclusion, appropriate analysis is essential if unbiased and correct conclusions are to be drawn from studies measuring analgesic consumption. In this work, we show that RTTE modeling of the timing of consecutive analgesic events is an appropriate way of handling these types of trial data. The approach allows for implementation of a model-based analysis that can be used to quantify the time changes in the probability of analgesic events in postoperative pain management. Not only are the data correctly described, but also the method holds the potential to quantify multiple time-varying variables and dynamical relations, such as exposure–response of opioid-sparing adjuvant analgesics. We consider the analysis of analgesic consumption by RTTE modeling a valuable tool and encourage its use in postoperative pain research.

This work was funded by the Department of Drug Design and Pharmacology, University of Copenhagen, Copenhagen, Denmark.

The authors declare no competing interests.