On the censored costeffectiveness analysis using copula information
 Charles Fontaine^{1}Email authorView ORCID ID profile,
 JeanPierre Daurès^{1} and
 Paul Landais^{1}
DOI: 10.1186/s1287401703059
© The Author(s) 2017
Received: 20 June 2016
Accepted: 2 February 2017
Published: 15 February 2017
Abstract
Background
Information and theory beyond copula concepts are essential to understand the dependence relationship between several marginal covariates distributions. In a therapeutic trial data scheme, most of the time, censoring occurs. That could lead to a biased interpretation of the dependence relationship between marginal distributions. Furthermore, it could result in a biased inference of the joint probability distribution function. A particular case is the costeffectiveness analysis (CEA), which has shown its utility in many medicoeconomic studies and where censoring often occurs.
Methods
This paper discusses a copulabased modeling of the joint density and an estimation method of the costs, and quality adjusted life years (QALY) in a costeffectiveness analysis in case of censoring. This method is not based on any linearity assumption on the inferred variables, but on a punctual estimation obtained from the marginal distributions together with their dependence link.
Results
Our results show that the proposed methodology keeps only the bias resulting statistical inference and don’t have anymore a bias based on a unverified linearity assumption. An acupuncture study for chronic headache in primary care was used to show the applicability of the method and the obtained ICER keeps in the confidence interval of the standard regression methodology.
Conclusion
For the costeffectiveness literature, such a technique without any linearity assumption is a progress since it does not need the specification of a global linear regression model. Hence, the estimation of the a marginal distributions for each therapeutic arm, the concordance measures between these populations and the right copulas families is now sufficient to process to the whole CEA.
Keywords
Costeffectiveness analysis Censored data Copulas Parametric models Subgroups analysisBackground
Due to the variety of treatments for a specific health problem and in conjunction with their increasing costs, costeffectiveness studies of new therapies is challenging. These studies could achieve to a statistical analysis since that the common practice in laboratories is to collect individual patient cost data in randomized studies. Furthermore, it is now possible to compute the incremental net benefit from the use of a new therapy in comparison with the commoninuse therapy.
In the last decades, the costeffectiveness analysis (CEA) of new treatments became an actual subject of work for statisticians. It is used in two particular designs: decision modelingbased CEA and trialbased CEA. The major difference between both approaches is that in the case of trialbased CEA, data are gathered at the patients level in particular studies and it may lead to overlearning from the study, which may lead to mistakes in interpretation when the results are inferred for populations. In contrast, decision modelingbased CEA are based on easily generalizable data.
where λ is the willingnesstopay for a unit of effectiveness.
In the literature, many articles propose ways to estimate these quantities. At first, Willan and Lin [1] proposed an approach which is based on the sample mean. It was applied directly to survival time years. In case of censored data, they proposed to estimate the survival function for each arm using the productlimit estimator [2] and then to estimate the survival time by integrating the survival functions until a time τ (i.e. μ _{ j } is the life expectancy until τ), the maximal time of followup. Concerning the estimation of costs in case of censoring, many estimators may be used. Zhao et al. [3] have shown the equivalence among them.
More recently, Willan et al. [6] proposed to realize the whole costeffectiveness analysis using linear regression methods. Let C _{ i } be the observed cost for patient i, then \(\mathbb {E}(C_{ji})=\beta _{C_{i}}^{T}Z_{C_{ji}}, i=1,2,\ldots,n_{j}, Z_{C_{ji}}\) is a vector of covariates affecting costs and \(\beta _{C_{j}}\) is a vector of unknown regression parameters. Then, using the inverse probability of censoring weighting (IPCW) methods, they proposed a way to estimate the second component of \(\beta _{C_{j}}\) which is the mean difference in costs between randomization groups adjusted to other covariates, \(\hat {\Delta }_{c}\), and its associate variance. The same methodology is done there to estimate the mean difference in mean survival between randomization groups adjusted for the other covariates, \(\hat {\Delta }_{e}\). Thus, they proposed to estimate the ICER adjusted for the quality of life by \(\hat {\Delta }_{c}/\hat {\Delta }_{e}\) and used this time again the Fieller’s theorem to find a 100(1−α) confidence interval. For the adjusted INB, they proposed to estimate it by \(\hat {b}_{\lambda }=\lambda \hat {\Delta }_{e}\hat {\Delta }_{c}\) with variance \(\hat {\sigma }_{\lambda }^{2}=\lambda ^{2}\hat {\sigma }^{2}_{\Delta _{e}}+\hat {\sigma }^{2}_{\Delta _{c}}2\lambda \hat {\sigma }_{\Delta _{c}\Delta _{e}}\). Thus, if the INB is positive, the therapy is costeffective and if it is negative, there is no costeffectiveness. One remarks that it is possible to determine the costeffectiveness under a significance level α using a statistic of test: if \(\hat {b}_{\lambda }/\hat {\sigma }_{\lambda }\) is greater than the test level z _{1−α }, at a level α, the therapy is costeffective compared to the standard.
From that linear regression approach, many variants exist either parametric or semiparametric [7]. The main problem of the estimator for costeffectiveness analysis based on the linear regression is that even if the inverseprobability censored weighting estimator is consistent, it is not efficient because if an individual is censored before or at time a _{ K+1}, the patient does not contribute to the sum that constitutes this estimator and the practitioner may lose some statistical significant information. Furthermore, when there is no censoring, it is equivalent to the solve ordinary least squares equations; for which the linearity assumption may lead to serious bias in case of nonlinearity. In the case of a semiparametric estimator, the main problem is that the estimation of the conditional moments may be hard to find accurately or may lead to some bias in function of the used estimator. For these reasons, we introduce a new costeffectiveness analysis methodology and a modeling of the joint density function between costs and QALY using parametric copulas. It is therefore based only on the dependence between covariates and the prior information on the variables distributions.
Methods
Copula function
Note that parametric and nonparametric copulas estimation exists. Also, one has to remark that given two marginal distributions, the copula that joins them is unique if and only if these margins are continuous. In the literature, few papers deal with copula and costs data. At first, Embrechts et al. [9] present an original article about the correlation measurement with those data. However, the copula concepts presented there are more about generalizations of the copula theory than about a modeling methodology. Secondly, Frees and Valdez [10] use copulas for costs data in the insurance area. Their work explains how to fit copulas and measure the mortality risks. Thirdly, there is Hougaard [11] who uses copula with costs data in a multivariate survival analysis context. Also, there is Zhao and Zhou who work with copula models using medical costs, but in a stochastic process context, which is not adapted to QALY and costs data when the information is limited.
In this paper, we assume the copula to be parametric, which means that the copula can be written in a particular way as a function of the chosen family, with one parameter who summarizes the dependence between the variables, and the marginal distributions to be continuous. For more about copulas, see Nelsen [12].
Model
Determination of QALY in terms of time and quality of life
 1.
