Modelbased estimation of measures of association for timetoevent outcomes
 Federico Ambrogi^{1}Email author,
 Elia Biganzoli^{1, 2} and
 Patrizia Boracchi^{1}
DOI: 10.1186/147122881497
© Ambrogi et al.; licensee BioMed Central Ltd. 2014
Received: 23 December 2013
Accepted: 25 July 2014
Published: 9 August 2014
Abstract
Background
Hazard ratios are ubiquitously used in time to event applications to quantify adjusted covariate effects. Although hazard ratios are invaluable for hypothesis testing, other adjusted measures of association, both relative and absolute, should be provided to fully appreciate studies results. The corrected group prognosis method is generally used to estimate the absolute risk reduction and the number needed to be treated for categorical covariates.
Methods
The goal of this paper is to present transformation models for timetoevent outcomes to obtain, directly from estimated coefficients, the measures of association widely used in biostatistics together with their confidence interval. Pseudovalues are used for a practical estimation of transformation models.
Results
Using the regression model estimated through pseudovalues with suitable link functions, relative risks, risk differences and the number needed to treat, are obtained together with their confidence intervals. One example based on literature data and one original application to the study of prognostic factors in primary retroperitoneal soft tissue sarcomas are presented. A simulation study is used to show some properties of the different estimation methods.
Conclusions
Clinically useful measures of treatment or exposure effect are widely available in epidemiology. When time to event outcomes are present, the analysis is performed generally resorting to predicted values from Cox regression model. It is now possible to resort to more general regression models, adopting suitable link functions and pseudo values for estimation, to obtain alternative measures of effect directly from regression coefficients together with their confidence interval. This may be especially useful when, in presence of time dependent covariate effects, it is not straightforward to specify the correct, if any, time dependent functional form. The method can easily be implemented with standard software.
Keywords
Survival analysis Transformation models Pseudovalues Link functions Numbers needed to treatBackground
Measures of disease frequency and measures of associations derived from them are among the basic building blocks of biostatistics and epidemiology. The appropriateness of the use of a specific measure of association may depend on the study objectives and design. Sometimes, however, the use of specific measures of association depends also on the statistical methods available for estimation. For example, in epidemiology, a debated subject concerns the use of odds ratios, estimated through logistic regression, in cohort studies of common outcomes [1, 2].
When timetoevent outcomes are analyzed, the presence of censoring calls for specific methods of analysis [3]. The evaluation of the effect of a treatment in a controlled trial can be performed through the graphical display and comparison of KaplanMeier curves at selected times, when adjustment is not required. Otherwise, the measure of effect generally considered is the adjusted hazard ratio estimated by means of Cox proportional hazard model, ([4, 5]).
However, the clinical literature in randomized controlled trials suggests the use of absolute measures of effect to assess the effects of a treatment, such as risk difference or the number needed to be treated, which are better suited than relative measures of effect for clinical decision support, see [6–10] among others. Schechtman highlights how relative measures are appropriate for summarizing the evidence while absolute measures for the concrete application in a clinical setting, [11].
The need for alternatives to hazard ratios, a relative measure of effect based on (instantaneous) incidence rates, is increasing in medical/epidemiological literature. In particular, the possibility to provide absolute measures of association computed using adjusted survival curves was explored in literature [12–14].
To precisely define the different measures of association in timetoevent applications, it is useful to distinguish between the risk of the event, F(t), i.e. the probability of a patient having the event over a defined followup time, and the event rate, λ(t), i.e. the number of events in a specified followup interval divided by the time at risk accumulated during the interval. The instantaneous hazard rate is obtained when the interval length approaches 0. The hazard rate at time t refers to the population survived until time t, while the risk refers to the whole population. The measures considered in this paper refer to ratios and differences between the risk of event of different groups of subjects. Let F_{1}(t) and F_{2}(t) be the event risk by t in two groups of subjects, (exposed and non exposed, standard treatment and new treatment), then we define the risk difference as R D(t)=F_{1}(t)−F_{2}(t). It is useful to translate R D(t) expressed as a percentage measure in a measure more sensible form a clinical perspective. To this end it is usual to use the number needed to be treated N N T(t)=1/R D(t) which is interpreted as the expected number of patients needed to be treated to avoid one additional death compared to the untreated. The measure has its roots in clinical trial literature and was extended in an epidemiological framework as the number needed to be exposed, N N E(t), i.e. the expected number of subjects to be exposed to have one additional event compared to the unexposed. In observational studies, an alternative definition of N N E(t) is the exposure effect among the unexposed, while the exposure impact number, E I N(t), describes the effect of removing the exposure among the exposed [15–17]. It is also interesting to define a relative risk, $\mathit{\text{RR}}\left(t\right)=\frac{{F}_{1}\left(t\right)}{{F}_{2}\left(t\right)}$, to be contrasted with the hazard ratio, $\mathit{\text{HR}}\left(t\right)=\frac{{\lambda}_{1}\left(t\right)}{{\lambda}_{2}\left(t\right)}$. The different measures of effect are in general timevarying. In certain situations, however, they are estimated as constant through followup, as it happens for example with the Cox proportional hazard model for the hazard ratio. When the measure is assumed to be constant during followup the time dependence is omitted (i.e., H R(t) is written as HR).
The purpose of this paper is to provide an outline of the methods generally adopted to estimate adjusted summary measures of associations, different from the hazard ratio, in timetoevent studies and to present a new method based on transformation models. The focus of this paper is not to provide guidance about which association measure should be used in different situations, but simply to provide an estimation method. Moreover, particular attention will be given to the estimate of adjusted R D(t). In fact, absolute measures of association are particularly advocated in survival analysis, to be combined with the generally used hazard ratio. A small simulation study is provided to show a preliminary evaluation of the properties of the different estimation procedures. Two examples are then developed. The first concerns literature data on a clinical trial on 506 prostate cancer patients [18]. The goal is to estimate the treatment R D(t) and N N T(t) with their confidence intervals. A comparison of the model based estimated R D(t) with that obtained with the classical corrected group prognosis method [14] is provided. The second application concerns an observational study of prognostic factors in primary retroperitoneal soft tissue sarcomas [19].
Methods
Computation of association measures
In many situations a researcher is interested in providing adjusted estimates of covariate associations with the outcome. In observational studies (involving no randomization) the exposure effect has to be adjusted for known confounders. Also in randomized controlled trials (RCT) the use of adjusted estimates is suggested for example to account for potential covariate imbalances or since prognostically relevant covariates were considered for a stratified randomization [20–22]. In these cases, Cox regression is widely used to adjust the estimated association between the covariate of interest and outcome for the other covariates.
where $\u015c\left(t\rightT=1,A=40,G=1)$ is the estimated survival probability for a male 40 years old in the experimental treatment group, given by ${\left[\phantom{\rule{0.3em}{0ex}}{\u015c}_{0}\right(t\left)\right]}^{\mathit{\text{exp}}(\widehat{\alpha}+\widehat{\beta}\phantom{\rule{2.77626pt}{0ex}}40+\widehat{\gamma})}$, while $\u015c\left(t\rightT=0,A=40,G=1)$ is the estimated survival probability for a male 40 years old, in the control group, given by ${\left[\phantom{\rule{0.3em}{0ex}}{\u015c}_{0}\right(t\left)\right]}^{\mathit{\text{exp}}(\widehat{\beta}\phantom{\rule{2.77626pt}{0ex}}40+\widehat{\gamma})}$. ${\u015c}_{0}\left(t\right)$ is the baseline survivor function from a Cox proportional hazards model estimated according to one of the available methods [23, 24].
In order to obtain adjusted measures of association, different from the hazard ratio, the Cox proportional hazard model is used to estimate adjusted survival curves [25] as outlined in the following paragraphs.
Average covariate method
The simplest approach for obtaining adjusted survival curves is the average covariate method. The mean values among the study patients of the covariates used for adjustment are plugged into the multivariable Cox proportional hazard model. Considering the above example, if the mean age of subjects in study is 45 and 30% males are included, the adjusted $\widehat{\mathit{\text{NNT}}}\left(t\right)$ for the treatment will be $\frac{1}{\u015c\left(t\rightT=1,A=45,G=0.3)\u015c(tT=0,A=45,G=0.3)}$. The average covariate method was once popular and largely adopted, due to its simplicity, but it was severely criticized [25, 26].
In fact it involves the averaging of categorical covariates, such as gender, which is difficult to understand. Moreover the method provides an estimate of the measure of effect for an hypothetical average individual and not a population averaged estimate.
Corrected group prognosis method and developments
An alternative idea is the corrected group prognosis method (CGPM), [14, 27, 28]. In the following, the CGPM to estimate R D(t), as described by Austin [14], is outlined:

A Multivariable Cox (or fully parametric) regression is used for the treatment and the covariates.

For each subject, the predicted survival probabilities, at the times of interest, are estimated using the multivariable model, assuming each subject is in the experimental treatment group; then the predictions are averaged;

the same predictions are obtained and averaged assuming each subject is in the control group.

the difference between the averaged predicted probabilities between experimental and control group is an estimate of the adjusted R D(t) for the experimental treatment at the specified times.
Pointwise confidence intervals of the obtained R D(t) estimates may be computed via bootstrap resampling [14]. For each bootstrap sample, i.e. a sample of the same size of the original one and randomly drawn with replacement from it, the R D(t) is computed according to the procedure outlined. A non parametric bootstrap 95% pointwise confidence interval is obtained resorting to the 2.5^{ t h } and 97.5^{ t h } percentiles of the obtained R D(t) bootstrap distribution.
The CGPM can be applied in principle to whatever regression model and an adequate model must be chosen. Considering, for example, the Cox regression model, in presence of time dependent covariate effects, an interaction of the covariates with a prespecified function of time should be specified, in order to estimate H R(t) varying during followup time. It is important to remark that it is not always easy to specify an adequate model in presence of time dependent covariate effects. In fact it is not always obvious how to model the time dependence itself. In general simple functions of time (linear or logarithm) or more flexible alternatives are used, [29].
To allow the estimation the data set must be augmented as it is done for true timedependent covariates [30]. It is to be remarked that, although the use of predicted values from regression models is simple from a practical point of view, the standard way to obtain summary measures of effect and their confidence interval is to use directly regression model coefficient estimates. The CGPM applied to the Cox model will be used for comparison with the method here proposed and described in the following section.
Laubender and Bender, [13], proposed different averaging techniques to estimate relevant impact numbers in observational studies using Cox model. For the purpose of illustration, let us consider the same example as before simply considering an exposure (E) instead of treatment. To obtain an estimate of N N E(t) it is possible to average predictions considering the subjects as if they were unexposed and as if they were exposed and taking the difference. As the distributions of the covariates used for adjusting are in general different in the exposed and unexposed groups, two different measures should be considered. Specifically, the estimate of the N N E(t) is obtained considering the unexposed subjects only, while E I N(t) is obtained considering the exposed subjects only. A comparison of the model based estimated R D(t) with that obtained through different averaging techniques, namely N N E(t) and E I N(t) [13, 15], is provided in the second example. However, the focus of the paper is not the comparison of different averaging techniques which are provided only for illustrative purposes. In particular, only the estimates obtained through the averaging performed over the whole population are compared with those based on transformation models methods.
Modelbased estimates of association measures
Adjusted modelbased estimates of measures of association can be obtained resorting to a general class of regression models used in Survival Analysis called transformation models [31].
Pseudo values
Considering the previous example, the transformation model can be written as g(S(tT,A,G))=g(S_{0}(t))+α T+β A+γ G.
where n is the sample size, $\u015c\left(t\right)$ is the survival probability based on the KaplanMeier estimator using the whole sample and ${\u015c}_{i}\left(t\right)$ is the survival probability obtained by deleting the i subject from the sample. When no censoring is present in the data, the pseudo values for subject i at time t is simply 1 if the subject is alive at t, while it is 0 if the event happened by t. Suppose to have an exposed male, 40 years old, which dies after 30 months of followup. The pseudo values computed at 12, 24 and 36 months are equal to 0, 0 and 1 respectively. The times at which the pseudovalues are computed are called pseudotimes.
When censoring is present in the data, pseudovalues are still defined for each subject (even those censored) and for each time, but the values may also be less than 0 or greater than 1 (See [33]; page 5310–11 for further details on the properties of pseudovalues).
In general, to allow inference on the entire survival curve, M (greater than 5) pseudo times are used, considering, for example, the quantiles of the unique failure time distribution. As M pseudo values are computed for each subject, an augmented data set is created with M observations for each subject.
Transformation models and association measures
The pseudovalues are then used as responses in a regression model for longitudinal data, where time is a covariate. As no explicit likelihood is available for pseudovalues, generalized estimating equations (GEE), [34], are used accounting for the correlation of the pseudovalues within each subject. The cluster robust variancecovariance is used for hypothesis testing using Wald tests. In general an independence working variancecovariance matrix can conveniently be used in the estimation process [32].
In order to model g(S_{0}(t)), the transformed baseline survival function, the standard procedure is to insert in the regression model indicator functions for each pseudotime. If all event times would be used to compute the pseudovalues, the insertion of indicator functions would result in a non parametric representation of the (transformed) baseline survival, as in the Cox model. In general only a small number of pseudotimes are used obtaining a parametric baseline representation. As an alternative, spline functions can be inserted in the regression model, as did [35] in a nonpseudovalues framework.
where B_{1}(t) and B_{2}(t) represent the first and second spline bases for time t. For example, if a restricted cubic spline basis is used with three knots at k_{1},k_{2},k_{3}, then B_{1}(t)=t and ${B}_{2}\left(t\right)={(t{k}_{1})}_{+}^{3}\frac{{(t{k}_{2})}_{+}^{3}({t}_{3}{t}_{1})}{({t}_{3}{t}_{2})}+\frac{{(t{k}_{3})}_{+}^{3}({t}_{2}{t}_{1})}{({t}_{3}{t}_{2})}$, where, for example, ${(t{k}_{1})}_{+}^{3}$ is equal to (t−k_{1})^{3} if t>k_{1}, otherwise is 0. Knots are chosen at quantiles of the failure time distribution. In the case of 3 knots the quantiles commonly suggested are 0.1, 0.5 and 0.9, [36]. To choose the complexity of the spline the QIC, [37], an information criterion proposed for generalized estimating equations, can be used. A less formal strategy is the graphical comparison between the KaplanMeier marginal survival probability and the marginal probability obtained from the transformation model without covariates. Such a procedure will be used in the examples.
The first part of the model, ϕ_{0}+ϕ_{1}B_{1}(t)+ϕ_{2}B_{2}(t), provides a parametric representation of the (transformed) baseline survival function, g(S_{0}(t)), during followup time.
The coefficients α, β and γ represent the covariate effects expressed as differences in the Survival probability, transformed by g associated with a unit increase in the covariates. Let us consider such an issue in detail. When g is the logit link function, a proportional odds model is estimated. Accordingly, α, β and γ represent the logarithm of the ratio of the odds of surviving associated with the change of one unit in the covariates. Such an effect is constant through followup times. The exponentiation of the parameter estimates represent therefore the ratio of the odds of surviving. Similarly, the logarithmic link produces a proportional risks model and the e x p(α), e x p(β) and e x p(γ) represent the ratio of the survival probabilities (Relative Risks, RR). The identity link produces a constant survival difference model: α, β and γ represent the adjusted differences in survival probabilities (risk differences, RD). A constant difference model through followup is often not practical as a model such that at the beginning of the followup the survival curves start at 1 and then, eventually, become different. However, it is to be noted that the first pseudotime is never placed at time 0, but later on the followup time scale. In Figure 1 an example of the model based RD estimate with pointwise confidence intervals, constant through time, is reported in the right panel. The constant model estimated RD can be used to obtain a constant estimate of NNT by inversion. In the case of treatment T: $\widehat{\mathit{\text{NNT}}}={\left[\widehat{\alpha}\right]}^{1}$. The value of 1 indicates the largest possible effect of NNT, while in correspondence of no covariate effect (RD=0) the NNT value is ±∞. The largest possible harmful effect is −1. Positive and negative values of NNT represent the expected number of patients needed to be treated for one additional patient to benefit and to be harmed, respectively.
In the case of the loglog link, g=l o g(−l o g(•)), e x p(α), e x p(β) and e x p(γ) are the ratio between cumulative hazard functions associated with the change of one unit in the covariates. This ratio is equal to that of hazard functions, only in the proportional hazard case.
The method allows to estimate the measures of effect also for continuous covariates. For example, the evaluation of a biomarker effect measured on a continuous scale, without cutoffs, is still possible with this methodology.
The use of different link functions to obtain a particular measure of effect, is an established technique in binomial regression, where the use of noncanonical links, such as the logarithm, allows to obtain adjusted measures of impact different from the odds ratio, [1, 38]. Wacholder, [39], is an excellent reference for deepen such aspects in the framework of logistic regression.
In such a case, the estimated gtransformed survival probability differences change during followup time. In order to show the effect, varying in time, of a dichotomous covariate, for example treatment T, it is useful to adopt a graphical display, where the time is put on the horizontal axis while the function e x p(α+γ_{1}B_{1}(t)+γ_{2}B_{2}(t)) is on the vertical axis (exponentiation is not used with the identity link; R D(t)=α+γ_{1}B_{1}(t)+γ_{2}B_{2}(t)). In this case the estimated N N T(t) is naturally varying through followup time and again obtain by inversion: R D(t)^{−1}.
For a continuous covariate, such as Age, A in the example, it is possible to use a surface plot, where Age and time are on the x and y axis, while the z axis reports the covariate effect with respect to a reference value. It would also be possible to model Age effect with spline bases. In this case, the interaction between Age and time is obtained through tensor product spline bases of Age and time.
When a large number of pseudotimes is used, spline functions allow to model parsimoniously the baseline risk compared to indicator functions. This is particularly important for the modelling of timedependent effects in connection to the different link functions. In principle when a covariate effect is constant using a specific link, it should be timevarying with the other links. No statistical evidence against a constant covariate effect for more than one link may only be due to lack of power. The problem can also be exacerbated by some multiple testing issue. Timedependent effects selection depends therefore on the link transformation used. As a consequence, the adjusted effect of a covariate may be constant using a link function, but timedependent using a different link.
Moreover, the fitted values of the different models selected for the different link functions are generally different, being equal only if the models are saturated. Traditionally, the strategy used in the application of transformation models such as (3) was to select the best fitting g transform, i.e., the transform where covariate effects are constant through time, see [40, 41] as examples. The approach considered here is different. The interest is in using the gtransform which is the most informative for the clinical or biological counterpart. Generally the best fitting link function and the one selected by the researcher are not the same. Time dependent effects should therefore be expected in the model.
Pointwise confidence intervals
where $\mathrm{st.error}\left(\widehat{\alpha}\right)$ is the estimated clusterrobust standard error for the model parameter α. When g is the log, logit or l o g−l o g link function, the 95% CI for the treatment effect (respectively an RR, OR or HR) is [ e x p.(l_{ lower }),e x p(l_{ upper })]. With the identity link, the 95% confidence intervals is [ l_{ lower },l_{ upper }], without additional transformations, and the corresponding interval for NNT is [ 1/l_{ upper },1/l_{ lower }].
A Clarification is necessary for the confidence interval of the NNT.
where V(•) stands for the clusterrobust variance, while C o v(•,•) stands for the clusterrobust covariance of two random variables. When the variances at the different times are calculated, the pointwise 95% CI can be computed as before.
Software implementation
The approach to censored data regression based on pseudo values was applied to regression models for the cumulative incidence functions in competing risks and for multistate modeling [43], for the restricted mean [44] and for the survival function at a fixed point in time [45]. Implementation details and software can be found in Klein et al., [32], and Andersen and Perme, [46].
Software is available to compute pseudo values (macro %pseudosurv in SAS and function pseudosurv in R package pseudo[32]) Standard GEE tools available in SAS or R can be used for regression. In SAS the proc genmod allows to change link functions using the instructions FWDLINK and INVLINK. In R, the package geepack can be conveniently used, see [32] for details.
where the variable pseudo contains the pseudo values and the variable tpseudo the pseudotimes according to the software reported in [32]. The R function rcs of the package rms, [47], is used to compute restricted cubic spline bases. Each subject is represented by multiple rows in the data, one for each pseudo time. The records for each subject are identified by means of the variable id which is used to estimate the robust standard error by the geese function. Using the identity link function, the estimated coefficients can be interpreted as the adjusted R D(t) estimates.
Results
The estimation of R D(t) through pseudovalues is evaluated through a simple simulation study. Moreover, two real examples are presented. The first example concerns a prostate cancer trial, well known in competing risks literature, to show a situation where proportional hazard fail to model the treatment effect during the whole followup. The second example regards the analysis of prognostic factors in an observational study of primary retroperitoneal soft tissue sarcoma patients. R software, [48], was used for the simulation and both the examples presented.
Simulation
Event times generated according to a Coxexponential model with a confounder X and an exposure status Z
BIAS  $\sqrt{\mathit{\text{MSE}}}$  

