Finding the right hazard function for time-to-event modeling: A tutorial and Shiny application.

Parametric time-to-event analysis is an important pharmacometric method to predict the probability of an event up until a certain time as a function of covariates and/or drug exposure. Modeling is performed at the level of the hazard function describing the instantaneous rate of an event occurring at that timepoint. We give an overview of the parametric time-to-event analysis starting with graphical exploration by Kaplan-Meier plotting for the event data including censoring and nonparametric hazard estimators such as the kernel-based visual hazard comparison for the underlying hazard. The most common hazard functions including the exponential, Gompertz, Weibull, log-normal, log-logistic, and circadian functions are described in detail. A Shiny application was developed to graphically guide the modeler which of the most common hazard functions presents a similar shape compared to the data in order to guide which hazard functions to test in the parametric time-to-event analysis. For the chosen hazard function(s), the Shiny application can additionally be used to explore corresponding parameter values to inform on suitable initial estimates for parametric modeling as well as on possible covariate or treatment relationships to certain parameters. Moreover, it can be used for the dissemination of results as well as communication, training, and workshops on time-to-event analysis. By guiding the modeler on which functions and what parameter values to test and compare as well as to assist in dissemination, the Shiny application developed here greatly supports the modeler in complicated parametric time-to-event modeling.

INTRODUCTION

Time‐to‐event analysis is an important tool in the pharmacometric toolbox to analyze event‐type data and predict the probability of an event up until at a certain time possibly accounting for covariates and/or treatment metrics. Event‐type data are one of the closest representations of clinical outcomes experienced by patients and therefore relevant and informative for the main objective of drug development. An event is any defined and delimited occurrence or experience at a specific timepoint, for example, a seizure, relapse or exacerbation, or death. Time‐to‐event analysis not only considers the occurrence of an event but also its timing. Time‐to‐event analysis is therefore favored over other analysis techniques such as logistic regression because of its incorporation of the event time as outcome. In addition, time‐to‐event analysis handles censoring of data appropriately. , Censoring occurs when the event is not observed because it is outside of the observation time (left‐ or right‐censoring) or between two observations (interval censoring). In clinical trials, this can happen because of end of trial or follow‐up, lost to follow‐up, or withdrawal from the study. Individuals are commonly right‐censored, that is, the event time is larger than the censor time. Left‐censoring, when the event time is less than or equal to the censor time, is less common, whereas interval‐censoring, when the event time is within a known time interval between two observations, occurs more often. It is important to be aware of competing risks for the event under analysis when considering censoring, , but this is outside of the scope of this tutorial.