\(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) and η _{ ji } are dependent,
 2.
\(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) and \(\phantom {\dot {i}\!}\eta _{{adj}_{ji}}\) are independent,
 3.
C _{ ji } and η _{ ji } are dependent.
Proofs are provided in the “Appendix” section.
Estimation of the parameters of the distributions
Determination of the parametric distributions
 1.
\(C_{j} \sim Normal \left (\mu _{C_{j}}, \sigma ^{2}_{C_{j}}\right),\)
 2.
\(\phantom {\dot {i}\!}C_{j} \sim {Gamma} (\mu _{C_{j}}, \rho _{C_{j}}),\)
 3.
\(C_{j} \sim Lognormal \left (\nu _{C_{j}},\tau ^{2}_{C_{j}}\right),\)
where ν _{ C } and \(\tau ^{2}_{C}\) are mean and variance of the logcosts, i.e. \(\nu _{C}=2log(\mu _{C})\frac {1}{2}log\left (\sigma ^{2}_{C}+\mu _{C}^{2}\right)\) and \(\tau _{C}^{2}=log\left (\sigma _{C}^{2}+\mu _{C}^{2}\right)2log(\mu _{C})\). Furthermore, ρ _{ C } is the shape parameter of the Gamma distribution, which is such that \(\rho _{C}=\mu _{C}^{2}/\sigma ^{2}_{C}\). Thus, once each modeling is achieved, a selection of the better parametric distribution has to be done using the deviance criteria. The best fit to the data corresponds to the distribution that has the lower deviance, which is minus two times the loglikelihood.
 1.
\(T_{{adj}_{j}} \sim Normal\left (\mu _{T_{adj}}, \sigma ^{2}_{T_{adj}}\right),\)
 2.