PV  PV  Cox  PV  PV  Cox  
(Z td)  (Z const)  CGPM  (Z td)  (Z const)  CGPM  
200  0.0112  0.0966  0.0170  0.0756  0.1042  0.0580 
400  0.0047  0.0681  0.0145  0.0701  0.0785  0.0523 
600  0.0081  0.0034  0.0094  0.0528  0.0392  0.0376 
Width  Coverage  
PV  PV  Cox  PV  PV  Cox  
(Z td)  (Z const)  CGPM  (Z td)  (Z const)  CGPM  
200  0.2857  0.1527  0.2143  0.9430  0.3150  0.9480 
400  0.2791  0.1527  0.1926  0.9560  0.5830  0.9340 
600  0.1949  0.1527  0.1418  0.9350  0.9400  0.9320 
R D(t) estimated through the CGPM using the Cox proportional hazard regression model (the model used to generate the simulation data) is used as the benchmark estimation method.
The method based on pseudovalues with Z time dependent appears effective especially in terms of bias. Confidence interval coverage is good, although the width of the confidence intervals with pseudovalues is fairly large. It is interesting to observe the results of pseudovalues with the covariate Z not time dependent. In this case the estimated risk difference is constant through time, namely RD, a situation which can result from lack of power to detect the timedependence of the R D(t). The simulation results appear very interesting for late followup times. At time 600, results are very similar to that of the Cox proportional model. However, the PV method with identity link and without time dependence estimating a constant RD leads to a strong undercoverage demonstrating that the estimation of a constant RD may be misleading.
Event times generated according to a Coxexponential model with a confounder X and an exposure status Z with timedependent effect
BIAS  $\sqrt{\mathit{\text{MSE}}}$  

