Voeltz et al.

BMC Medical Research

BMC Medical Research Methodology (2024) 24:48
[Link] Methodology

RESEARCH Open Access

A parametric additive hazard model

for time‑to‑event analysis
Dina Voeltz1,2*† , Annika Hoyer1† , Amelie Forkel2, Anke Schwandt3   and Oliver Kuß4   

Abstract
Background In recent years, the use of non- and semi-parametric models which estimate hazard ratios for analysing
time-to-event outcomes is continuously criticized in terms of interpretation, technical implementation, and flexibility.
Hazard ratios in particular are critically discussed for their misleading interpretation as relative risks and their non-
collapsibility. Additive hazard models do not have these drawbacks but are rarely used because they assume a non-
or semi-parametric additive hazard which renders computation and interpretation complicated.
Methods As a remedy, we propose a new parametric additive hazard model that allows results to be reported
on the original time rather than on the hazard scale. Being an essentially parametric model, survival, hazard and prob-
ability density functions are directly available. Parameter estimation is straightforward by maximizing the log-likeli-
hood function.
Results Applying the model to different parametric distributions in a simulation study and in an exemplary applica-
tion using data from a study investigating medical care to lung cancer patients, we show that the approach works
well in practice.
Conclusions Our proposed parametric additive hazard model can serve as a powerful tool to analyze time-to-event
outcomes due to its simple interpretation, flexibility and facilitated parameter estimation.
Keywords Additive hazard, Parametric modeling, Survival analysis, Time-to-event model

Background criticized in recent years. Most relevant points of con-

Regression models for time-to-event outcomes are pref- cern were that the hazard ratio (i) is interpretable only
erentially fitted by the proportional hazard model [9]. when mistaken as a relative risk [25], (ii) has a built-in
This is surprising because the hazard ratio, which is the selection or left truncation bias even in randomized trials
generic effect estimate of that model, has been repeatedly [2, 14, 24], and (iii) is non-collapsible [12], meaning that
adjusting for a covariate that is associated with the event
†
will in general change the hazard ratio, even if this covar-
Dina Voeltz and Annika Hoyer contributed equally to this work.
iate is not associated with the exposure. To conquer these
*Correspondence: shortcomings, parametric survival models have been rec-
Dina Voeltz
[Link]@[Link] ommended as they are simpler, more informative, more
1
Biostatistics and Medical Biometry, Medical School OWL, Bielefeld robust, and have hazard, survival and probability density
University, Universitätsstr. 25, Bielefeld 33615, Germany function directly available [18, 21].
2
Department of Statistics, Ludwig-Maximilians-University München,
München, Germany A hazard-based model that does not suffer from the
3
Nuvisan GmbH, Neu‑Ulm, Germany above mentioned problems is the additive hazard model
4
Institute for Biometrics and Epidemiology, German Diabetes Center, introduced by Aalen [1]. Results from additive hazard
Leibniz Center for Diabetes Research at Heinrich Heine University
Düsseldorf, Düsseldorf, Germany models can be translated to a relative survival scale which
is, for example, routinely used in cancer epidemiology

© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the
original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or
other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this
licence, visit [Link] The Creative Commons Public Domain Dedication waiver ([Link]
[Link]/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 2 of 11

when a cohort of cancer patients is compared to the gen- regression coefficients that measure the additive impact
eral population in terms of expected survival [11]. More- of covariates.
over, recent works showed that parameters of the additive Using the well-known relations between hazard, prob-
hazard model do not suffer from non-collapsibility [19], ability density function (pdf ) (f(t)) and survival function
at least when defined in continuous time [23]. How- S(t), the corresponding pdf of model (1), fx (t), can be
ever, the Aalen model is rarely used in applied research expressed in terms of the baseline pdf f0 (t), the baseline
because it assumes a non-parametric hazard as well as survival function S0,θ (t), and the covariates by
time-dependent covariates. Of course, these properties
give large flexibility in modeling, but also complicate f0 (t) + xβS0,θ (t)
fx (t) = . (2)
parameter estimation considerably. In addition, Bradburn exp(txβ)
et al. [4] state that “[...] the model coefficients are not
The complete derivation of the additive hazard
easy to understand, and as they change repeatedly over
model equation in terms of the pdf can be found in
time, can offer no single quantifiable effect size”. As a par-
Additional file 1.
tial solution for that problem, Lin et al. [17] proposed an
The corresponding survival function Sx (t) is given by
additive hazard model for (constant,) time-independent
covariates, but still use a semi-parametric hazard which S0,θ (t)
may require lengthy computations for some statistical Sx (t) = = 1 − Fx (t), (3)
exp(txi β)
standard software [22].
As a remedy, in this article we combine the additive with Fx (t) denoting the cumulative distribution function
hazard model with time-independent (, i.e., constant) (cdf ).
covariates and a parametric assumption for the base- Using an additive hazard regression model allows to
line hazard. This will result in a number of advantages in estimate relative survival instead of hazard ratios [27].
terms of interpretation, possible model extensions, and Contrarily to hazard ratios, which quantify the aver-
also enables parameter estimation using every software age or constant effect of covariates on the hazard func-
that allows maximizing a hand-coded likelihood function. tion, relative survival measures the cumulative effect of
In the Methods section, we proceed by introducing the covariates on the relative survival probability. Within
formal notation and deriving the model’s equations. Then, this context and in case x = 0 represents the absence
we report the settings of our simulation study. We illus- of disease, the relative survival probability can be
trate the model using a study which investigated provi- interpreted as the observed survival probability of the
sion of medical care to lung cancer patients in the eastern population studied, divided by the expected survival
part of Germany [3] and present all results in Results sec- probability if the population was free of the disease of
tion. Finally, we conclude with a discussion (Discussion interest. Using (3), we can confirm the relative survival
section). interpretation of our additive hazard model:
S0,θ (t)
Methods Sx (t) exp(txβ) 1
(4)
= = .
Parametric additive hazard model S0,θ (t) S0,θ (t) exp(txβ)
We start with the assumption that the hazard hx (t) with
covariate vector x at time t can be expressed as It should be noted that this relative survival interpre-
tation is independent of the baseline distribution.
hx (t) = h0,θ (t) + xβ, (1) The likelihood function for the model can be derived
by accounting for the fact that observations with an
with a parametric baseline hazard function h0,θ (t) which
event contribute the logarithm of the pdf, and censored
is independent of the covariates. The parameter θ denotes
observations the logarithm of the survival function [7].
the distribution parameters, which differ in terms of
The contribution of a single observation i ( i = 1, ..., N )
number and commonly used notation depending on the
with covariate vector xi and observation time ti to the
choice of the baseline distribution. The parameters β are
log-likelihood function ℓi is therefore

ℓi = (1 − δi ) log(fxi (ti )) + δi log(Sxi (ti ))

f0 (ti ) + xi βS0,θ (t)(ti ) S0,θ (t)(ti )
= (1 − δi ) log + δi log
exp(ti xi β) exp(ti xi β)
= (1 − δi )(log(f0 (ti ) + xi βS0,θ (t)(ti )) − ti xi β) + δi (log(S0,θ (t)(ti )) − ti xi β).
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 3 of 11

δi is the censoring indicator with δi = 1 if an obser- Accordingly, we focused on simulating a single binary
vation is censored, and δi = 0 if an event has been covariate which can be interpreted for example as expo-
observed. Parameter estimation is straightforward by sure in an observational study. Survival time is meas-
maximizing the log-likelihood function with respect to ured beginning with the day of diagnosis of lung cancer.
the unknown regression coefficients β and the parame- As true underlying models used for data generation we
ters of the assumed baseline distribution. Each software assumed (i) the Weibull additive hazard model and (ii)
that allows coding and maximizing such function, as the Log-Logistic additive hazard model as shown in
for example SAS via the NLMIXED procedure or R via equation (5) and (6). The true parameters of the distribu-
the optim-function, can be used to this task. tions equaled their estimates from the HALLUCA study
For practical application, it is possible to assume a wide
range of baseline distributions, including for example the – with aWB = 0.86, bWB = 1.77 for the Weibull distribu-
Exponential, Weibull, Gamma, Gompertz, Log-Normal tion, and
and Log-Logistic distribution. For instance, assuming the – with aLL = 1.06, bLL = 1.14 for the Log-Logistic
Weibull distribution as baseline distribution and includ- distribution.
ing a covariate xi , the cdf is then given by:
  aWB  In addition, we varied the number of participants per
ti
Fxi (ti ) = 1 − Sxi (ti ) = 1 − exp − − t i xi β . (5) study which were set to 50 or 200. Moreover, we distin-
bWB
guish between a smaller number of observed events per
Alternatively, using a Log-Logistic distribution as base- study of 60% (40% censoring) and a higher number of
line distribution, this leads to 80% observed events (20% censoring). Further, we eval-
−1 uated the true effect of the binary covariate for β = 0,
β = 0.8 and β = 1.6.
 
bLL
ti
+ 1
aLL (6)
Fxi (ti ) = 1 − Sxi (ti ) = 1 − .
exp(ti xi β) Data generation
Combining all parameters led to 24 different settings, i.e.,
aWB , bWB , aLL and bLL denote the distribution-specific 12 settings for each of (i) and (ii), for which data were
parameters. generated. For every setting, we simulated 1000 data
Results from applying the different baseline distribu- sets. Participants of each study were randomly allocated
tions can be compared via model selection criteria as the into two groups, as indicated by the binary covariate to
AIC or BIC. Parameters of main interest are finally the which the finally estimated regression coefficient β cor-
estimated regression coefficients and suitable transfor- responds. We assumed that our two groups are of equal
mations of the distribution parameters that have more size. Survival times for each participant were generated
intuitive interpretations, as the baseline mean or median from the true underlying distribution using inverse trans-
of the assumed distribution. form sampling. In case the covariate takes the value of
zero, the cdf Fxi (ti ) from which we had to sample simply
equals the Weibull or Log-Logistic distribution. Contra-
Simulation study
rily, if the covariate takes a value of one, the procedure is
To evaluate the parametric additive hazard model,
less trivial. In that case, it is not possible to invert the cdf
we conducted a simulation study. For comparison, we
analytically and we used numerical inversion to obtain
included the semi-parametric additive hazard model of
random numbers from the respective function [5]. The
Lin et al. [17]. The simulation study was implemented
likelihood of whether a study participant experienced an
using R (version 4.1.2), the full code is publicly available
event or if the event is censored was generated from a
on Zenodo [26].
Bernoulli distribution with respective success probability.
To guarantee uninformative censoring, we multiplied the
Setting original survival time with a uniformly distributed ran-
Our parameter settings were motivated by the Halle dom number in case the observation was censored.
Lung Carcinoma (HALLUCA) study which we also use
as exemplary application in the Example section (Exam- Estimation methods and outcomes
ple: the HALLUCA study) [3]. The data from the HAL- For parameter estimation we used the Weibull and Log-
LUCA study has been used in previous work proposing Logistic additive hazard model. The corresponding likeli-
an extension of relative survival models for clustered hood functions were manually coded and implemented
responses [16] and are also suitable for our purpose. in the statistical software R. Optimization was done using
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 4 of 11

the optim-function. True parameter values served as Bias

starting values for this procedure. As described above, we Appendix II of the Additional file 1 provides the
also applied the semi-parametric additive hazard model numerical results in terms of bias. Figure 1 illustrates
of Lin et al. [17] as comparative model. For this aim, we the corresponding estimates. If the estimated addi-
used the R-package addhazard [15] where the Lin and tive hazard model is consistent with the true underly-
Ying model is already implemented. Outcome of primary ing distribution, we observe the best performance, i.e.,
interest was the estimated regression coefficient β as this smallest median bias, with some exceptions. Further,
parameter is estimated by the parametric as well as by the also considering that the estimated additive hazard
semi-parametric model. The estimated parameters were model is consistent with the true underlying distribu-
compared in terms of bias, mean squared error (MSE) and tion, the additive hazard model slightly outperforms
empirical coverage. Moreover, we identified the number of the semi-parametric model by Lin and Ying. This holds
converged simulation runs to assess numerical robustness. true for all estimated models, leading to a bias between
−0.36 and 0.17. For most settings, the treatment effect
Results is slightly underestimated which leads to a negative
Simulation study bias. The observed underestimation also holds for the
In reporting the results of the simulation study, we con- Lin-Ying model which results in underestimation in
fine ourselves to the outcome of main interest, the esti- more than 65% of all cases. We observe a change in the
mated regression coefficient β. In the following, we give sign of the bias for settings where the Weibull-model
a brief overview of the results concerning bias, MSE, is estimated using the new additive hazard model and
empirical coverage and numerical robustness. The com- relying on data that was generated via a Log-Logistic
plete simulation results can be found in Appendix II to V distribution. This underlines that in these settings the
of the Additional file 1. Weibull distribution is not an appropriate choice for

Fig. 1 Box plots of estimated median bias for the regression coefficient β over 1000 simulations for each setting. Y-axis denotes each setting
with an ID consisting of the abbreviated true model (WBAH = Weibull additive hazard, LLAH = Log-Logistic additive hazard), the true β, the number
of participants per study and number of events (e.g., Weibull additive hazard model, true β = 0, number of patients = 50, number of events = 60%
results in “WBAH_0_50_0.6”). Left-most plot shows results for the Weibull (WB) additive hazard model, middle plot shows results for the Log-Logistic
(LL) additive hazard model and right plot shows results for the Lin-Ying (LY) model
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 5 of 11

modeling data that follows a Log-Logistic distribution, the true value of β , where generally the bias is closer
and an alternative choice of the underlying distribution to zero and its variability is reduced for smaller values
should be made. The most extreme positive bias of 0.17 of β . In most cases, the bias is smaller in settings with
is observed for the estimated Weibull model when the more events (80% compared to 60%). With regards to
Log-Logistic distribution, a true β of 0.8, 20% censor- the Weibull model, this only holds true if the estimated
ing and 200 patients were assumed for data generation. and true distribution are consistent.
The most extreme negative bias of −0.35 is observed
for the estimated Weibull and Log-Logistic model in MSE
the setting with data generated from the Weibull addi- Appendix III of the Additional file 1 shows the detailed
tive hazard model, a true β of 1.6, 40% censoring and results with regards to the MSE. Figure 2 visualizes
200 observations per sample. Generally, when data was the findings. Generally, there are hardly any differ-
generated assuming a Weibull distribution, all three ences for the MSE across the three estimated models
models perform similarly. Though, if a Log-Logistic dis- for the respective settings. In some settings we observe
tribution is assumed as true, the estimated Log-Logistic an MSE of 0.00. These results indicate that the model
and Lin-Ying model perform similarly, while results works very well at minimum in these contexts. Contra-
from the Weibull model slightly differ. The number of rily, we observe a maximum error of 0.45 for the Lin-Ying
participants modeled per study minorly influences the model in one setting (true model based on Log-Logistic
median bias. The variability of the bias is smaller for distribution, true β = 1.6, number of observations = 50,
higher numbers of observations (200 versus 50). Fur- number of events 60%). As for the bias, in most cases
thermore, the estimated models are sensitive towards the MSE is slightly smaller when the estimated model

Fig. 2 Box plots of estimated median mean squared error (MSE) for the regression coefficient β over 1000 simulations for each setting. Y-axis
denotes each setting with an ID consisting of the abbreviated true model (WBAH = Weibull additive hazard, LLAH = Log-Logistic additive hazard),
the true β, the number of participants per study and number of events (e.g., Weibull additive hazard model, true β = 0, number of patients = 50,
number of events = 60% results in “WBAH_0_50_0.6”). Left-most plot shows results for the Weibull (WB) additive hazard model, middle plot shows
results for the Log-Logistic (LL) additive hazard model and right plot shows results for the Lin-Ying (LY) model
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 6 of 11

is consistent with the true underlying distribution. For and Log-Logistic model mostly outperform the Lin-Ying
almost any setting, it holds true that the median MSE model, respectively. Overall, the empirical coverage
of the Weibull and Log-Logistic model is better than the often falls below 95%. Precisely, this is the case for 75%
MSE of the respective Lin-Ying model. The MSE seems to of the settings for the Weibull, for 71% of settings for the
be negatively correlated with the number of participants Log-Logistic and for 71% of all settings for the Lin-Ying
per study, as for settings with 200 observations per sam- model. The highest coverage of 100% is observed for the
ple (versus 50 observations) the MSE is smaller. Similar Log-Logistic model when data was generated from a Log-
results are found with regards to the assumed number of Logistic distribution with β equal to 0, 60% events and
events: the MSE decreases when more events (80% com- 200 observations. On the other hand, the smallest cov-
pared to 60%) are observed. Furthermore, for larger val- erage of 57.7% is observed for the Log-Logistic for data
ues of β, the MSE and its variability increase for each of generated from a Weibull distribution with β equal to
the estimated models. 0, 80% events and 200 observations. Similar to bias and
MSE, results depend on the true value of β. Coverage
Empirical coverage results for the estimated models seem to improve when
Appendix IV of the Additional file 1 provides detailed β is closer to zero. With regards to the number of events,
results about the models’ empirical coverage. Figure 3 results correlate positively with the event probability.
illustrates the corresponding estimates. In terms of Precisely, coverage results improve in settings where the
empirical coverage (on the 95% level), the performances event probability is higher (80% versus 60%). The number
of the three estimated models are relatively similar with of participants per sample seems to have a slight influ-
some exceptions for a handful of settings. For the major- ence on coverage results. However, the correlation of
ity of the settings, the models again perform best when coverage and sample size is less clear for all three models
the estimated model and the true underlying distribution compared to the correlation of the bias and study size.
are equal. If the distributions are consistent, the Weibull
Numerical robustness
Appendix V of the Additional file 1 shows the results
concerning the models’ convergence. The numerical
robustness is very satisfactory for all estimated time-to-
event models. The Lin-Ying model performs most sta-
ble in terms of numerical robustness, always returning
1000 converged simulation runs. With respect to the
Weibull model, the results are also very satisfactory.
Only for two settings where data was generated from
a Log-Logistic distribution, the model converges for
999 instead of 1000 simulation runs. The Log-Logistic
model performs slightly less numerically robust with
a minimum of 983 converged runs. In 75% of the set-
tings, the Log-Logistic model achieves convergence for
all 1000 generated data sets.

Example: the HALLUCA study

For illustration purposes, we used data from the HAL-
LUCA study, an epidemiological study that investigated
medical care for lung cancer patients in the region of
Fig. 3 Bar plot showing relative frequency of coverage Halle (Saale) in the eastern part of Germany [3, 26]. In
for the regression coefficient β over 1000 simulations for each setting. the HALLUCA study, a total of 1696 lung cancer patients
Y-axis denotes each setting with an ID consisting of the abbreviated were observed between April 1996 and September 2000.
true model (WBAH = Weibull additive hazard, LLAH = Log-Logistic 1349 patients (79.5%) died until the end of follow-up.
additive hazard), the true β, the number of participants per study
and number of events (e.g., Weibull additive hazard model, true
The median survival time in the population was 284 days
β = 0, number of patients = 50, number of events = 60% results (0.78 years). For the analysis reported here, we focused
in “WBAH_0_50_0.6”). Left-most plot shows results for the Weibull on a dichotomized version of the TNM-scale for the clas-
(WB) additive hazard model, middle plot shows results sification of malignant tumors as a predictor, where 739
for the Log-Logistic (LL) additive hazard model and right plot shows patients had a TNM-scale of IIIb or smaller (TNM < IV),
results for the Lin-Ying (LY) model
and 621 patients had a TNM-scale of IV (TNM IV). The
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 7 of 11

336 remaining patients for whom the TNM-scale was not (95% confidence interval [1.56; 1.81]). The observed
reported, were deleted. The starting point for the survival median survival time is 1.10 years [0.96; 1.24]. The fig-
definition relates to the day of diagnosis of lung cancer, ure shows that the estimated survival curves from the
while the end point of a patient’s survival time relates to additive hazard model fit well with the Kaplan-Meier
death or the study end. curves. For an increasing time t there are greater differ-
The Kaplan-Meier estimates for the two TNM classes ences between the curves, especially for the Exponential,
and the estimated survival probabilities using the para- Weibull, Gamma and Gompertz model. However, for the
metric additive hazard model with various distributions Log-Normal and Log-Logistic distribution the fit is also
are given in Fig. 4. For the reference group (TNM < IV), quite well for larger survival times.
for which we can estimate the parametric baseline dis- In Table 1, the results of the additive hazard model
tribution, the observed mean survival time is 1.69 years for six different baseline distributions are given. All

Fig. 4 Estimated survival probabilities from the parametric additive hazard model using different distributions compared to the Kaplan-Meier
curves
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 8 of 11

Table 1 Estimates for the HALLUCA data set with 95% confidence intervals
Distribution β TNMIV Baseline Mean Baseline Median No. of distribution BIC
parameters

Exponential 0.84 [0.72; 0.97] 1.77 [1.62; 1.92] 1.22 [1.12; 1.33] 1 2441.9
Weibull 0.78 [0.65; 0.91] 1.91 [1.70; 2.11] 1.15 [1.03; 1.26] 2 2427.1
Gamma 0.79 [0.67; 0.92] 1.85 [1.66; 2.03] – 2 2431.8
Gompertz 0.84 [0.72; 0.97] – 1.22 [1.12; 1.33] 2 2449.3
Log-Normal 0.80 [0.68; 0.92] 3.44 [2.75; 4.13] 1.07 [0.95; 1.19] 2 2450.3
Log-Logistic 0.78 [0.65; 0.91] 7.90 [3.39; 12.40] 1.06 [0.94; 1.18] 2 2417.6

models include a single binary covariate TNMIV, Figure 5 compares the relative survival probability from
with (TNM < IV) being the reference group. Thus, we the parametric additive hazard model assuming different
model the additive hazard of being in TNM stage IV. baseline distributions. The relative survival curves, taken
For evaluation purposes, we report the baseline mean as the ratio of survival in a group of individuals with
and baseline median which can be compared to the TNM stage IV in comparison to the survival of a corre-
empirical mean and median reported in the original sponding population with TNM < IV (see (4)), are rather
HALLUCA study. To compare the model fit, we also similar between the various baseline distributions. It can
report the Bayesian information criterion (BIC), where easily be seen that survival strictly decreases over time,
a smaller value indicates a better fit. The estimated with a relative survival probability almost equal to zero
regression coefficients for TNMIV show that being in after five years.
stage IV increases the hazard by approximately 0.8.
Being in an additive setting, the estimated β-value is Discussion
added to the baseline distribution. That means being in In this article, we propose a new parametric additive haz-
the worse tumour group, the risk of dying rises because ard model for time-to-event outcomes. Generally, para-
our model describes by how much the risk of dying metric survival models are advisable as they are simpler,
changes in additive manner, not multiplicative. Alter- more informative and robust than non- or semi-paramet-
natively, the coefficient could be interpreted in abso- ric models [8, 10, 20]. In addition, they have hazard, sur-
lute terms, saying that when observing 100 people for vival and probability density functions directly available.
one year, there are an additional 80 deaths among can- We showed that the model is valid from a theoretical
cer patients in TNM class 4 compared to those with a point of view. Further, in a simulation study as well as for
lower TNM class. The results change slightly depend- an example from lung cancer research, we demonstrated
ing on the different distributions. Overall however, esti- that the model showed that the model works in practice.
mates from the parametric models are similar. In terms Convincingly, our proposed model can be implemented
of the BIC, the Log-Logistic, Weibull and Gamma using any standard software that allows coding and maxi-
model fit best, returning almost identical parameter mizing a likelihood function. The corresponding SAS
estimates and confidence intervals for βTNMIV . With and R codes of our applications are publicly available on
respect to the baseline location measures, the differ- Zenodo [26].
ences between the distributions are rather large, rang- The facilitated interpretation of our model is one of its
ing from 1.77 to 7.90 for the baseline mean, and from most convincing advantages and is important for exam-
1.06 to 1.22 for the baseline median. For example, for ple from a clinical point of view. For instance, allowing
the Log-Logistic approach, the model with the best to communicate outcomes of the additive hazard model
BIC, we get a large baseline mean of 7.90 with a broad in relative and absolute terms increases the comprehen-
95%-confidence interval of [3.39; 12.40]. With refer- sibility of its results. Moreover, including time-independ-
ence to Fig. 4, this can be explained by the flat survival ent (constant) covariates, our approach is advantageous
curves at the end of the prediction range. Therefore, over the additive hazard model of Aalen, since there,
the baseline mean is necessarily larger with a broader model coefficients often are difficult to understand as
confidence interval. However, regarding the baseline they change repeatedly over time and cannot be sum-
median, the Log-Logistic model yields an estimate of marised easily. Furthermore, by modeling regression
1.06 which is in line with the empirical median of 1.10 coefficients and distribution parameters that have more
from the Kaplan-Meier approach. The same holds true intuitive interpretations, our approach overcomes inter-
for the Log-Normal and Weibull model which achieve pretational, mathematical and technical problems that
estimates that are also close to the empirical median. arise when estimating hazard ratios as generic effect
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 9 of 11

Fig. 5 Relative survival probabilities from the parametric additive hazard model using different distributions. The gray shaded area depicts
estimated confidence intervals

estimates. For practical application, several different small sample sizes and also rather moderate number of
distributions such as the Exponential, Weibull, Gamma, events. Consequently, it is not unlikely that the results
Gompertz, Log-Normal and Log-Logistic distribution of all models investigated in the simulation are biased
can be implemented and compared via model selection in a certain extent and no optimal solution was found.
criteria as the AIC or BIC. Moreover, parameter esti- For practical application, it should be considered posi-
mation of our model by the maximum likelihood prin- tively that the treatment effect is more likely to be
ciple is straightforward and can be accomplished with underestimated (versus overestimated). Thus, slightly
standard statistical software that allows for maximizing biased inference in practical research may be less prob-
a likelihood function. Additionally, the proposed model lematic than in a case of overestimation. Relying on
is highly flexible in terms of possible extensions such as the asymptotical properties and a correct implemented
modeling correlated data by including random effects in maximum likelihood estimation, we are confident that
the linear predictor, modeling non-linear covariate effects with an increased number of events and larger sample
by splines, or by specifying baseline distributions that sizes, i.e., larger number of participants per study, the
have more than two parameters. Maximum flexibility bias and its variability would asymptotically approach
with respect to the baseline distribution can be obtained 0. Further, in practical research, it would evidently be
by using a piecewise-constant model, i.e., by dividing the an option to tune the parameters, settings and start-
observation period a-priori into intervals and assuming ing values of software functions used in the estimation
an exponential distribution in each of these intervals. to increase robustness and ensure convergence. For
Evidently, there are some limitations with regards to our study, we refrained from a comparison to the Cox
our simulation study. The simulation was motivated model and instead used the Lin-Ying model for model
by a real-world data example and thus offers a realistic evaluation. This is due to the fact that (i) the Lin-Ying
setting, which we varied in a certain range of potential model also relies on the additive hazards assump-
scenarios. However, these scenarios consider relatively tion and is thus a more direct competitor and (ii) our
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 10 of 11

proposed model and the Cox model are inherently dif- parametric procedures require starting values for the
ferent and in a certain way incomparable. For instance, optimization which may sometimes be complicated
converting the additive into the multiplicative effect or problematic to define in practice. In case a distri-
estimate is impossible. Further, the Cox model assumes bution can be confidently specified to the data, para-
multiplicative effects. Therefore, if data are generated metric models will usually be more informative than
using an additive model, it is obvious that the addi- semi- or non-parametric approaches. However, if this
tive model will perform superior to the Cox model in is not the case and the assumed distribution is false,
terms of the effect estimator and its corresponding results and conclusions are likely to be biased. For
confidence interval. Vice versa, when generating data practical application, we recommend that the choice of
using the Cox model, we presume that our proposed the baseline distribution should be based on evidence
parametric model would perform worse as it relies on from literature and previous research. However, if this
the additivity assumption. However, it is left open for is not applicable, and selection criteria such as the BIC
future research to assess whether p-values are con- or AIC is used for a data-based decision, this must be
sistent for the different models, even if they are simu- acknowledged in the reporting of the study and should
lated from the “wrong” model. Overall, the choice of be critically discussed. In that context, users should
using the Cox model or generally any additive model be cautious with any interpretation and should be are
is dependent on the primary effect measure and study aware that confidence intervals for instance are poten-
question. If the hazard ratio is of interest, then the Cox tially estimated too narrow. Besides, and by defini-
model would be preferable. If that is not the case, an tion, our model assumes that the covariates act in an
additive model might be the better choice as it offers additive way to the baseline hazard function. Appar-
more easily interpretable effect measures. Further, the ently, this assumption may be doubtful in some appli-
choice must be made under consideration of whether cations and in that case, other models may be more
the user wants to assume multiplicativity or additivity. appropriate.
There is a number of methods for checking the assump- With regards to the framework of phases of methodo-
tions and model fit for the Cox model, however, not all logical research Heinze et al. [13], the proposed additive
can be extended to the additive model [6]. For instance, hazard model currently belongs to early stages of meth-
we are not aware of any methods to confirm the valid- odological development, i.e., phase I or II. The aim here
ity of the additivity assumption or of the multiplicativ- was to introduce a new idea, as well as to demonstrate its
ity assumption, which would be the pendant specified validity and its potential to improve on existing methods.
in the Cox model. Although the direct comparison with Into the bargain, we derived the new methodological idea
a Cox model is generally questionable and outside the while providing, logical reasoning and proofs of empirical
scope of this article, it nevertheless may remain inter- evidence through a real-world data example and a simu-
esting in case of insufficient knowledge about the effect lation study in a (yet) relatively narrow but suitable target
of the covariates and whether to choose an additive setting. Carefully planned method comparison studies
or multiplicative model. We are unaware of an “over- that investigate the model in future works could advance
all”, combining model that contains both, additive and the proposed additive hazard model to later stages of
multiplicative models as a special case, and that has methodological development and would be of great value
an additional parameter which measures the type of for future users and the scientific community as a whole.
relationship (i.e., additive, multiplicative or something Therefore, there is a need for future studies that explore
in between). The development of such a framework the empirical properties of our model in a wider range
or model may be an interesting and valuable topic for of problems, highlight its advantages and limitations,
future research. and possibly uncover previously unknown behavior (e.g.,
With regards to the weaknesses of the proposed in simulations with wide range of scenarios and differ-
additive hazard model, our method suffers from the ent outcome types and realistic or complex comparative
drawbacks that limits parametric versus semi- or example data analyses).
non-parametric models. Non-parametric approaches
generally require fewer assumptions about the data.
Therefore, these may prove better when the true dis- Conclusion
tribution is unknown and/or cannot be easily approxi- To summarize, the proposed parametric additive hazard
mated. Vice versa, parametric methods are inherently model can serve as a powerful tool to analyze time-to-
dependent on the distribution chosen for estima- event outcomes. By definition, it simplifies interpretation,
tion and on the assumption that this distribution is facilitates parameter estimation and permits greater flex-
correctly specified. Further, it must be noted that ibility than most existing and commonly used methods.
Voeltz et al. BMC Medical Research Methodology (2024) 24:48 Page 11 of 11

Supplementary Information 9. Cox DR. Regression models and life-tables. J R Stat Soc Ser B (Methodol).
1972;34(2):187–202.
The online version contains supplementary material available at [Link]
10. Crowther MJ, Lambert PC. A general framework for parametric survival
org/​10.​1186/​s12874-​024-​02180-y.
analysis. Stat Med. 2014;33(30):5280–97.
11. Dickman PW, Sloggett A, Hills M, Hakulinen T. Regression models for rela-
Supplementary material 1. tive survival. Stat Med. 2004;23(1):51–64.
12. Greenland S. Absence of confounding does not correspond to collaps-
ibility of the rate ratio or rate difference. Epidemiology. 1996;7(5):498–501.
Acknowledgements 13. Heinze G, Boulesteix AL, Kammer M, Morris TP, White IR, Simulation
Not applicable. Panel of the STRATOS initiative. Phases of methodological research
in biostatistics-Building the evidence base for new methods. Biom J.
Authors’ contributions 2022;66(1):2200222.
DV is the corresponding author of this work and, as such, had full access to all 14. Hernán MA. The hazards of hazard ratios. Epidemiology (Camb Mass).
data in the study and takes responsibility for the integrity of the data and the 2010;21(1):13.
accuracy of the analysis. AH and OK derived the methodology. AH supervised 15. Hu J, Breslow N, Chan, G. addhazard: Fit Additive Hazards Models for Sur-
the study and, together with DV, developed the theoretical concept and vival Analysis. 2017. R package version. R package version 1.1.0. [Link]
outline of the present study. With support of AS and AF, AH and DV performed CRAN.R-​proje​ct.​org/​packa​ge=​addha​zard. Accessed 19 Feb 2024.
the formal analysis. DV and AH wrote the manuscript. All authors discussed 16. Kuß O, Blankenburg T, Haerting J. A relative survival model for clustered
the results, reviewed and edited the manuscript and the statistical analyses. responses. Biom J J Math Methods Biosci. 2008;50(3):408–18.
17. Lin DY, Ying Z. Semiparametric analysis of the additive risk model. Biom-
Funding etrika. 1994;81(1):61–71.
Open Access funding enabled and organized by Projekt DEAL. The authors 18. Lindsey J. A study of interval censoring in parametric regression models.
received no financial support for the research, authorship, and/or publication Lifetime Data Anal. 1998;4(4):329–54.
of this article. 19. Martinussen T, Vansteelandt S. On collapsibility and confounding bias in
Cox and Aalen regression models. Lifetime Data Anal. 2013;19(3):279–96.
Availability of data and materials 20. Nardi A, Schemper M. Comparing Cox and parametric models in clinical
The datasets generated and/or analysed during the current study are available studies. Stat Med. 2003;22(23):3597–610.
in the Zenodo repository, [Link] [26]. 21. Reid N. A conversation with Sir David Cox. Stat Sci. 1994;9(3):439–55.
22. Schaubel DE, Wei G. Fitting semiparametric additive hazards models
using standard statistical software. Biom J. 2007;49(5):719–30.
Declarations
23. Sjölander A, Dahlqwist E, Zetterqvist J. A note on the noncollapsibility of
rate differences and rate ratios. Epidemiology. 2016;27(3):356–9.
Ethics approval and consent to participate
24. Stensrud MJ. Interpreting Hazard Ratios: Insights from Frailty Models.
Not applicable.
2017. arXiv preprint arXiv:170106014. [Link]
Accessed 19 Feb 2024.
Consent for publication
25. Sutradhar R, Austin PC. Relative rates not relative risks: addressing a wide-
Not applicable.
spread misinterpretation of hazard ratios. Ann Epidemiol. 2018;28(1):54–7.
26. Voeltz D, Hoyer A. Parametric additive hazard model for time-to-event
Competing interests
analysis. Zenodo. 2022. [Link]
The authors declare no competing interests.
27. Zahl PH, Aalen OO. Adjusting and comparing survival curves by means of
an additive risk model. Lifetime Data Anal. 1998;4(2):149–68.
Received: 9 May 2023 Accepted: 14 February 2024
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in pub-
lished maps and institutional affiliations.

References
1. Aalen OO. A linear regression model for the analysis of life times. Stat
Med. 1989;8(8):907–25.
2. Aalen OO, Cook RJ, Røysland K. Does Cox analysis of a randomized
survival study yield a causal treatment effect? Lifetime Data Anal.
2015;21(4):579–93.
3. Bollmann A, Blankenburg T, Haerting J, Kuss O, Schütte W, Dunst J,
et al. Survival of patients in clinical stages I-IIIb of non-small-cell lung
cancer treated with radiation therapy alone. Strahlenther Onkol.
2004;180(8):488–96.
4. Bradburn MJ, Clark TG, Love SB, Altman DG. Survival analysis part II:
multivariate data analysis-an introduction to concepts and methods. Br J
Cancer. 2003;89(3):431–6.
5. Brinks R, Landwehr S, Fischer-Betz R, Schneider M, Giani G. Lexis diagram
and illness-death model: simulating populations in chronic disease
epidemiology. PLoS ONE. 2014;9(9):e106043.
6. Cao HL. A comparison between the additive and multiplicative risk mod-
els. Université Laval; 2005. [Link]
sne/​CaoFi​nal.​pdf. Accessed 12 Feb 2024.
7. Collett D. Modelling survival data in medical research. CRC press; 2015.
[Link]
alysis/​Colle​tt_1_​3.​PDF. Accessed 19 Feb 2024.
8. Columbia. 2022. [Link]
ation-​health-​metho​ds/​time-​event-​data-​analy​sis#. Accessed 8 Jan 2024.