\(\phantom {\dot {i}\!}T_{{adj}_{j}} \sim Gamma (\mu _{T_{adj}}, \rho _{T_{adj}}).\)
Inference on Kendall’s tau
To further compute the copula parameter for each tested copula in the selection process, one has to get the global dependence parameter: Kendall’s tau. The idea is that performing inference on Kendall’s tau instead of directly on the copula parameter for each tested copula leads to only one estimation process instead of as many inferences as the quantity of tested models. Surely, the maximum likelihood estimator can be used in order to get the copula parameter. However, in this paper, we suggest the use of Kendall’s tau for nonrandomized data in reason of the easiness of the method. Furthermore, as shown by Genest et al. [18, 19], the difference between both inference methods is slightly significant compared to the gain in computational efficiency. Note that in case of randomized data, the exchangeability phenomenon among individuals occurs and concordance measure may be biased. Hence, in this case, we suggest the use of a standard maximum likelihood estimation for the copula parameter.
Bayesian selection of the copula
They proposed to use the Monte Carlo approximations to \(\mathbb {E}[\!D(\Phi _{k})\boldsymbol {x},\mathcal {M}_{k}]\) and \(\mathbb {E}[\!\Phi _{k}\boldsymbol {x},\mathcal {M}_{k}]\) which are respectively \(L^{1}\sum _{l=1}^{L} D\left (\Phi _{k}^{l}\right)\) and \(L^{1}\sum _{l=1}^{L} \Phi _{k}^{l}\). Here, one supposes that \(\left \{ \Phi ^{(1)}_{k},\ldots,\Phi _{k}^{(L)} \right \}\) is a sample from the posterior distribution \(f(\Phi _{k}\boldsymbol {x},\mathcal {M}_{k})\). Then, one chooses the copula model in all the chosen range with the smaller DIC.
One remarks that such a selection process requests parametric copulas with a limited quantity of parameters. Hence, one suggests the use of archimedean and elliptic copulas families for these modelings. In fact, the range of dependance structure models provided by archimedean copulas is enough large to cover the dependance parametrisation in a costeffectiveness analysis. For example, a Clayton copula may represent a study where a small QALY is directly linked to small costs, but a huge QALY is independent of costs. Also, one remarks that if costs and QALY are independents, then the independence copula is selected and one gets directly the values that constitute the ICER.
Incremental costeffectiveness ratio
In this formula, z _{1−α/2} represents the 100(1−α/2) percentile of the standard normal distribution. Furthermore, \(\hat {\sigma }^{2}_{\Delta _{T_{adj}}}\) represents the variance of the distribution for effectiveness where \(\Delta _{T_{adj}}\) is the difference between \(\phantom {\dot {i}\!}T_{{adj}_{j=1}}\) and \(\phantom {\dot {i}\!}T_{{adj}_{j=0}}\). The same scheme arises for \(\hat {\sigma }^{2}_{\Delta _{C}}\). For \(\hat {\sigma }_{\Delta _{C} \Delta _{T_{adj}}}\), it is nothing more than the estimated covariance of the differences, between j=0 and j=1, in costs and in quality adjusted life years.
The reason to use Fieller’s theorem is to avoid standard ways of using bootstrap. As shown by Siani and Moatti [31], Fieller’s method is often as robust as bootstrap methods are (both parametric and nonparametric ways), even in problematic cases. However, for noncommon situations, an approach based on the ICER graphical plan (as proposed by Bang and Zhao [32]) is recommended. Indeed, in such a case, a graphical approach minimize the use of biased or nonsufficient statistics in a parametric confidence interval model.
Incremental net benefit
The adjusted INB(λ) is estimated by \(\widehat {INB}=\lambda (\widehat {\mathbb {E}(T_{{adj}_{j=1}})}\widehat {\mathbb {E}(T_{{adj}_{j=0}})})(\widehat {\mathbb {E}(C_{j=1})}\widehat {\mathbb {E}(C_{j=0})})\) with variance \(\hat {\sigma }^{2}_{\lambda }=\lambda ^{2}\hat {\sigma }^{2}_{\Delta _{T_{adj}}}+\hat {\sigma }^{2}_{\Delta _{C}}2\lambda \hat {\sigma }_{\Delta _{C} \Delta _{T_{adj}}}\) where λ is the willingnesstopay for a unit of effectiveness.
Subgroups analysis
It is possible to accomplish a cohort analysis using that procedure. The main idea here is to perform a costeffectiveness analysis while achieving a discrimination between two or more subgroups. The principle is that there exists a baseline variable Z _{ k },k∈{1,2,…,d}, even for costs than for QALY, which is in fact a categorical variable (dichotomous or multichotomous) and for which one should determine a marginal INB. Such subgroups have to be based on clinically important covariates. Since that these subgroups are in the therapeutic arms, it is not possible to assume that they are balanced. As shown in Nixon and Thompson [33], and Tsai and Peace [34], a naive use of these subgroups information without any adjustment may lead to serious bias.
Let \(Z_{jki}^{C}, Z_{jki}^{T_{adj}}\) be some population attribute indicator covariates (e.g. sex, smoking status, etc.) for costs and QALY on an individual i belonging to the clinical arm j. For the sake of illustrating the concept here, let say that one tests the therapeutic effect on smokers. Therefore, there will be four subgroups: smokers in the treated group, nonsmokers in the treated group, and the same for the control groups.
Here, the terms \(cov(\overline {T}_{{adj}_{j=1}}, \overline {C}_{j=1})\) and \(cov(\overline {T}_{{adj}_{j=0}}, \overline {C}_{j=0})\) are computed in a standard way similarly to the randomized case. However, for the crossedarms covariance terms, \(cov(\overline {T}_{{adj}_{j=1}}, \overline {C}_{j=0})\) and \(cov(\overline {T}_{{adj}_{j=0}}, \overline {C}_{j=1})\), the approach that we suggest is to estimate the cumulative joint distributions \(F(\overline {E}_{j=1}, \overline {C}_{j=0})\) by \(C_{\hat {\theta }}(\hat {F}_{\bar {E}}(y_{j=1}),\hat {F}_{\bar {C}}(\chi _{j=0}))\) and \(F(\overline {T}_{{adj}_{j=1}}, \overline {C}_{j=0})\) by \(C_{\hat {\theta }}(\hat {F}_{\bar {T_{adj}}}(y_{j=1}),\hat {F}_{\bar {C}}(\chi _{j=0}))\) and F \((\overline {T}_{{adj}_{j=1}}, \overline {C}_{j=0})\) by \(C_{\hat {\theta }}(\hat {F}_{\bar {T_{adj}}}(y_{j=1}),\hat {F}_{\bar {C}}(\chi _{j=0}))\) according to the methodology shown in this paper, and then using the covariance definition \(cov(\overline {C},\overline {T_{adj}})=\mathbb {E}[\!\overline {T_{adj}},\overline {C}] \mathbb {E}[\!\overline {T_{adj}}]\mathbb {E}[\overline {C}]\), compute the desired covariance from the estimated joint density function and the estimated marginal density functions.
Results and discussion
This section gives an illustration of the performance of the copula which provides the best estimate of the true copula and the cumulative joint distribution of costs and QALY respectively, according to the method presented above. Let the exact copula be \(C_{\theta }^{(\star)}(F_{T_{adj}}(y\Phi _{T_{adj}}),F_{C}(\chi \Phi _{C}))\) and its estimate be \(C_{\hat {\theta }}^{(i)}(\tilde {F}_{T_{adj}}(y\hat {\Phi }_{T_{adj}}),\tilde {F}_{C}(\chi \hat {\Phi }_{C}))\) where (i) is the copula model selected among all the tested models and (⋆) the exact copula model. Furthermore, \(\tilde {F}\) is the chosen distribution for F. Then, the objective of these simulations is to show that the bias generated by the approximation of θ by \(\hat {\theta }, \Phi =\Phi _{C} \cup \Phi _{T_{adj}}\) by \(\hat {\Phi }\), the selection of C ^{(i)} as the copula model and \(\tilde {F}_{C}, \tilde {F}_{T_{adj}}\) as the marginal parametric models, is relatively weak.
Scheme of the 27 data generated cases
Generating copula  Costs distribution  Censoring level  DGP 