Because of the special nature of the data, time‐to‐event analysis or survival analysis relies on two main functions. The two main functions in time‐to‐event analysis are the survival (S) function, quantifying an individual's probability of reaching a timepoint without an event, and the hazard (h) function, quantifying the instantaneous rate of the event happening at that timepoint. Although the survival function is bounded to [0,1], the hazard function's lower boundary is 0, but it has no upper boundary ([0,∞>). The survival function relates to the hazard function as shown in Equation (1), exponentiating the negative cumulative hazard function. The product of the hazard function and the survival function yields the probability density function of the event times. Analysis of event time data can be done empirically, semiparametrically, and parametrically. The simplest approach is empirical and descriptive in nature. The nonparametric Kaplan–Meier estimator is a common example. The impact of predictive covariates such as demographics or treatment group on the event time as outcome has to be assessed by analysis of the corresponding subset of the data. Semiparametric analyses such as the Cox proportional hazard model do not define the baseline hazard function and only define the event time in terms of covariates and corresponding estimated coefficient/parameter, from which a hazard ratio can be calculated. Although semiparametric time‐to‐event analysis does not require assuming a hazard function, other assumptions have to be made on, for example, proportional hazards, proportional odds, or accelerated failure time where the hazards, odds, or failure times are assumed to be related to each other by a constant factor such as the hazard ratio. In addition, semiparametric analyses are limited to the observed timepoints in the data and will have output that corresponds to those timepoints, such as the stepwise Kaplan–Meier curve. In contrast, a fully parametric approach defines the underlying baseline hazard function in addition to the covariate coefficient and the baseline hazard function and can be used for simulations in contrast to semiparametric approaches. Therefore, it is a more complete method to characterize the survival and hazard functions. A comprehensive introduction to time‐to‐event or survival analysis can be found elsewhere as well as a tutorial to its application on pharmacometrics. To build on that tutorial, the scope of this work is describing the most commonly used hazard functions in parametric time‐to‐event pharmacometric analyses and supporting model development using those hazard functions.

The hazard function characterizing the outcome distribution is at the core of parametric time‐to‐event analysis. Especially when analytical solutions are not feasible, for example, with time‐varying components, modeling takes place at the hazard function level. Also, coefficients quantifying the effect of covariates and treatment are modeled at the hazard level. It is therefore convenient to graphically explore the observed data and model diagnostics at the hazard level. Recent methods have been proposed to characterize the shape of the hazard profile to assist in time‐to‐event model development. These nonparametric methods use estimators for the hazard based on histogram or smoothing functions. , , , The hazard‐based visual predictive check (VPC) estimates the hazard over time based on the observed data nonparametrically using a binned hazard estimator or a kernel method with either local or global bandwidth. The nonparametric hazard estimates reflecting the observed data are subsequently overlaid with (parametric) simulations from the model to diagnose the model fit. The kernel‐based visual hazard comparison (kbVHC) in contrast is a simulation‐free method that uses the nonparametric kernel estimator of the hazard to calculate a 95% confidence interval of the hazard over time based on the observed data. The parametric model‐predicted hazard over time can subsequently be overlaid with the confidence interval reflecting the observed data. Both methods have been proposed as model diagnostics; however, they can also be used to explore the (raw) data graphically as an important first step in pharmacometric model development. Here, we use the kbVHC because its confidence interval reflects the variability in the observed data, but the hazard‐based VPC can be used interchangeably.

Theoretically, any user‐defined function can be used to describe the hazard, and hazard functions can range from a simple, one‐parameter constant to multiparameter, nonmonotonic, or periodic functions. The choice of an appropriate hazard function can be challenging. Also, more complex hazard functions can be sensitive to the initial estimates for the parameters provided to the estimation algorithm. In addition, as time‐to‐event analyses are very often collaborations with clinicians or other scientists with less experience in modeling, it would be helpful to visualize hazard functions to make the abstract concept more graspable. The objective of this tutorial is therefore to support time‐to‐event model building by developing a Shiny application to find the right hazard function(s), including their appropriate initial estimates, to explore parameter sensitivity to inform on, for example, covariate or treatment relationships, and to assist in dissemination.

TIME‐TO‐EVENT ANALYSIS WORKFLOW

The workflow of a typical time‐to‐event analysis is shown in Figure 1. Generally, the first step in data analysis is exploration of the data either graphically or numerically. One way to graphically explore the observed individual event times is by showing the follow‐up time including occurring event(s) per subject over time (Figure 2a). In addition, the survival probability over time in the population including censoring can be explored by a Kaplan–Meier plot (Figure 2b). A more informative way to graphically explore the observed data regarding the underlying hazard is by using a nonparametric hazard estimation methods based on the histogram or smoothing functions mentioned previously (Figure 2c). , Similar to, for example, assessing one‐ or two‐compartmental behavior in pharmacokinetic modeling by log‐linear graphics, the nonparametric hazard estimation informs the modeler of the shape of the hazard over time and which parametric functions to test. Thus, a graphical exploration of the hazard can be taken forward to the choice of the hazard functions to test which shape corresponds to the hazard estimation most. In the case of Figure 2c, the hazard shows an increasing and subsequently decreasing pattern over time, so hazard functions capable of that dynamic, such as log‐normal or log‐logistic, should be selected to be included in model development. The most common hazard functions are discussed in more detail in the “Hazard Functions for Parametric Time‐To‐Event Analysis” section. To start model development for the selected hazard functions, the initial estimates for the hazard function parameters need to be chosen, which can be challenging for abstract and less‐intuitive hazard functions. Here, a Shiny application was developed that can be used to that aim to explore which parameter values of, for example, the log‐normal or log‐logistic functions lead to a shape similar to that of the nonparametric hazard estimation (Figure 2c). Similarly, the sensitivity of the hazard over time to certain parameters can be explored. This can inform the modeler on the modeling strategy, which parameter, for example, might be best to start testing covariates or treatment effects on, especially with a nonparametric hazard estimation of the corresponding subset of the data (e.g., male/female, treatment/placebo). This is discussed in more detail in the “Exploration of Hazard Function Parameter Values” section. With the right hazard function(s) and corresponding initial estimates to test, model development using designated nonlinear mixed‐effects modeling software can be performed in a well‐informed manner. Finally, the results of the analysis will be communicated and disseminated, for which the developed Shiny application can be used as well (“Dissemination Using the Shiny Application” section).

FIGURE 1

Time‐to‐event analysis workflow. White box represents the input data, green boxes represent actions that can be supported by the proposed Shiny application, and blue boxes represent actions by other algorithms (e.g., nonparametric hazard estimators, nonlinear mixed‐effects modeling software)

FIGURE 2

Graphical exploration of a typical time‐to‐event data set. The data set was simulated with a log‐logistic hazard function (λ = 0.002, γ = 1.1) over 365 days for 1000 individuals (IDs) with 10% random right‐censoring. (a) Follow‐up time (horizontal line), event (circle symbol), and censoring (plus symbol) for 100 IDs randomly drawn from the simulated data and ordered on event time; (b) Kaplan–Meier plot of the survival probability including censoring (plus symbol) based on the full simulated data set; and (c) kernel‐based visual hazard comparison showing the nonparametric kernel‐based hazard estimate (dashed line) and its 95% confidence interval (shaded area) based on the full simulated data set

HAZARD FUNCTIONS FOR PARAMETRIC TIME‐TO‐EVENT ANALYSIS

The hazard functions are described in detail next and shown in Table 1 and Figure 3. The more complicated hazard functions commonly contain a shape parameter (e.g., γ) modifying the curve of the function and a scale parameter (e.g., λ).

TABLE 1

Hazard functions for the most commonly used models in time‐to‐event analysis

Model	Hazard function
Exponential	ht=λ
Gompertz	ht=λ·e^γ·t
Weibull	ht=λ·γ·λ·t^γ−1 Alternatively: ht=λ·γ·t^γ−1
Log‐normal	ht=σ·t·2π^−1·e^−12·logt−μσ^21−Φlogt−μσ
Log‐logistic	ht=λ·γ·λ·t^γ−11+λ·t^γ Alternatively: ht=λ·γ·t^γ−11+λ·t^γ
Generalized gamma	St=1−Γλ^−2e^−βt^λσλ^−2ifλ>0Γλ^−2e^−βt^λσλ^−2ifλ<0 ft=λσtΓλ^−2λ^−2e^−βt^λσ^λ^−2e^−λ^−2e^−βt^λσ ht=ftSt
Circadian	ht=λ·1+amp·Sin2πperiod·t+phase

FIGURE 3

Hazard over time for the exponential model (solid line: λ = 0.1; dotted line: λ = 0.05; dashed line: λ = 0.01), the Gompertz model (λ = 0.02, solid line: γ = −0.01; dotted line: γ = −0.001; dashed line: γ = 0.002), the Weibull model (λ = 0.02, solid line: γ = 0.4; dotted line: γ = 0.8; dashed line: γ = 1.2; alternative parameterization in gray), the log‐normal model (μ = 5, solid line: σ = 1.5; dotted line: σ = 1; dashed line: σ = 0.75), the log‐logistic model (λ = 0.02, solid line: γ = 0.9; dotted line: γ = 1.1; dashed line: γ = 1.3; alternative parameterization in gray), the generalized gamma model (solid line: μ = 1.1, Q = 1.5, σ = 1.2; dotted line: μ = 1.2, Q = 0.4, σ = 0.5; dashed line: μ = 0.9, Q = 1.2, σ = 1; dot‐dashed line: μ = 0.3, Q = 0.5, σ = 1; Prentice parameterization), and the circadian model (λ = 0.02, period = 360, solid line: amplitude = 0.1, phase = 0; dotted line: amplitude = 0.2, phase = 180; dashed line: amplitude = 0.5, phase = 0)

Exponential model

The constant hazard function is the simplest hazard function with no change in the hazard over time (Table 1, Figure 3). The survival function is described by an exponential function of minus λ times time, and the probability density function of the events time, which equals the product of the hazard function and the survival function, shows an exponential distribution. Because of this, the constant hazard function is also referred to as the exponential model in time‐to‐event analysis. The constant hazard function is more restrictive than the semiparametric Cox proportional hazard function because it assumes a constant (baseline) hazard function, whereas the Cox model only assumes the proportionality of hazards when (undefined) hazard functions do not need to be constant as long as they are proportional. It has, for example, been applied to characterize the time to positivity in a tuberculosis liquid culture assay, time to acute idiopathic pulmonary fibrosis exacerbations, time to recrudescent visceral leishmaniasis infection, and overall survival of patients with acute myeloid leukemia.

Gompertz model

The Gompertz model is a two‐parameter proportional hazard function (Table 1, Figure 3). It is similar to the Cox proportional hazard model but has a defined baseline hazard function with a shape parameter γ following the Gompertz distribution. When γ < 0, the hazard decreases over time, whereas y > 0 characterizes an increasing hazard with time. The Gompertz model reduces to the exponential model for γ = 0. Its application in pharmacometrics include the characterization of the time to recurrent venous thromboembolism, acute deep vein thrombosis, and pulmonary embolism ; the time to analgesic remedication after surgery, , , which can be expected to decrease with time after surgery; and the time to recurrent Clostridium difficile infection.

Weibull model

The Weibull model is a two‐parameter model with flexibility to characterize both increasing and decreasing hazard over time, similar to the Gompertz model, but for the shape parameter >1 or <1, respectively. , The Weibull model is a more recent addition to time‐to‐event analysis compared with the Gompertz model and shows more flexibility in, for example, characterizing an increasing hazard of which the slope decreases with time. If the shape parameter does not vary with covariates, both the accelerated failure time and proportional hazard assumptions hold, which is unique to this model. , This makes the Weibull model attractive for intuitive interpretation by nonmodelers regarding the hazard ratio between populations and a treatment or covariate effect on the survival time.

The Weibull model is commonly found in pharmacometric analyses in, for example, time to relapse in multiple sclerosis, time to epileptic seizure, time to onset of drug toxicity in preclinical experiments with both left‐ and right‐censoring, and overall survival in patients with metastatic colorectal cancer, and therefore it is also used in more theoretical exercises. , , An alternative parameterization (see Table 1, Figure 3) has been observed in the tutorial by Holford and the handbook by Kleinbaum and Klein as well as in applications to time to adverse drug events, time to anemia, and mortality time in patients with amyotrophic lateral sclerosis. The two functions differ by the factor λ and are thus equal at λ = 1 (and λ = 0). They collapse for λalternative = λγ. For λ > 1, parameterization 1 > parameterization 2 for t > 0, whereas 0 < λ < 1 results in parameterization 1 < parameterization 2 for t > 0. The alternative parameterization 2 seems to be more unstable and sensitive to initial estimates in our experience, running into numerical difficulties with the integration routine. With respect to reproducibility, it is recommended to include the exact formula of the hazard function in the methods in addition to its name.

Log‐normal model

The log‐normal distribution is an accelerated failure time model, but not a proportional odds or propotional hazard model (Table 1, Figure 3). The hazard increases to a maximum and then decreases over time (nonmonotically), resulting in a log‐normal distribution of the event time. Although it is a two‐parameter model (𝜇 and 𝜎 as the mean and standard deviation) similar to the Gompertz and Weibull models, its parameterization is more complex, including the standard normal cumulative distribution function Φ. This complex parameterization makes an intuitive choice of initial estimate challenging. The log‐normal model has, for example, been applied to time to undetectable levels of hepatitis C virus and overall survival and progression‐free survival of patients with metastatic gastrointestinal stromal tumors.

Log‐logistic model

The log‐logistic distribution allows the hazard to increase for a shape parameter γ > 1, while it decreases over time for shape parameter γ ≤ 1 (Table 1, Figure 3). If the accelerated failure time assumption holds, the log‐logistic distribution is a proportional odds model, that is, the odds ratio is constant over time. The event time approximates a log‐logistic distribution. Similar to the Weibull model, the log‐logistic model has an alternative parameterization. An advantage of the log‐logistic model over the log‐normal model is that is can be solved analytically, which decreases run times. However, finding initial estimates remains challenging. The log‐logistic model has been applied in characterizing the time to pain response in preclinical experiments where it showed better performance compared with the Weibull model in run times and predictive performance as well as in modeling dropout in a metastatic colorectal cancer clinical trial and progression‐free survival in anaplastic lymphoma kinase‐positive non‐small cell lung cancer.

Generalized gamma model

The generalized gamma model with three parameters (scale σ, shape λ, and location β) allows for even more flexibility than the previous models (Table 1). The flexibility of the generalized gamma model allows its hazard function to take the most common distributions: increasing over time, decreasing over time, decreasing then increasing over time (U shape or bathtub), and increasing then decreasing over time (inverted U shape). The location parameter accelerates the function with time, whereas the scale and shape parameter determine if the distribution increases, decreases, or has an (inverted) U shape over time as demonstrated by Cox et al. The generalized gamma model can collapse into simpler models, such as the Weibull (for λ = 1), log‐normal (for λ = 0), or exponential (for λ = σ = 1) model. , The generalized gamma model was reparameterized by Prentice to facilitate modeling that is incorporated in statistical software including R. , It has been applied to characterize the time to unassisted breathing or death in patients who are critically ill. It is, however, more commonly used as a tool to find the most appropriate time‐to‐event model by checking the Weibull or log‐normal model assumptions than used as a final model itself because the generalized gamma model can reduce to the Weibull and log‐normal models. , One reason is its complexity, including the standardized cumulative gamma distribution function (Γ, mean and variance γ > 0).

Circadian or periodical model

A more mechanistic characterization of the instantaneous rate of an event over time can be captured by a circadian or more general periodic function, where the hazard changes with the period of the day, month, or year (Table 1, Figure 3). Especially for clinical trials that cover an extended period of time and can be influenced by seasonality, such models can become relevant. The circadian model consists of three additional parameters that cover the amplitude, period, and phase of the periodic relationship in the hazard and can be used as multiplicative factor for any baseline hazard functions (shown in Table 1 for the exponential model).

EXPLORATION OF HAZARD FUNCTION PARAMETER VALUES

Finding the right hazard functions to test is a first step in time‐to‐event model development (see Figure 1), but model algorithms require initial estimates. The more complex hazard functions are sensitive to the initial estimates the modeler starts with and can easily run into numerical issues. This was even observed for different parameterizations of, for example, the Weibull model. The descriptions of the aforementioned models give a first idea on the constraints of the parameter values given a certain shape (e.g., λ > 1 for an increasing hazard with time in the Weibull model), but not on the order of magnitude of the parameters. It is beneficial to modelers to familiarize themselves with the parameter space in relation to the exact shape of the hazard over time of the observed data. In addition, it is important to be aware of the sensitivity of the hazard over time to the different parameters. Exploring the parameter space can inform on which covariate relationships should be tested, especially when a graphical exploration based on a nonparametric hazard estimation was performed on subsets of the data for the covariate(s) of interest. The same holds for testing the treatment effect. If, for example, the peak of the hazard over time in Figure 2c occurs later in the treatment group compared with the placebo group, that effect might be captured best by modeling the treatment at the γ parameter, whereas if the peak decreases, it might be better captured on the λ parameter.

To explore the hazard function parameter values and their influence on the hazard over time, we have developed a Shiny application (available at https://pqp‐uu.shinyapps.io/HazardFunctionsInParametricTTE/ and https://github.com/rcvanwijk/HazardFunctionsInParametricTTE; Figure 4 and Supplementary Material) using the Shiny package (Version 1.6.0) in R (Version 4.0.4) through the RStudio (Version 1.4.1106) interface. The application simulates the six most common hazard functions (exponential, Gompertz, Weibull, log‐normal, log‐logistic, circadian; Table 1) in a reactive manner, meaning it responds automatically to changes the user makes in the input parameters. The impact of covariates/treatment on the hazard over time can be explored by an additional set of inputs. A covariate/treatment can affect the hazard function proportionally. For each function, a dedicated coefficient input is available in the application. Under the proportional hazard assumption, a hazard ratio can be calculated by exponentiating the coefficient. Alternatively, a covariate/treatment can affect the hazard function's individual parameters for which the application has a relative input per parameter. Both slider and numerical Shiny input widgets are built into the application for function parameters and covariate/treatment effect as well as for the simulation time to give the user full control. For each function, a figure based on the reactive simulation is assigned to a reactive function that is rendered to display in the application and made available in tiff format through a download handler. The Shiny application was code reviewed and tested independently by two researchers.

FIGURE 4

Shiny application with (a) slider panel for the hazard function parameters, (b) slider panel for covariate/treatment effect on the hazard function parameters, and (c) the corresponding hazard over time profiles, here illustrated for log‐logistic and log‐normal functions

Figure 4 shows an overview of the Shiny application. Sliders or numerical input for the hazard function parameters (Figure 4a) and the covariate/treatment effect on those parameters (Figure 4b) can be used to simulate the hazard over time profile in the graphical panel (Figure 4c). With the nonparametric hazard estimation based on the observed events, the right functions showing a similar shape can be readily selected. Based on Figure 2c, either the log‐normal or log‐logistic model would be best to describe the data. In addition, the input parameters can be updated to improve the shape toward the observed data. These values can then be used as initial estimates for the nonlinear mixed‐effects modeling software. For the log‐normal and log‐logistic models, the parameter values were adjusted (Figure 4a) to approximate the nonparametric hazard shape (compare Figure 4c with Figure 2c), and these values can now be used as initial estimates into the parametric time‐to‐event algorithm. Changing the parameters using the sliders will also illustrate the sensitivity of the shape of the function to different parameters. In addition, the covariate/treatment effect sliders can even be used to explore which parameter should be tested for covariate and/or treatment effects. Overall, a better understanding of the meaning and impact of the parameters will result from the exploration of the hazard functions using the Shiny application developed here, which will improve the following model development steps.

DISSEMINATION USING THE SHINY APPLICATION

Time‐to‐event analysis in general and the hazard function specifically are abstract concepts that are difficult to intuitively understand by nonexperts. Dissemination of the time‐to‐event projects can be supported by informative graphics. One example is the Kaplan–Meier VPC, which overlays the time‐to‐event observations with a prediction interval based on a high number (e.g., 500 or 1000) of simulated time‐to‐event data sets using the final model. The advantage of the Kaplan–Meier VPC is that it shows the model performance with a graphic commonly used and understood. However, the simulation‐based technique is not without limitations, including a longer run time and less sensitivity to more complex hazard functions. In principle it is better to communicate at the level where modeling decisions take place on, for example, treatment and covariate effects. When developing covariate models in pharmacokinetics, the covariate relationship to the pharmacokinetic parameter such as clearance is presented in addition to the concentration‐time profile. The same holds for the hazard function, albeit complex. To aid with that complexity, the Shiny application developed here can be used. In presentations or workshops, the different hazard functions can be visualized as well as the meaning and impact of the parameters on their shape. The Shiny application also includes the possibility to export the created figures that can be used in further communication. In addition to specific project dissemination, the Shiny application can be used in educational programs teaching time‐to‐event analysis.

DISCUSSION AND CONCLUSION

We have developed a Shiny application to find the right hazard function(s) to test in time‐to‐event modeling, explore the influence of the parameter values on the hazard over time to inform on initial estimates and possible covariate or treatment relationships to test, and support the dissemination of project results.

Time‐to‐event modeling is an important method in the toolbox of the pharmacometrician. It is an excellent tool to characterize not only the occurrence of a clinical (or preclinical) event but also its time as well as the impact of censored data. That quantitative relationship is essential for interpretation of treatment and covariate effects and translation to unknown phases in drug development, clinical situations, or populations. As time‐to‐event modeling and quantification of treatment and covariate effects takes place at the hazard level, the hazard functions at the pharmacometrician's disposal must be well understood. Combining recently developed methods to graphically explore the time‐to‐event data by hazard estimators , with knowledge of the theoretical hazard functions is essential for time‐to‐event model development. The Shiny application developed here assists in that understanding.

The hazard function is an essential element in parametric time‐to‐event modeling, but other elements are important to address as well. Dropout that is nonrandom between, for example, the treatment and the control group, and therefore informative on the treatment effect, should be considered. The Shiny application developed here is limited to visualize informative dropout, which can be detected more readily in, for example, Kaplan–Meier curves displaying dropout as right‐censoring. Another element that complicates time‐to‐event modeling is time‐varying elements such as covariates or treatment. Although the Shiny application developed here could visualize periodic hazard functions through, for example, the circadian function, it was not designed to detect time‐varying covariate relationships in parametric time‐to‐event analysis. A combination of tools is at the disposal of the pharmacometrician to elucidate the more complex relationships in drug development to which we add this Shiny application.

In conclusion, we described in detail the most common hazard functions including mathematical and graphical representation and developed a Shiny application to support time‐to‐event modelers to start model development in a well‐informed manner as well as offer graphics to aid in the dissemination of results.

CONFLICT OF INTEREST

The authors declared no competing interests for this work.

Appendix S1

Click here for additional data file.