PV  Cox  Cox  PV  Cox  Cox  
r c s(t)  log(t)  r c s(t)  r c s(t)  log(t)  r c s(t)  
200  0.0290  0.0036  0.0151  0.0792  0.0689  0.1403 
400  0.0248  0.0108  0.0188  0.0730  0.0657  0.1119 
600  0.0185  0.0134  0.0242  0.0730  0.0657  0.1119 
Width  Coverage  
PV  Cox  Cox  PV  Cox  Cox  
r c s(t)  log(t)  r c s(t)  r c s(t)  log(t)  r c s(t)  
200  0.2899  0.2700  0.5162  0.9270  0.9530  0.9120 
400  0.2747  0.2528  0.4211  0.9470  0.9500  0.9170 
600  0.2158  0.1860  0.3256  0.9410  0.9410  0.9300 
In this simulation the Cox model with the correct specification of the time dependent effect of Z, that is log(t), is used as benchmark estimation. When the time dependence of Z is modelled using the restricted cubic spline, the performance of the CGPM is less appealing, compared to the benchmark, regarding all the parameters considered into the simulation. The pseudovalue model is really a competitor in this situation. It is in particular interesting to observe the 95% confidence interval width. The transformation model using pseudovalues with identity link is a valuable alternative to the CGPM when the time dependent effect in the Cox model is unknown and modelled using a flexible method. Model checking is therefore very important and pseudovalues can be of help also in this case, see the work of Anderson and Perme, [46].
Prostate cancer
Literature data on 502 prostate cancer patients, publicly available at the web site http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets, (Byar & Greene prostate cancer data), treated with different doses of diethylstilbestrol in a randomized clinical trial, [18], were used to estimate the adjusted treatment effect (high versus low dose) on overall mortality. Seven covariates were used for adjustment, namely: age (0, < 75 years; 1, 75−80 years; 2, ≥ 80 years), weight index (0, ≥ 100;1, 80−99; 2, < 80), performance rating (0, normal; 1, limitation of activity), history of cardiovascular disease (0, no; 1, yes), serum haemoglobin (0, ≥ 12 g/100 ml; 1, 912 g/100 ml; 2, < 9 g/100 ml), size of primary lesion (0, < 30 c m^{2}; 1, ≥ 30 c m^{2}), and Gleason stage category (0, ≤ 10; 1, > 10). 483 patients with complete information on the seven covariates available were considered. 344 patients died: 149 for cancer; 139 for cardiovascular causes; 56 for other causes.
The estimated R D(t) is based on a proportional hazard model. In fact there is no evidence for timedependent treatment effect in Cox model according to Schoenfeld residuals. However, Kay, [18], carefully investigated the fit of the Cox model, dividing the time axis into three time interval: [0−13]; (13−32]; (32−∞). The log HR for treatment has a positive sign in the first period (0.09), then becomes negative (−0.40 and −0.31). Comparing by likelihood ratio the model fitted on three intervals with an overall survival model, evidence is found against the proportional HR assumption. In fact, cardiovascular deaths are more frequent than cancer deaths earlier during followup while, later on, cancer deaths are prevalent. Accordingly, the beneficial effect of treatment appears evident only after the first year of followup. The estimates obtained with the CGPM appear therefore distorted due to the use of a proportional hazard model.
As a model with jumps in covariate effects at specific times is not biologically plausible, a Cox regression model with an interaction between the treatment and a function of time is used (specifically a Bspline with 1 interior knot at the median of the failure time distribution). According to the work of Hess, [29], splines are used to model flexibly the possible treatment time dependent effect. The estimated R D(t) according to this flexible model is reported in Figure 2. The harmful initial treatment effect is not yet captured.
A backward selection procedure is used to select the time dependent effects for each covariate. The complete model has a total of 45 degrees of freedom, including the 8 covariates and their time dependent effects. The selected model exhibit a time dependent effect of history of cardiovascular disease, size of primary lesion and Gleason stage. There is no evidence for a time dependent treatment effect. The constant estimated RD is 2.2% with 95% confidence interval [ −3.5%;7.9%]. The corresponding constant estimate of NNT is about 45 patients to be treated for one patient to benefit with 95% confidence interval [ −∞ to −28;13 to ∞].
R D ( t ) estimated with different methods at months 13, 32 and 60
RD  