Gaussian copula  \(F_{C} \sim \mathcal {N}(\mu _{C}=1500, \sigma _{C}=400)\)  15%  DGP 1 
θ≈0.809  30%  DGP 2  
70%  DGP 3  
F _{ C }∼Γ(shape _{ C }=12,scale _{ C }=125)  15%  DGP 4  
30%  DGP 5  
70%  DGP 6  
\(F_{C} \sim log\mathcal {N}(\nu _{C}=7.30,\tau _{C}=0.25)\)  15%  DGP 7  
30%  DGP 8  
70%  DGP 9  
Clayton copula  \(F_{C} \sim \mathcal {N}(\mu _{C}=1500, \sigma _{C}=400)\)  15%  DGP 10 
θ=3  30%  DGP 11  
70%  DGP 12  
F _{ C }∼Γ(shape _{ C }=12,scale _{ C }=125)  15%  DGP 13  
30%  DGP 14  
70%  DGP 15  
\(F_{C} \sim log\mathcal {N}(\nu _{C}=7.30,\tau _{C}=0.25)\)  15%  DGP 16  
30%  DGP 17  
70%  DGP 18  
Gumbel copula  \(F_{C} \sim \mathcal {N}(\mu _{C}=1500, \sigma _{C}=400)\)  15%  DGP 19 
θ≈0.809  30%  DGP 20  
70%  DGP 21  
F _{ C }∼Γ(shape _{ C }=12,scale _{ C }=125)  15%  DGP 22  
30%  DGP 23  
70%  DGP 24  
\(F_{C} \sim log\mathcal {N}(\nu _{C}=7.30,\tau _{C}=0.25)\)  15%  DGP 25  
30%  DGP 26  
70%  DGP 27 
The simulation of linearly dependent and uncensored covariates for costs leads to a bias in our advantage for the computation of the mean and the variance, compared to the estimation performed according to a standard clinical scheme. We decided therefore to challenge our method using the KaplanMeier mean estimate of the survival function and its associated variance (which is the appropriate approach in absence of covariates of interest) instead of the presented procedure based on the covariates. Then, we applied the following steps: selecting a parametric distribution for the margins, selecting a parametric copula using the information criterion and finally, looking for the copula parameter. We replicated this procedure 500 times for n=1000 data, then we collected the provided information on the frequency of successful procedures for the inference on the margins, the estimated Kendall’s tau and the choice of copula, respectively.
Inference on Kendall’s tau
Information on the estimation of the Kendall’s tau for each censoring level
Censoring level  Mean \((\hat {\tau }_{K})\)  Var \((\hat {\tau }_{K})\)  Min \((\hat {\tau }_{K})\)  Max \((\hat {\tau }_{K})\) 

Censoring =0%  0.6002  0.00019  0.5488  0.6476 
Censoring =15%  0.6011  0.00024  0.5416  0.6539 
Censoring =30%  0.6030  0.00035  0.5319  0.6648 
Censoring =70%  0.6089  0.00146  0.4624  0.7257 
Inference on the marginal distribution for the costs
Inference on the copula models
Frequency of choice of a copula for each DGP given 500 iterations for the three main copulas
DGP  Gaussian  Clayton  Gumbel 

copula  copula  copula  
DPG 1  460  0  40 
DGP 2  454  0  46 
DGP 3  25  94  381 
DGP 4  415  0  85 
DGP 5  429  0  71 
DGP 6  140  21  339 
DGP 7  473  0  27 
DGP 8  490  0  10 
DGP 9  261  10  229 
DGP 10  0  500  0 
DGP 11  0  500  0 
DGP 12  0  499  1 
DGP 13  18  482  0 
DGP 14  11  489  0 
DGP 15  68  371  61 
DGP 16  2  498  0 
DGP 17  6  493  1 
DGP 18  177  531  62 
DGP 19  2  0  498 
DGP 20  9  0  491 
DGP 21  1  9  490 
DGP 22  123  0  377 
DGP 23  109  0  391 
DGP 24  36  0  467 
DGP 25  80  0  420 
DGP 26  84  0  416 
DGP 27  121  0  379 
Frequency of choice of copula for each DGP on 500 iterations beyond the three main copulas and intermediate copulas
DGP  Gaussian  Student  Clayton  Gumbel  Frank  Joe 

copula  copula  copula  copula  copula  copula  
DPG 1  7  45  6  25  415  2 
DGP 2  2  44  4  25  425  1 
DGP 3  4  54  4  17  421  0 
DGP 4  32  297  0  56  112  3 
DGP 5  26  291  0  59  120  4 
DGP 6  23  284  0  79  113  1 
DGP 7  47  343  1  45  38  26 
DGP 8  55  344  0  43  33  25 
DGP 9  55  328  0  56  37  24 
DGP 10  0  12  364  0  124  0 
DGP 11  0  17  372  0  111  0 
DGP 12  0  5  385  0  110  0 
DGP 13  0  427  56  0  17  0 
DGP 14  2  433  55  0  10  0 
DGP 15  1  431  58  0  10  0 
DGP 16  1  477  19  2  1  0 
DGP 17  2  476  18  0  4  0 
DGP 18  4  472  19  0  4  1 
DGP 19  0  18  0  110  196  176 
DGP 20  0  20  0  107  214  159 
DGP 21  0  17  0  109  202  172 
DGP 22  13  236  0  153  67  31 
DGP 23  19  231  0  157  58  35 
DGP 24  23  161  0  204  75  37 
DGP 25  32  198  0  177  13  80 
DGP 26  39  213  0  148  13  87 
DGP 27  37  210  0  167  17  69 
According to the structure of the marginal distributions, the choice of the copula could only be done in a limited spectrum of families. Therefore, when one seeks to find the best structural copula family, it appears essential to include the most known copula families covering as many dependence states as possible. Thus, in harmony with Table 4, a selected copula which is not the generating one is only performing a better adequacy to the dependence between margins structures than the original one.
One remarks that an high censoring level does not really impact the issue of the copula estimation. Indeed, for all DGPs with 70% censoring level, there is only one case where the result changes: DGP 24. In fact, the Gumbel copula is leftskewed, as Student copula may be; which can explains the wrong estimation of the copula in this case.
Example: acupuncture for chronic headache in primary care data
We used the acupuncture for chronic headache in primary care data from Vickers et al. [35–37] containing migraine and chronic tension headache of 401 patients aged 18 to 65 years old who reported an average of at least two headaches per month. Subjects were recruited in the general practice context in England and Wales and they were allocated to receive until 12 acupuncture treatments for a period of three months. For the sake of the study, acupuncture intervention was provided in the community by the United Kingdom National Health Service (NHS). The study starts in 2002 with a time horizon of 12 months and was registered ISRCTN96537534.
Information gained in the analysis process for costs and QALY in both arms
Modelisation process  Control arm  Acupuncture arm 