13  32  60  
R D(t): Cox prop  3.7%; (−1.2%;7.7%)  5.7%; (−1.2%;12.3%)  5.8%; (−1.2%;12.7%) 
R D(t): Cox TD  0.1%; (−6.8%;6.5%)  6.5%; (−1.5%;14.1%)  7.7%; (−0.5%;15.6%) 
RD: PV identity Z const  2.2%; (3.5%;−7.9%)  2.2%; (3.5%;−7.9%)  2.2%; (3.5%;−7.9%) 
R D(t): PV identity Z TD  −1.6%; (−8.5%;5.3%)  4.6%; (−3.2%;12.4%)  6.2%; (−1.4%;13.7%) 
In the right panel of Figure 2 are also reported the averaged survival probabilities obtained by the CGPM using the Cox model and the transformation model both with the timedependent effect of treatment. The plot of the averaged survival probabilities is an important completion of the R D(t) plot. In fact R D(t), as usual for the effect measures, is reducing two numbers to a single number.
The prostate cancer data outline a scenario in which the proportional hazard assumption for the treatment effect is not tenable during all followup times. Data on prostate cancer should be actually analyzed accounting for competing risks. When a non competing risks survival analysis is performed, CGPM applied to the Cox model without a time dependent treatment effect gives an estimate of the R D(t) increasing during followup until reaching a plateau. At the same time, the constant estimate RD obtained using PV provides a distorted estimate constant through followup. CGPM applied to the Cox model with a time dependent treatment effect provides an R D(t) estimate not yet capturing the initial harmful treatment effect. The timedependent estimate obtained from the pseudovalue model is instead effective in describing the treatment effect during followup time: harmful at the beginning, when cardiovascular deaths are more frequent, beneficial later on when cancer deaths are more frequent
Primary retroperitoneal soft tissue sarcoma
Discussion
In survival analysis the adjusted measure of association everywhere adopted is the hazard ratio. Although the efficiency of the hazard ratio makes it attractive for hypothesis testing, it may not carry the most useful information for clinicians/biostatisticians. Schechtman, [11], suggests using also absolute measures in conjunction with relative measures of covariateoutcome association. To provide adjusted measures of association different from the hazard ratio, a simple strategy is going through the calculation of the predicted probabilities of event for an “average” subject, the so called “average covariate method”. Such a procedure estimates the measures of effect for an hypothetical average subject and not population averaged estimates. An alternative idea to provide adjusted summary measures of effect is the corrected group prognosis method [27, 28]. Extending this idea, using the concept of counterfactuals, Laubender and Bender [13] proposed methods for computing such population measures taking into account the confounder distribution. Bootstrap is then used to obtain pointwise confidence intervals. Such approaches are particularly appealing as they may adopt the Cox regression model which is widely accepted in medical literature.
However the Cox model may easily not be the best regression procedure to be applied, simply because of the assumption of proportional hazards. In fact, in presence of time dependent effects Cox regression may be less appealing. This is demonstrated here through a simple simulation. When timedependent effects in the Cox model are specified using wellknown flexible methods, [29], without committing to a specific functional form, the estimates are not optimal, especially in terms of efficiency. In these circumstances the results of the simple simulation presented here, suggest that the use of the pseudovalue model may represent a valid alternative to the CGPM. The simulation is not exhaustive and more work is needed to fully understand the properties and the relationships among the different estimation methods.
Considering a similar problem in the context of logistic regression, Gehrmann and colleagues, [49], concluded that the CGPM applied to logistic regression is the preferred method to estimate RD and NNT adjusted for covariates compared to binomial, Poisson and linear regression methods that directly estimates the RD (similar to pseudovalues with identity link) even if the fitted response function differs from the true response function. The context of timetoevent outcomes is more complex than that of logistic regression especially for the problem of time dependent effects. Whether similar results hold for the Cox model has therefore to be further explored thorough a series of simulation studies. In any case, when using the CGPM, the basic model used to obtain event probabilities during followup has to be adequate. This means, for example, that the proportional hazard Cox regression should not be applied if the proportional hazard assumption is not satisfied.
In clinical literature, results of statistical analysis are commonly reported in terms of regression coefficients and their confidence intervals. Applied survival analysis resorts entirely to the Cox model which is a particular case of transformation models. Transformation models include also the accelerated failure time models providing a variety of measures of effects to be considered.
It is to be noted that, additive and multiplicativeadditive hazard regression models, [31] are not comprised in this class. The estimated coefficients are differences of hazard rates rather than ratios. These models were mainly proposed to improve the fitting where the proportional hazards model is not adequate or to check for proportional hazard assumptions. Moreover, the measures of impact provided are still based on hazards.
In this work a simple approach to obtain point and interval estimates of association measures, by using transformation models through suitable link functions, is presented. The general technique of estimation based on pseudovalues proposed by Andersen and colleagues [43] is used as it is simply implementable with standard software.
Other techniques could have been considered to estimate the transformation models. In this context, it is of particular interest the estimate of the baseline survival function to model time dependent effects. Therefore, semiparametric techniques (see for example [50]) are not of interest here. Maximum likelihood (ML) estimation can however be conveniently used. Royston and Parmar proposed ML for transformation models with cloglog and logit links, [51], which are, however, the links of less interest here. Moreover as the pseudovalue is defined between the first and the last pseudo times (which are not the first and last event times) it still makes sense to have a constant RD model. In fact considering the whole followup time interval, and especially the beginning of the followup, R D(t) estimates should instead always be timedependent.
From the methodological viewpoint, the only difference introduced in the presented examples with respect to standard applications of the pseudovalue model is the use of spline functions to estimate the transformed baseline survival. This modification is without practical efforts, considered the wide availability of software to compute spline bases. From the modeling point of view, care must be paid to the monotonicity of the estimated transformed baseline survival function. In general, provided the number of knots is limited, no problems of non monotonicity were observed.
Another issue concerns the possible simultaneous use of different association measures estimated through different link functions in a transformation model. In such a case, practitioners must be aware that, due to lack of power, not all time dependent effects can be correctly specified, and likely the different models cannot hold simultaneously. In such a case, if different measures are of interest, different models can be used simultaneously only if the results are in agreement with each other.
From a theoretical point of view, cloglog and logit links guarantee that the estimated probabilities are within the range 01 but this is not guaranteed if log or identity links are used with standard software. Research is in progress to face this issue.
Conclusions
The use of different link functions in transformation models has been studied by several authors simply investigating the goodness of fit of the different links [40, 41].
The alternative perspective considered here, evaluates the use of different links with the goal of providing suitable measures of association between covariates and outcome. As a consequence, when a specific link function is chosen, it should be expected the need to include time dependent covariate effects into the model.
Transformation models estimated through pseudovalues appear an easily implemented alternative to the available approaches mainly based on Cox proportional hazard model to obtain adjusted measures of association eventually timedependent also for continuous covariates.
Declarations
Acknowledgments
The study was partly supported by the AIRC grant IG201213420. The Authors wish to thank Dr. A. Gronchi for providing data about retroperitoneal soft tissue sarcoma patients.
Authors’ Affiliations
References
 Greenland S: Modelbased estimation of relative risks and other epidemiologic measures in studies of common outcomes and in casecontrol studies. Am J Epidemiol. 2004, 160 (4): 301305.View ArticlePubMed
 McNutt LA, Wu C, Xue X, Hafner JP: Estimating the relative risk in cohort studies and clinical trials of common outcomes. Am J Epidemiol. 2003, 157 (10): 940943.View ArticlePubMed
 Marubini E, Valsecchi MG: Analysing Survival Data from Clinical Trials and Observational Studies. Statistics in Practice. 2004, Chichester: Wiley
 Cox D: Regression models and lifetables (with discussion). J R Stat Soc B. 1972, 34: 557565.
 Cox D: Partial likelihood. Biometrika. 1975, 62: 269276.View Article
 Laupacis A, Sackett DL, Roberts RS: An assessment of clinically useful measures of the consequences of treatment. N Engl J Med. 1988, 318 (26): 17281733.View ArticlePubMed
 Boracchi P, Mezzanotte G, Mariani L, Valagussa P, Marubini E: Clinically useful measures to assess the effectiveness of treatments: hints for a proper choice with special emphasis on cancer research. Tumori. 1990, 76 (1): 19.PubMed
 Jaeschke R, Guyatt G, Shannon H, Walter S, Cook D, Heddle N: Basic statistics for clinicians: 3. Assessing the effects of treatment: measures of association. CMAJ. 1995, 152 (3): 351357.PubMed CentralPubMed
 DiCenso A: Clinically useful measures of the effects of treatment. Evid Based Nurs. 2001, 4 (2): 3639.View ArticlePubMed
 Case LD, Kimmick G, Paskett ED, Lohman K, Tucker R: Interpreting measures of treatment effect in cancer clinical trials. Oncologist. 2002, 7 (3): 181187.View ArticlePubMed
 Schechtman E: Odds ratio, relative risk, absolute risk reduction, and the number needed to TreatWhich of these should we use?. Value Health. 2002, 5 (5): 431436.View ArticlePubMed
 Austin PC, Laupacis A: A tutorial on methods to estimating clinically and policymeaningful measures of treatment effects in prospective observational studies: a review. Int J Biostat. 2011, 7 (1): 132.View Article
 Laubender RP, Bender R: Estimating adjusted risk difference (RD, and number needed to treat (NNT, measures in the Cox regression model. Stat Med. 2010, 29 (7–8): 851859.View ArticlePubMed
 Austin PC: Absolute risk reductions and numbers needed to treat can be obtained from adjusted survival models for timetoevent outcomes. J Clin Epidemiol. 2010, 63 (1): 4655.View ArticlePubMed
 Bender R, Kuss O: Methods to calculate relative risks, risk differences, and numbers needed to treat from logistic regression. J Clin Epidemiol. 2010, 63 (1): 78.View ArticlePubMed
 Bender R, Kuss O, Hildebrandt M, Gehrmann U: Estimating adjusted nnt measures in logistic regression analysis. Stat Med. 2007, 26 (30): 55865595.View ArticlePubMed
 Austin PC: Different measures of treatment effect for different research questions. J Clin Epidemiol. 2010, 63 (1): 910.View ArticlePubMed
 Kay R: Treatment effects in competingrisks analysis of prostate cancer data. Biometrics. 1986, 42 (1): 203211.View ArticlePubMed
 Ardoino I, Miceli R, Berselli M, Mariani L, Biganzoli E, Fiore M, Collini P, Stacchiotti S, Casali PG, Gronchi A: Histologyspecific nomogram for primary retroperitoneal soft tissue sarcoma. Cancer. 2010, 116 (10): 24292436.PubMed
 Hauck WW, Anderson S, Marcus SM: Should we adjust for covariates in nonlinear regression analyses of randomized trials?. Control Clin Trials. 1998, 19: 249256.View ArticlePubMed
 Senn SJ: Covariate imbalance and random allocation in clinical trials. Stat Med. 1989, 8: 467475.View ArticlePubMed
 International Conference on Harmonization: Harmonised tripartite guideline: statistical principles for clinical trials. Stat Med. 1999, 18 (15): 19051942. Cited By (since 1996) 214,
 Fleming TH, Harrington DP: Nonparametric estimation of the survival distribution in censored data. Comm Stat. 1984, 13: 24692486.View Article
 Kalbfleisch JD, Prentice RL: The Statistical Analysis of Failure Time Data (Wiley Series in Probability and Statistics). 2002, New York: WileyInterscienceView Article
 Ghali WA, Quan H, Brant R, van Melle G, Norris CM, Faris PD, Galbraith PD, Knudtson ML: Comparison of 2 methods for calculating adjusted survival curves from proportional hazards models. JAMA. 2001, 286 (12): 14941497.View ArticlePubMed
 Nieto FJ, Coresh J: Adjusting survival curves for confounders: a review and a new method. Am J Epidemiol. 1996, 143 (10): 10591068.View ArticlePubMed
 Chang IM, Gelman R, Pagano M: Corrected group prognostic curves and summary statistics. J Chronic Dis. 1982, 35 (8): 669674.View ArticlePubMed
 Makuch RW: Adjusted survival curve estimation using covariates. J Chronic Dis. 1982, 35 (6): 437443.View ArticlePubMed
 Hess KR: Assessing timebycovariate interactions in proportional hazards regression models using cubic spline functions. Stat Med. 1994, 13 (10): 10451062.View ArticlePubMed
 Therneau TM, Grambsch PM: Modeling Survival Data: Extending the Cox Model. 2010, New York: Springer
 Martinussen T, Scheike TH: Dynamic Regression Models for Survival Data (Statistics for Biology and Health). 2006, New York: Springer
 Klein JP, Gerster M, Andersen PK, Tarima S, Perme MP: SAS and R functions to compute pseudovalues for censored data regression. Comput Methods Programs Biomed. 2008, 89 (3): 289300.PubMed CentralView ArticlePubMed
 Perme MP, Andersen PK: Checking hazard regression models using pseudoobservations. Stat Med. 2008, 27 (25): 53095328.PubMed CentralView ArticlePubMed
 Liang KY, Zeger SL: Longitudinal data analysis using linear models. Biometrika. 1986, 73: 1322.View Article
 Royston P, Parmar MK: Flexible parametric proportionalhazards and proportionalodds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Stat Med. 2002, 21 (15): 21752197.View ArticlePubMed
 Harrell FE: Regression Modeling Strategies. 2001,View Article
 Pan W: Akaike’s information criterion in generalized estimating equations. Biometrics. 2001, 57: 120125.View ArticlePubMed
 Hardin JW, Hilbe J: Generalized Linear Models and Extensions. 2001, College Station, Texas: Stata Press
 Wacholder S: Binomial regression in GLIM: estimating risk ratios and risk differences. Am J Epidemiol. 1986, 123 (1): 174184.PubMed
 Bennett S: Analysis of survival data by the proportional odds model. Stat Med. 1983, 2 (2): 273277.View ArticlePubMed
 Wei LJ: Testing goodness of fit for proportional hazards model with censored observations. J Am Stat Assoc. 1984, 79 (387): 649652.View Article
 Altman DG: Confidence intervals for the number needed to treat. BMJ. 1998, 317 (7168): 13091312.PubMed CentralView ArticlePubMed
 Andersen PK, Klein JP, Rosthöj S: Generalised linear models for correlated pseudoobservations, with applications to multistate models. Biometrika. 2003, 90 (1): 1527.View Article
 Andersen PK, Hansen MG, Klein JP: Regression analysis of restricted mean survival time based on pseudoobservations. Lifetime Data Anal. 2004, 10 (4): 335350.View ArticlePubMed
 Klein JP, Logan B, Harhoff M, Andersen PK: Analyzing survival curves at a fixed point in time. Stat Med. 2007, 26 (24): 45054519.View ArticlePubMed
 Andersen PK, Perme MP: Pseudoobservations in survival analysis. Stat Methods Med Res. 2010, 19 (1): 7199.View ArticlePubMed
 Jr FEH: Rms: Regression modeling strategies. 2013, R package version 3.63. [http://CRAN.Rproject.org/package=rms],
 R Core Team: R: A Language and Environment for Statistical Computing. 2012, Vienna: R Foundation for Statistical Computing, [http://www.Rproject.org/]
 Gehrmann U, Kuss O, Wellmann J, Bender R: Logistic regression was preferred to estimate risk differences and numbers needed to be exposed adjusted for covariates. J Clin Epidemiol. 2010, 63 (11): 12231231.View ArticlePubMed
 Cheng S, Wei L, Ying Z: Analysis of transformation models with censored data. Biometrika. 1995, 82 (4): 835845.View Article
 Royston P, Parmar MK: Flexible parametric proportionalhazards and proportionalodds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Stat Med. 2002, 21 (15): 21752197.View ArticlePubMed
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14712288/14/97/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.