Kendall’s tau \((\hat {\tau }_{K})\)  −0.1065  −0.1232 
QALY distribution  \(T_{adj} \sim \mathcal {N}(\hat {\mu }_{T_{adj}}=0.7083, \hat {\sigma }_{T_{adj}}=0.1118)\)  \(T_{adj} \sim \mathcal {N}(\hat {\mu }_{T_{adj}}=0.7268, \hat {\sigma }_{T_{adj}}=0.1190)\) 
Costs statistics  \(\hat {\mu }_{C}=217.20\)  \(\hat {\mu }_{C}=403.40 \) 
\(\hat {\sigma }_{C}=486.00\)  \(\hat {\sigma }_{C}=356.59\)  
Costs distribution  \(C \sim logN(\hat {\nu }_{C}=4.4844, \hat {\tau }_{C}=1.3390)\)  \(C \sim logN(\hat {\nu }_{C}=5.7111, \hat {\tau }_{C}=0.7600)\) 
Selected copula family  Gaussian  Student (t) 
Copula parameter \((\hat {\theta })\)  −0.1664  −0.1923 
Impact of an artificially created censoring on the acupuncture example
Information gained in the analysis process QALY in both arms, where an artificial censoring around 30 percents has been created
Modelisation process  Control arm  Acupuncture arm 

Kendall’s tau \((\hat {\tau }_{K})\)  −0.1388  −0.1011 
QALY distribution  \(T_{adj} \sim \mathcal {N}(\hat {\mu }_{T_{adj}}=0.7133, \hat {\sigma }_{T_{adj}}=0.1026)\)  \(T_{adj} \sim \mathcal {N}(\hat {\mu }_{T_{adj}}=0.7304, \hat {\sigma }_{T_{adj}}=0.1005)\) 
Selected copula family  Gaussian  Gaussian 
Copula parameter \((\hat {\theta })\)  −0.2163  −0.1582 
Hence, with a certain level of censoring, the estimation of ICER loses accuracy, and the confidence interval gets larger. However, the estimated ICER in case of censoring stay relatively close to the estimated ICER with original data, and stays in his confidence interval.
Synthesis on the acupuncture example
Firstly, let take a look at the conclusions of Wonderling et al. [36], who were the first to work on these original data. Without totally detailing their methodology, they used a linear regression adjusted on covariates of interest (age, sex, diagnosis, severity of headache at baseline, number of years of headache disorder, baseline SF36 results and geographical site) to evaluate the mean difference for costs and effectiveness (in terms of QALY). Hence, from a linear model, they have got an \(\hat {ICER}\) equals to 9180 £with a mean health gain for acupuncture treatment of 0.021 QALY.
Using the copulabased methodology presented in this paper, we get, with these original data, a mean health gain for acupuncture treatment of 0.026 QALY and when we apply an exponentially distributed censoring around 30% on QALY variable, we get a mean health gain for acupuncture treatment of 0.021 QALY. The major difference between both approaches is on the estimated ICER value. However, we remark that the value of 9180 £keeps in the confidence interval of the copulabased ICER.
Conclusion
One motivation for this work was generated by the limitations of the at standard regression models applied to the costeffectiveness analysis where the dependence structure between costs and utility along with time was not taken into account. We provided a simple stepbystep procedure to find INB and ICER and their confidence intervals, even in case of censoring.
The methodology presented here can easily be implemented on computational software. Indeed, on R, using the packages copula [38] and CDVine [39], one can easily apply the whole process to a dataset, either censored or not. On SAS, the use of a PROC COPULA will be enough to fit a certain copula on a whole dataset.
Appendix
Proposition 1
The random variables \(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) and η _{ ji } are dependent.
Proof
where f is a function of η(ω). Thus, T _{ adj }(ω) is dependent of η(ω). □
Proposition 2
The random variables \(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) and \(\phantom {\dot {i}\!}\eta _{{adj}_{ji}}\) are independent.
Proof
where H is an invertible Borel function (hence monotone). Since that T(ω) and η(ω) are independent, consequently H[ T(ω)] and H[ η(ω)] are independent. Thus, \(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) and \(\phantom {\dot {i}\!}\eta _{{adj}_{ji}}\) are independent. □
Proposition 3
The random variables C _{ ji } and η _{ ji } are dependent.
Proof
Thus, assuming that T(ω) is independent from η(ω), one sees that η(ω) and C(ω) are dependent. □
Abbreviations
 CEA:

Costeffectiveness analysis
 DGP:

Data generating process
 DIC:

Deviance information criterion
 ICER:

Incremental costeffectiveness ratio
 INB:

Incremental net benefit
 IPCW:

Inverse probability censoring weighted
 QALY:

Quality adjusted life years
Declarations
Acknowledgments
The authors wish to thank the National Health Society (UK), through the courtesy of Andrew J. Vickers, for allowing us to use the data on acupuncture for chronic headache in primary care data in this paper.
Funding
This research received no specific grant from any funding agency in the public, commercial, or notforprofit sectors.
Availability of data and materials
Since the acupuncture data are not the property of the authors, they are not available through this article. However, they are available via the supplemental material of Vickers [37] on the link: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1489946/bin/17456215715S1.xls.
Authors’ contributions
CF performed all simulations, statistical analysis, interpretations and drafted the article. Both CF and JPD designed the new methodology presented and the theoretical extensions. PL critically revised the article. All authors have read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Data are deidentified and authors got the authorization from Pr. A.J. Vickers for a secondary analysis of the data.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.
Authors’ Affiliations
References
 Willan A, Lin D. Incremental net benefit in randomized clinical trial. Stat Med. 2001; 20:1563–74.View ArticlePubMedGoogle Scholar
 Kaplan E, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958; 53:457–81.View ArticleGoogle Scholar
 Zhao H, Bang H, Wang H, Pfeifer PE. On the equivalence of some medical cost estimators with censored data. Stat med. 2007; 26(24):4520–30.View ArticlePubMedGoogle Scholar
 Weinstein MC, Stason WB. Foundations of costeffectiveness analysis for health and medical practices. New England J Med. 1977; 296(13):716–21.View ArticleGoogle Scholar
 Willan A, Chen E, Cook R, Lin D. Incremental net benefit in clinical trials with qualityadjusted survival. Stat Med. 2003; 22:353–62.View ArticlePubMedGoogle Scholar
 Willan A, Lin D, Manca A. Regression methods for costeffectiveness analysis with censored data. Stat Med. 2005; 24:131–45.View ArticlePubMedGoogle Scholar
 Pullenayegum E, Willan A. Semiparametric regression models for costeffectiveness analysis: Improving the efficiency of estimation from censored data. Stat Med. 2007; 26:3274–99.View ArticlePubMedGoogle Scholar
 Sklar M. Fonctions de répartition à n dimensions et leurs marges. Publ Inst Statist Univ Paris. 1959; 8:229–31.Google Scholar
 Embrechts P, McNeil A, Straumann D. Correlation and dependence in risk management: properties and pitfalls. Risk management: value at risk and beyond. 2002;:176–223.
 Frees EW, Valdez EA. Understanding relationships using copulas. North Am Actuar J. 1998; 2(1):1–25.View ArticleGoogle Scholar
 Hougaard P. Analysis of multivariate survival data: Springer Science & Business Media; 2012.
 Nelsen RB. An introduction to copulas. Springer Series in Statistics. New York: Springer; 2006.Google Scholar
 Thompson S, Nixon R. How sensitive are costeffectiveness analyses to choice of parametric distributions?Med Decision Making. 2007; 4:416–23.Google Scholar
 Stamey J, Beavers D, Faries D, Price K, Seaman J. Bayesian modeling of costeffectiveness studies with unmeasured confounding: a simulation study. Pharm Stat. 2014; 13:94–100.View ArticlePubMedGoogle Scholar
 Lin D. Linear regression analysis of censored medical costs. Biostatistics. 2000; 1:35–47.View ArticlePubMedGoogle Scholar
 Buckley J, James I. Linear regression with censored data. Biometrika. 1979; 66:429–36.View ArticleGoogle Scholar
 Miller R, Halpern J. Regression with censored data. Biometrika. 1982; 693:521–31.View ArticleGoogle Scholar
 Genest C, Rivest LP. Statistical inference procedures for bivariate archimedean copulas. J Am Stat Assoc. 1993; 88:1034–43.View ArticleGoogle Scholar
 Genest C, Ghoudi K, Rivest LP. A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika. 1995; 82:543–52.View ArticleGoogle Scholar
 Brahimi B, Necir A. A semiparametric estimation of copula models based on the method of moments. Stat Method. 2012; 9(4):467–77.View ArticleGoogle Scholar
 Kendall M. A new measure of rank correlation. Biometrika. 1938; 30:81–9.View ArticleGoogle Scholar
 Oakes D. A concordance test for independence in the presence of censoring. Biometrics. 1982; 38:451–5.View ArticlePubMedGoogle Scholar
 Oakes D. On consistency of kendall’s tau under censoring. Biometrika. 2008; 954:997–1001.View ArticleGoogle Scholar
 LakhalChaieb L. Copula inference under censoring. Biometrika. 2010; 972:505–12.View ArticleGoogle Scholar
 Dos Santos Silva R, Freitas Lopes H. Copula, marginal distributions and model selection: a bayesian note. Stat Comput. 2008; 18:313–20.View ArticleGoogle Scholar
 Genest C, Neslehova J, Ben Ghorbal N. Estimators based on kendall’s tau in multivariate copula models. Aust NZ J Stat. 2011; 53:157–77.View ArticleGoogle Scholar
 El Maache H, Lepage Y. Spearman’s rho and kendall’s tau for multivariate data sets. Math Stat Appl, Lect NotesMonograph Series. 2003; 42:113–30.Google Scholar
 Fieller EC. Some problems in interval estimation. J Stat Royal Soc Ser B. 1954; 16:175–85.Google Scholar
 Willan A, O’Brien B. Confidence intervals for costeffectiveness ratios: an application of fieller’s theorem. Health Econ. 1996; 5:297–305.View ArticlePubMedGoogle Scholar
 Chaudhary N, Stearns S. Confidence intervals for costeffectiveness ratios: an example from randomized trial. Stat Med. 1996; 15:1447–58.View ArticlePubMedGoogle Scholar
 Siani C, Moatti JP. Quelles méthodes de calcul des régions de confiance du ratio coûtefficacité incrémental choisir?: Universites d’AixMarseille II et III; 2002.
 Bang H, Zhao H. Medianbased incremental costeffectiveness ratio (icer). J Stat Theory Prac. 2012; 6(3):428–42.View ArticleGoogle Scholar
 Nixon R, Thompson S. Methods for incorporating covariate adjustment, subgroup analysis and betweencentre differences into costeffectiveness evaluations. Health Econ. 2005; 14:1217–29.View ArticlePubMedGoogle Scholar
 Tsai KT, Peace KE, et al. Analysis of subgroup data in clinical trials. Journal of Causal Inference. 2013; 1(2):193–207. De Gruyter.View ArticleGoogle Scholar
 Vickers AJ, Rees RW, Zollman CE, McCarney R, Smith CM, Ellis N, et al.Acupuncture for chronic headache in primary care: large, pragmatic, randomised trial. Bmj. 2004; 328(7442):744.View ArticlePubMedPubMed CentralGoogle Scholar
 Wonderling D, Vickers A, Grieve R, McCarney R. Cost effectiveness analysis of a randomised trial of acupuncture for chronic headache in primary care. British Med J. 2004; 328:747–9.View ArticleGoogle Scholar
 Vickers AJ. Whose data set is it anyway? sharing raw data from randomized trials. Trials. 2006; 7(1):15.View ArticlePubMedPubMed CentralGoogle Scholar
 Kojadinovic I, Yan J, et al.Modeling multivariate distributions with continuous margins using the copula r package. J Stat Softw. 2010; 34(9):1–20.View ArticleGoogle Scholar
 Brechmann E, Schepsmeier U. Modeling dependence with cand dvine copulas: The rpackage cdvine. J Stat Softw. 2013; 52(3):1–27.View ArticleGoogle Scholar