# On the censored cost-effectiveness analysis using copula information

Charles Fontaine, Jean-Pierre Daurès, Paul Landais

► To cite this version:

Charles Fontaine, Jean-Pierre Daurès, Paul Landais. On the censored cost-effectiveness analysis using copula information. BMC Medical Research Methodology, 2017, 17 (1), pp.27. 10.1186/s12874-017-0305-9. hal-01896706

HAL Id: hal-01896706

https://hal.science/hal-01896706

Submitted on 27 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

# On the censored cost-effectiveness analysis using copula information

Charles Fontaine* (ID), Jean-Pierre Daurès and Paul Landais

## 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 cost-effectiveness analysis (CEA), which has shown its utility in many medico-economic studies and where censoring often occurs.

**Methods:** This paper discusses a copula-based modeling of the joint density and an estimation method of the costs, and quality adjusted life years (QALY) in a cost-effectiveness 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 cost-effectiveness 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 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:** Cost-effectiveness analysis, Censored data, Copulas, Parametric models, Subgroups analysis

## Background

Due to the variety of treatments for a specific health problem and in conjunction with their increasing costs, cost-effectiveness 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 common-in-use therapy.

In the last decades, the cost-effectiveness analysis (CEA) of new treatments became an actual subject of work for

statisticians. It is used in two particular designs: decision modeling-based CEA and trial-based CEA. The major difference between both approaches is that in the case of trial-based 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 modeling-based CEA are based on easily generalizable data.

The cost-effectiveness analysis is used to measure the incremental cost-effectiveness ratio (ICER) and the incremental net benefit (INB). The ICER is defined as the ratio:

$$ICER = \frac{E(C_1) - E(C_0)}{E(T_1) - E(T_0)}$$

<footnote>*Correspondence: charles fontaine@usherbrooke.ca
UPRES EA2415-Institut Universitaire de Recherche Clinique, Université de Montpellier, 641, Av. du doyen G.-Giraud, Montpellier, France</footnote>
where C₁ is the cost of the tested therapeutic, C₀ is the cost of the control group which is usually measured in term of a given monetary unit, T₁ is the effectiveness of the tested therapeutic which is usually measured in term of survival life years, T₀ is the effectiveness of the control group and, E(•) is the expectation function. Therefore, it is an indicator of the monetary cost of using a new therapy in terms of survival time. On the other hand, the INB is defined as the following difference:

$$INB = \lambda(T_1 - T_0) - (C_1 - C_0)$$

where λ is the willingness-to-pay 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 product-limit estimator [2] and then to estimate the survival time by integrating the survival functions until a time τ (i.e. μⱼ is the life expectancy until τ), the maximal time of follow-up. Concerning the estimation of costs in case of censoring, many estimators may be used. Zhao et al. [3] have shown the equivalence among them.

The quality adjusted life years (QALY) concept was introduced in 1977 by Weinstein and Stason [4], and is still actually one of the most important notions in the cost-effectiveness theoretical framework. In the paper of Willan et al. [5], the concept of quality of life adjusted to the survival time is defined as follows. Let qⱼᵢ be the quality adjusted survival for the period of interest for the patient i which follows the treatment j and let qⱼⱼᵢ be the observed quality of life for the patient i receiving treatment j, j = 0, 1 during the interval of time aₖ. In fact, this is a representation of the standard survival times scaled down by the quality of life experienced by patients. One note that the duration of interest of a study (0, τ] is divided in K (arbitrary, relative to the data) sub-intervals [aₖ, aₖ₊₁) where k = 1, 2, ..., K and where 0 = a₁ < a₂ < ... < aₖ₊₁ = τ. Thus, one can determine the value of qⱼⱼᵢ in the following path: let a patient i be on treatment j with Bⱼᵢ quality of life measured at times tⱼᵢ₁, tⱼᵢ₂, ..., tⱼᵢⱼᵢ with respective scores Qⱼᵢ₁, Qⱼᵢ₂, ..., Qⱼᵢⱼᵢ. These scores are nothing more than the utility values. Therefore, qⱼⱼᵢ = ∫ₐₖ aₖ₊₁ Q(t)dt is a weighted sum of times spent in the different quality of life states where:

$$Q(t) = \begin{cases} Q_{ji1} & \text{if } 0 \le t < t_{ji1}; \\ Q_{jih} + \frac{(Q_{ji,h+1} - Q_{jih})(t_{ji,h+1} - t_{jih})}{t_{ji,h+1} - t_{jih}} & \text{if } t_{jih} \le t < t_{ji,h+1}; \\ Q_{jim_ji} & \text{if } t_{jim_ji} \le t < X_{ji}; \\ 0 & \text{if } \ge X_{ji}, \end{cases}$$

and where Xⱼᵢ = min(Tⱼᵢ, ηⱼᵢ), δⱼᵢ = 1[Tⱼᵢ < ηⱼᵢ] where η represents the censoring random variable. Furthermore, let Yⱼⱼᵢ = 1(Xⱼᵢ ≥ aₖ and [Xⱼᵢ ≥ aₖ₊₁ or δⱼᵢ = 1]), which indicates if a patient i on the treatment arm j is alive at time aₖ and is not censored on [aₖ, aₖ₊₁) and let Yⱼⱼ = Σₖᵢ=₁⁻ᵏ Yⱼᵢ. If one notes ȳⱼⱼ = Σₖᵢ=₁⁻ᵏ (Yⱼᵢqⱼᵢ)/Yⱼⱼ, then using a known expression of variance [5], one obtains the estimation of μⱼ, the expected value of effectiveness adjusted to QALY, with

$$\hat{\mu}_j = \sum_{k=1}^{K} \hat{S}_j(a_k) \bar{q}_{jk}.$$

More recently, Willan et al. [6] proposed to realize the whole cost-effectiveness analysis using linear regression methods. Let Cᵢ be the observed cost for patient i, then E(Cⱼᵢ) = βₐᵢTᵢ Zₐⱼᵢ, i = 1, 2, ..., nⱼ, Zₐⱼᵢ is a vector of covariates affecting costs and βₐⱼ 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 βₐⱼ which is the mean difference in costs between randomization groups adjusted to other covariates, Δ̂ₐ, 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, Δ̂ₐ. Thus, they proposed to estimate the ICER adjusted for the quality of life by Δ̂ₐ/Δ̂ₐ 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 ȳₐ = λΔ̂ₐ - Δ̂ₐ with variance ȳ̂ₐ² = λ²ȳ̂ₐ² + ȳ̂ₐ² - 2λȳΔ̂ₐΔ̂ₐ. Thus, if the INB is positive, the therapy is cost-effective and if it is negative, there is no cost-effectiveness. One remarks that it is possible to determine the cost-effectiveness under a significance level α using a statistic of test: if ȳₐ/ȳₐ is greater than the test level z₁-α, at a level α, the therapy is cost-effective compared to the standard.

From that linear regression approach, many variants exist either parametric or semi-parametric [7]. The main problem of the estimator for cost-effectiveness analysis based on the linear regression is that even if the inverse-probability censored weighting estimator is consistent, it is not efficient because if an individual is censored before or at time aₖ₊₁, 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 non-linearity. In the case of a semi-parametric 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 cost-effectiveness 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

To introduce the theory beyond the new proposed estimator for cost-effectiveness analysis, it is crucial to present the concept of modeling the dependence among two or more variables, namely the copula function concept. The idea of the copula started with Sklar [8] who formulated his famous theorem as follows: a d-dimensional multivariate distribution $H(x_1, x_2, \dots, x_d) = \mathbb{P}(X_1 \le x_1, X_2 \le x_2, \dots, X_d \le x_d)$ from a random vector $(X_1, X_2, \dots, X_d)^T$ with marginal distributions $F_i(x) = \mathbb{P}(X_i \le x)$ can be written as:

$$H(x_1, x_2, \dots, x_d) = C(F_1(x_1), F_2(x_2), \dots, F_d(x_d))$$

where C is the cumulative distribution function of the copula. In fact, it is a cumulative distribution function from $[0, 1]^d$ to $[0, 1]$ with uniform margins over the unit interval $[0, 1]$ such that: $C(u_1, u_2, \dots, u_d) = \mathbb{P}(F_1(X_1) \le u_1, F_2(X_2) \le u_2, \dots, F_d(X_d) \le u_d)$. Therefore, the copula density can be written as follows:

$$c(u_1, u_2, \dots, u_d) = \frac{\partial^d C(u_1, u_2, \dots, u_d)}{\partial u_1, \partial u_2, \dots, \partial u_d}.$$

Thus, for a bivariate model, the joint density function of stochastic variables $X_1$ and $X_2$ is:

$$f(x_1, x_2) = c(F_1(x_1), F_2(x_2))f_1(x_1)f_2(x_2).$$

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

If the survival time adjusted for the quality of life is already measured, one should directly estimate the parameters of its distribution. However, most of the time, practitioners only have two variables: time and quality of life. As shown in the previous section, the classical adjustment method for time on quality of life is given by $T_{adj}(w) = \int_0^{T(w)} Q(v(t))dt$ where $Q(v(t))$ is the adjustment of quality of life scores on the interval of time of interest. Since that the function $H(t) = \int_0^t Q(v(y))dy$ is monotonically increasing, it is possible to rewrite the cumulative distribution function of $T_{adj}$ as a composition of functions such as:

$$F_{adj}(y) = F \circ H^{-1}(y)$$

where $H^{-1}(\cdot)$ is the generalized inverse function, and the probability density function such that:

$$f_{adj}(y) = f[H^{-1}(y)] \frac{1}{Q[v(H^{-1}(y))]}$$

where $f_{adj}$ is the density function of $T_{adj}$ and $f$ is the density function of $T$. Therefore, for an individual $i$ on treatment $j$, the practitioners have the following measures, $E_{adjji}$, which is such that:

$$E_{adjji} = inf_i [T_{adjji}, \eta_{adjji}]$$

where $\eta_{adjji}$ represents the censoring adjusted on quality of life and $T_{adjji}$ the survival time having the same adjustment. Also, let $C_{ji}$ the cumulative cost for individual $i$ on arm $j$. Thus, we get the following dependency evidences:

1. $T_{adjji}$ and $\eta_{ji}$ are dependent,

2. $T_{adjji}$ and $\eta_{adjji}$ are independent,

3. $C_{ji}$ and $\eta_{ji}$ are dependent.

Proofs are provided in the “Appendix” section.
## Estimation of the parameters of the distributions

Even if, to begin, the right distributions for costs and QALY are unknown, it is possible to infer their two main parameters: mean and variance. In fact, we will consider here each arm of the trial as a distinct random variable with distinct mean and variance but with the same probability distribution. Furthermore, we will assume that non-administrative censoring exists as the main consideration of the estimation. Let $Z_{ji}^C$ be the d-vector of covariates that affect costs for arm j, j = 0, 1, for the grouped population, and $Z_{ji}^E$ be the one for QALY. Then, as proposed by Thompson and Nixon [13] and Stamey et al. [14], the costs mean function, on an arm j, is defined as:

$$\mu_j^C = \alpha_0 + \alpha_1 z_{1j}^C + \dots + \alpha_d z_{dj}^C,$$

and, the QALY variable mean function given costs is defined as:

$$\mu_j^{T_{adj}} = \beta_0 + \beta_1 z_{1j}^{T_{adj}} + \dots + \beta_d z_{dj}^{T_{adj}}.$$

As these models are, in fact, linear regression models with censoring on covariates, using the method of Lin [15], one can estimate the regression coefficients vector $\alpha_C$ by the sum over the k periods of time of interest $\hat{\alpha}_C = \sum_{k=1}^K \hat{\alpha}_{C_k}$ using an inverse probability weighting method (IPCW) such that, for an individual i belonging to arm j,

$$\hat{\alpha}_{C_k} = \left( \sum_{i=1}^{n} \frac{\delta_{jki}^{*}}{\hat{G}(X_{jki}^{*})} Z_{j}^{C} (Z_{j}^{C})^{t} \right)^{-1} \sum_{i=1}^{n} \frac{\delta_{jki}^{*} C_{jki}}{\hat{G}(X_{jki}^{*})} Z_{j}^{C}$$

where $X_{jki}^* = \min(X_{ji}, a_{k+1})$, $\hat{G}(\bullet)$ is the Kaplan-Meier estimator of $G(\bullet)$, $\delta_{jki}^* = \delta_{ji} + (1 - \delta_{ji}) \mathbb{1}(X_{ji} \ge a_{k+1})$ and, as described in “Background” section, $X_{ji}$ is the minimum between time from randomization to death and time from randomization to censoring, and $\delta_{ji} = \mathbb{1}(T_{ji} \le \eta_{ji})$. The same approach is used to find $\beta$, the vector of regression coefficients for QALY. Thus, from the inference on coefficients, it is possible to determine the adjusted mean on survival.

In terms of variance, we propose the use of the result of Buckley and James [16] (see also Miller and Halpern [17]), which is a generalization of the IPCW techniques.

Thus, the approximate variance of the cost distribution on a given arm is:

$$\hat{\sigma}_{C_j}^2 = \frac{1}{\sum_{l=1}^{n} \delta_l - 2} \sum_{i=1}^{n} \delta_i \left( \hat{e}_i^0 - \frac{1}{\sum_{l=1}^{n} \delta_l} \sum_{j=1}^{n} \delta_j \hat{e}_j^0 \right)^2$$

where $\hat{e}_0^i$ is an error term such that $\hat{e}_0^i = C_i - Z_j^C \hat{\alpha}_j$. A similar approach is done for QALY.

## Determination of the parametric distributions

To model costs, three common distributions are frequently used: Gamma, Normal and Lognormal distributions. Their parametrization is easily done given the mean and the variance of the distribution. Let $\mu_C$ be the mean and $\sigma_C^2$ be the variance of costs for any clinical arm j. Then, the parametric distribution choice will be one of the following:

1. $C_j \sim Normal(\mu_{C_j}, \sigma_{C_j}^2)$,

2. $C_j \sim Gamma(\mu_{C_j}, \rho_{C_j})$,

3. $C_j \sim Lognormal(v_{C_j}, \tau_{C_j}^2)$,

where $v_C$ and $\tau_C^2$ are mean and variance of the log-costs, i.e. $v_C = 2\log(\mu_C) - \frac{1}{2}\log(\sigma_C^2 + \mu_C^2)$ and $\tau_C^2 = \log(\sigma_C^2 + \mu_C^2) - 2\log(\mu_C)$. Furthermore, $\rho_C$ is the shape parameter of the Gamma distribution, which is such that $\rho_C = \mu_C^2/\sigma_C^2$. 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 log-likelihood.

In the case of QALY, the choice may be any symmetric or skew-symmetric distribution. We propose two options here, but the classical way here is to consider only a gaussian distribution following the work of Thompson and Nixon [13]. One notes that, even if it may look strange to use a real defined function for a real positive distribution, the distribution of $T_{adj}$ usually has a high mean with a low standard error such that the negative part of the fitted distribution is in fact negligible. Then, the proposed options are:

1. $T_{adj} \sim Normal(\mu_{T_{adj}}, \sigma_{T_{adj}}^2)$,

2. $T_{adj} \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 non-randomized 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.

The Kendall's tau inversion method to infer a copula parameter has been shown to be consistent under specific assumptions, which holds here [20]. In fact, for every parametric copula family, there exists a direct relationship between the copula parameter and the Kendall's tau. For example, with Clayton copula, one has $\hat{\theta} = 2\tau/(1 - \tau)$ and for Gumbel copula, $\hat{\theta} = 1/(1 - \tau)$. These relationships are clearly given in almost all the literature about parametric copulas [12]. Let consider the couple of random variables $(C, T_{adj})$ on a fixed therapeutic arm $j, j = 0, 1$. Furthermore, let consider $(C^{(1)}, T_{adj}^{(1)})$ and $(C^{(2)}, T_{adj}^{(2)})$ two independent joint observations of the couple $(C, T_{adj})$. Then, the pair is said concordant if $(C^{(1)} - C^{(2)})(T_{adj}^{(1)} - T_{adj}^{(2)}) > 0$ and discordant elsewhere. The Kendall's tau, which is in fact a concordance measure, is defined in Kendall [21] by

$$ 
\begin{aligned}
\tau_K &= \mathbb{P}\left[ (C^{(1)} - C^{(2)})(T_{adj}^{(1)} - T_{adj}^{(2)}) > 0 \right] \\
&\quad - \mathbb{P}\left[ (C^{(1)} - C^{(2)})(T_{adj}^{(1)} - T_{adj}^{(2)}) < 0 \right] \\
&= 2 \cdot \mathbb{P}\left[ (C^{(1)} - C^{(2)})(T_{adj}^{(1)} - T_{adj}^{(2)}) > 0 \right] - 1 \\
&= \mathbb{E}\left[ \left( 2 \cdot 1[C^{(1)} - C^{(2)} > 0] - 1 \right) \right. \\
&\qquad \left. \times \left( 2 \cdot 1[T_{adj}^{(1)} - T_{adj}^{(2)} > 0] - 1 \right) \right] \\
&= \mathbb{E}[a_{12}b_{12}]
\end{aligned}
$$

where $\mathbb{E}$ is the expectation, $1$ is the indicator function, $a_{12} = 2 \cdot 1[C^{(1)} - C^{(2)} > 0] - 1$ and $b_{12} = 2 \cdot 1[T_{adj}^{(1)} - T_{adj}^{(2)} > 0] - 1$. In a more general frame, one has the couples $(C^{(1)}, T_{adj}^{(1)})$, $(C^{(2)}, T_{adj}^{(2)})$, ..., $(C^{(n)}, T_{adj}^{(n)})$ where all the values of $C^{(r)}, T_{adj}^{(r)}, r = 1, ..., n$ are unique. Thus, one can write $a_{rs} = 2 \cdot 1[C^{(r)} - C^{(s)} > 0] - 1$ and $b_{rs} = 2 \cdot 1[T_{adj}^{(r)} - T_{adj}^{(s)} > 0] - 1$ where $r$ and $s$ are

the index of the independent replications. In absence of censoring, the estimation of $\tau$ is given by its sample version:

$$ \hat{\tau}_K = \binom{n}{2}^{-1} \sum_{1 \le r < s \le n} a_{rs} b_{rs}, $$

where $n$ is the sample size. In fact, it is simply the $n(n - 1)/2$ pairs of bivariate observations that could be constructed, multiplied by the subtraction of the number of discordant pairs to the number of concordant pairs. Under censoring, the approach of Oakes [22] propose to add an uncensoring indicator $L_{rs} = 1$ [$\min(C^{(r)}, C^{(s)}) < \min(U_C^{(r)}, U_C^{(s)}), \min(T_{adj}^{(r)}, T_{adj}^{(s)}) < \min(U_{T_{adj}^{(r)}}, U_{T_{adj}^{(s)}})$] to that equation such that:

$$ \hat{\tau}_K = \binom{n}{2}^{-1} \sum_{1 \le r < s \le n} L_{rs} a_{rs} b_{rs}, $$

where $U^{(r)}$ and $U^{(s)}$ represent the censoring variables under each independent replication. The problem with that estimator is the lack of consistency on a high-dependent scheme. Therefore, one advocates the use of the renormalized Oakes' estimator [23] for which consistency has been shown. Therefore, the estimator

$$ \hat{\tau}_K = \frac{\sum_{1 \le r < s \le n} L_{rs} a_{rs} b_{rs}}{\sum_{1 \le r < s \le n} L_{rs}} $$

is simply the ratio of the subtraction of the number of uncensored discordant pairs to the number of uncensored concordant pairs over the total of uncensored pairs.

## Bayesian selection of the copula

Inference here may be done on any consistent parametric copula. To select the right copula, few alternatives are proposed in the literature. For complete data distributions, Genest and Rivest [18] proposed a procedure based on a probability integral transformation. Furthermore, based on the goodness-of-fit, Lakhal-Chaieb [24] proposed a procedure for censored data when the distributions are estimated by the survival functions. However, when available, a prior knowledge of the distribution of the margins in a copula is not a negligible information. It should be taken into account for the inference of the copula model when the latter is unknown, to minimize the risk of errors. In their paper, Dos Santos Silva et al. [25] proposed a
method of selection based on the information criterion. Let note $F_{T_{adj}}(y)$ and $F_C(\chi)$ the cumulative distribution functions for QALY and costs for a given randomization arm. Thus, one gets:

$$f(y, \chi | \Phi) = c(F_{T_{adj}}(y | \Phi_{T_{adj}}), F_C(\chi | \Phi_C) | \Phi_{\theta}) \\ \times f_{T_{adj}}(y | \Phi_{T_{adj}}) f_C(\chi | \Phi_C)$$

where $\Phi_C$ stands for a parameter vector comprising parameters of the distribution of costs, $\Phi_{T_{adj}}$ is the parameter vector for QALY distribution, $\Phi_{\theta}$ is the dependence parameter which is functional of the Kendall's tau and $\Phi = \Phi_C \cup \Phi_{T_{adj}} \cup \Phi_{\theta}$ is the union of all these vectors of parameters. Also, $f$ and $F$ respectively stand for the density probability function and the cumulative probability function. Note that a similar writing could be done for a multivariate modeling. As shown in Genest et al. [26], one can find the copula parameter of any parametric copula while just having the Kendall's tau [21] measure, even in multivariate models [27].

Let $\mathbf{x}$ be a bivariate sample of size $n$ of that density function. Also, let $\mathcal{M}_k$ a copula model for $k = 1..m$ where $m$ is the quantity of models one wants to test. Therefore, the likelihood function is given by:

$$L(\mathbf{x}|\Phi, \mathcal{M}_k) = \prod_{j=1}^{n} \left[ c(F_{T_{adj}}(y_j|\Phi_{T_{adj}}, \mathcal{M}_k), F_C \times (\chi_j|\Phi_C, \mathcal{M}_k)|\Phi_{\theta}, \mathcal{M}_k) f_{T_{adj}}(y_j|\Phi_{T_{adj}}, \mathcal{M}_k) f_C(\chi_j|\Phi_C, \mathcal{M}_k) \right]^{\delta_j} \times \left[1-C(F_{T_{adj}}(y_j|\Phi_{T_{adj}}, \mathcal{M}_k), F_C(\chi_j|\Phi_C, \mathcal{M}_k)|\Phi_{\theta}, \mathcal{M}_k)\right]^{1-\delta_j} \times \left[c(F_{\eta_{adj}}(y_j|\Phi_{\eta_{adj}}, \mathcal{M}_k), F_C(\chi_j|\Phi_C, \mathcal{M}_k)|\Phi_{\theta}, \mathcal{M}_k) \times f_{\eta_{adj}}(y_j|\Phi_{\eta_{adj}}, \mathcal{M}_k) f_C(\chi_j|\Phi_C, \mathcal{M}_k)\right]^{1-\delta_j} \times \left[1-C(F_{\eta_{adj}}(y_j|\Phi_{\eta_{adj}}, \mathcal{M}_k), F_C(\chi_j|\Phi_C, \mathcal{M}_k)|\Phi_{\theta}, \mathcal{M}_k)\right]^{\delta_j},$$

where $\delta_j$ indicates if an individual is censored or not. Then, using the deviance function which is $D(\Phi_k) = -2ll(\mathbf{x}|\Phi, \mathcal{M}_k)$ where $ll$ stands for the log-likelihood function, Dos Santos Silva et al. [25] proposed to use the deviance information criterion (DIC) which is:

$$DIC(\mathcal{M}_k) = 2\mathbb{E}[D(\Phi_k)|\mathbf{x}, \mathcal{M}_k] - D(\mathbb{E}[\Phi_k|\mathbf{x}, \mathcal{M}_k]).$$

They proposed to use the Monte Carlo approximations to $\mathbb{E}[D(\Phi_k)|\mathbf{x}, \mathcal{M}_k]$ and $\mathbb{E}[\Phi_k|\mathbf{x}, \mathcal{M}_k]$ which are respectively $L^{-1} \sum_{l=1}^{L} D(\Phi_k^l)$ and $L^{-1} \sum_{l=1}^{L} \Phi_k^l$. Here, one supposes that $\{\Phi_k^{(1)}, \dots, \Phi_k^{(L)}\}$ is a sample from the posterior distribution $f(\Phi_k|\mathbf{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 dependence parametrisation in a cost-effectiveness 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 cost-effectiveness ratio

From the estimated joint densities $f(y_j, \chi_j), j \in \{0, 1\}$, one writes, for costs:

$$\mathbb{E}[C_j] = \int_{\mathbb{D}_{C_j}} \int_{\mathbb{D}_{T_{adj_j}}} \chi_j f(\chi_j, y_j) dy_j d\chi_j \\ \approx \int_{\mathbb{D}_{C_j}} \int_{\mathbb{D}_{T_{adj_j}}} \chi_j c_{\theta}^{(i)}(\tilde{F}_{T_{adj}}(y|\hat{\Phi}_{T_{adj}}), \tilde{F}_C(\chi|\hat{\Phi}_C)) \\ \times \tilde{f}_C(y|\hat{\Phi}_C) \tilde{f}_{T_{adj}}(\chi|\hat{\Phi}_{T_{adj}}) dy_j d\chi_j$$

where $\mathbb{D}_{C_j}$ and $\mathbb{D}_{T_{adj_j}}$ are respectively the domain of definition for the random variables $C$ and $T_{adj_j}$ for arm $j$. For QALY, one has the following:

Thus, given the expected costs and survival time adjusted to quality of life, the adjusted ICER is estimated by

$$\widehat{\text{ICER}} = \frac{\mathbb{E}(C_{j=1}) - \mathbb{E}(C_{j=0})}{\mathbb{E}(T_{adj_j=1}) - \mathbb{E}(T_{adj_j=0})}$$

and using Fieller's theorem (Fieller [28], Willan and O'Brien [29], Chaudhary and Stearns [30]), one gets the $100(1 - \alpha)\%$ confidence interval such that
$$ \widehat{ICER} \left( \frac{\left( 1 - z_{1-\alpha/2}^2 \hat{\sigma}_{\Delta_C \Delta_T_{adj}} \right) \pm z_{1-\alpha/2} \sqrt{\hat{\sigma}_{\Delta_T_{adj}}^2 + \hat{\sigma}_{\Delta_C}^2 - 2\hat{\sigma}_{\Delta_C \Delta_T_{adj}}^2 - z_{1-\alpha/2}^2 (\hat{\sigma}_{\Delta_T_{adj}}^2 \hat{\sigma}_{\Delta_C}^2 - \hat{\sigma}_{\Delta_C \Delta_T_{adj}}^2)}}{1 - z_{1-\alpha/2}^2 \hat{\sigma}_{\Delta_T_{adj}}^2} \right). $$

In this formula, $z_{1-\alpha/2}$ represents the 100(1 - $\alpha$/2) percentile of the standard normal distribution. Furthermore, $\hat{\sigma}_{\Delta_T_{adj}}^2$ represents the variance of the distribution for effectiveness where $\Delta_{T_{adj}}$ is the difference between $T_{adj=1}$ and $T_{adj=0}$. The same scheme arises for $\hat{\sigma}_{\Delta_C}^2$. 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 non-parametric ways), even in problematic cases. However, for non-common 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 non-sufficient statistics in a parametric confidence interval model.

## Incremental net benefit

The adjusted INB($\lambda$) is estimated by $\widehat{\text{INB}} = \lambda(\mathbb{E}(\widehat{T_{adj=1}}) - \mathbb{E}(\widehat{T_{adj=0}})) - (\mathbb{E}(\widehat{C_{j=1}}) - \mathbb{E}(\widehat{C_{j=0}}))$ with variance $\hat{\sigma}_{\lambda}^2 = \lambda^2 \hat{\sigma}_{\Delta_T_{adj}}^2 + \hat{\sigma}_{\Delta_C}^2 - 2\lambda \hat{\sigma}_{\Delta_C \Delta_T_{adj}}$ where $\lambda$ is the willingness-to-pay 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 cost-effectiveness analysis while achieving a discrimination between two or more subgroups. The principle is that there exists a baseline variable $Z_k$, $k \in \{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,

non-smokers in the treated group, and the same for the control groups.

Let $T_{adj=1,k=1,i}$ be the effectiveness for smokers individuals $i$ on the treated arm, $T_{adj=1,k=0,i}$ be the effectiveness for non-smokers individuals $i$ on the treated arm; and the same for the control arm. One does a similar writing for costs. Let $\mathbb{E}(\widehat{T_{adj}}) = \mathbb{E}(T_{adj_{k=1}}) - \mathbb{E}(T_{adj_{k=0}})$ and the same for costs. Then, the interest in this discrimination is on the incremental net benefit marginalized to the smokers cohort, which is $\widehat{\text{INB}}(\lambda) = \lambda(\mathbb{E}(\widehat{T_{adj=1}}) - \mathbb{E}(\widehat{T_{adj=0}})) - (\mathbb{E}(\widehat{C_{j=1}}) - \mathbb{E}(\widehat{C_{j=0}}))$. Since that subgroups are inside clinical arms, the main issue is to establish the expression of variance. Adjusting Fieller's method to the subgroups context, one has

$$ \begin{aligned} \mathbb{V}ar(\widehat{\text{INB}}(\lambda)) &= \lambda^2 \sigma_{\Delta_T_{adj}}^2 + \sigma_{\Delta_C}^2 - 2\lambda \sigma_{\Delta_T_{adj}} \Delta_C \\ &= \lambda^2 \mathbb{V}ar(\widehat{T_{adj=1}} - \widehat{T_{adj=0}}) \\ &\quad + \mathbb{V}ar(\widehat{C_{j=1}} - \widehat{C_{j=0}}) \\ &\quad - 2\lambda \mathbb{cov}(\widehat{T_{adj=1}} - \widehat{T_{adj=0}}, \widehat{C_{j=1}} - \widehat{C_{j=0}}) \end{aligned} $$

where variances are computed in the standard path. For the covariance term, $\sigma_{\Delta_T_{adj} \Delta_C}$, there are two possible scenarios. Firstly, when the assumption that subgroups in treated and control arms are randomized is possible, one has

$$ \begin{aligned} \operatorname{cov}(\widehat{T}_{adj=1} - \widehat{T}_{adj=0}, \widehat{C}_{j=1} - \widehat{C}_{j=0}) \\ = \operatorname{cov}(\widehat{T}_{adj=1}, \widehat{C}_{j=1}) + \operatorname{cov}(\widehat{T}_{adj=0}, \widehat{C}_{j=0}) \end{aligned} $$

which can be found easily using standard techniques. Secondly, when the randomization assumption between subgroup is not possible because the cohorts are unbalanced in the clinical arms, then the covariance is:

$$ \begin{aligned} \operatorname{cov}(\widehat{T}_{adj=1} - \widehat{T}_{adj=0}, \widehat{C}_{j=1} - \widehat{C}_{j=0}) \\ = \operatorname{cov}(\widehat{T}_{adj=1}, \widehat{C}_{j=1}) + \operatorname{cov}(\widehat{T}_{adj=0}, \widehat{C}_{j=0}) \\ - \operatorname{cov}(\widehat{T}_{adj=1}, \widehat{C}_{j=0}) - \operatorname{cov}(\widehat{T}_{adj=0}, \widehat{C}_{j=1}). \end{aligned} $$

Here, the terms $\operatorname{cov}(\widehat{T}_{adj=1}, \widehat{C}_{j=1})$ and $\operatorname{cov}(\widehat{T}_{adj=0}, \widehat{C}_{j=0})$ are computed in a standard way similarly to the randomized case. However, for the crossed-arms covariance terms, $\operatorname{cov}(\widehat{T}_{adj=1}, \widehat{C}_{j=0})$ and $\operatorname{cov}(\widehat{T}_{adj=0}, \widehat{C}_{j=1})$, the approach that we suggest is to estimate the cumulative
joint distributions $F(\bar{E}_{j=1}, \bar{C}_{j=0})$ by $C_{\hat{\theta}}(\hat{F}_{\hat{E}}(y_{j=1}), \hat{F}_{\hat{C}}(\chi_{j=0}))$ and $F(\bar{T}_{adj=1}, \bar{C}_{j=0})$ by $C_{\hat{\theta}}(\hat{F}_{\bar{T}_{adj}}(y_{j=1}), \hat{F}_{\bar{C}}(\chi_{j=0}))$ and $F(\bar{T}_{adj=1}, \bar{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(\bar{C}, \bar{T}_{adj}) = E[\bar{T}_{adj}, \bar{C}] - E[\bar{T}_{adj}] E[\bar{C}]$, compute the desired covariance from the estimated joint density function and the estimated marginal density functions.

In the spirit of the test of Willan et al. [6], one can test the equality of the INB between cohorts, which can be rejected at a level $\alpha$ if

$$ \frac{|\overline{INB}(\lambda)|}{\sqrt{\overline{Var}(\overline{INB}(\lambda))}} $$

is greater than the $z_{1-\alpha/2}$ percentile of a standard gaussian distribution.

## 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_{\hat{\theta}}^{(*)}(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 $(\star)$ 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 $\theta$ 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.

In order to evaluate the performance of the proposed method in non-trivial cases, we performed Monte-Carlo simulations on 27 different simulation schemes. The methodology was to simulate bivariate data (representing costs and QALY) from three specific copulas. For each copula, we simulated the three possible configurations for the marginal distributions (costs are either normally distributed or lognormally distributed or follow a gamma distribution, while QALY is normally distributed). Then, we applied three different levels of randomly censoring on the marginal distribution of QALY (15, 30 and 70%). Hence, there were nine possible copula configurations; for these nine, we also challenged our methodology without any censoring to get a point of comparison for the inference of the Kendall's tau part of the model. For all the simulations, we assumed that $T_{adj}$ follows a normal distribution. Alternatively, we could have set $T_{adj}$ following a Gamma distribution. However, in the present core, the goal relatively to margins was to check that the determination of the marginal distributions criterion was adequate,

which may be tested on only one margin (to simplify the text). The censoring followed an exponential distribution such that $\lambda_{s=15} = 0.041$, $\lambda_{s=30} = 0.090$ and $\lambda_{s=70} = 0.308$ where $s$ represents the censoring percentage simulated. For all data generating processes (DGP), the Kendall's tau was identical and represented an intermediate level of dependence between marginal distributions to be fair with the reality: $\tau_K = 0.60$. Then, we computed the right copula parameter, for each copula, based on this Kendall's tau. We also used a relatively standard mean and variance in cost-effectiveness analysis, following parameters used by Thompson and Nixon [13], such that $\mu_C = 1500$, $\sigma_C = 400$; $\mu_{T_{adj}} = 4$, $\sigma_{T_{adj}} = 0.75$, and we parametrized each marginal distribution to keep close to these values. For the choice of generating copulas, we selected the three most-known ones with the biggest difference: Gaussian, Clayton and Gumbel copulas. The DGPs schemes is shown on Table 1.

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 Kaplan-Meier 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

The proposed way to infer Kendall's tau under censoring gives results close to the real $\tau_K$ measured on data just before applying the censoring. On Fig. 1 and Table 2, one can see the dispersion of the computed Kendall's tau. We remind the reader that the theoretical value of tau used for simulations is 0.6, and that we performed simulations without censoring to give an idea of the Kendall's tau measure on complete data. Thus, firstly, all simulations were performed without any censoring. Secondly, a 15% rate of censoring were applied on data. Thirdly, for all the simulations, we used a 30% rate of censoring and, for the last case, we applied a 70% rate of censoring. Therefore, the vector that contains all the inferred Kendall's tau for a (un)censoring level has a length of 4500. One observes that the estimated values of $\hat{\tau}_K$ are relatively close to the theoretical value of 0.6, but one remarks that the dispersion of values for some $\hat{\theta}$ is not close to the expected value
<table><caption>**Table 1** Scheme of the 27 data generated cases</caption><thead><tr><th>Generating copula</th><th>Costs distribution</th><th>Censoring level</th><th>DGP</th></tr></thead><tbody><tr><td rowspan="6">Gaussian copula<br>θ ≈ 0.809</td><td>F<sub>C</sub> ~ N(μ<sub>C</sub> = 1500, σ<sub>C</sub> = 400)</td><td>15%</td><td>DGP 1</td></tr><tr><td></td><td>30%</td><td>DGP 2</td></tr><tr><td></td><td>70%</td><td>DGP 3</td></tr><tr><td>F<sub>C</sub> ~ Γ shape<sub>C</sub> = 12, scale<sub>C</sub> = 125)</td><td>15%</td><td>DGP 4</td></tr><tr><td></td><td>30%</td><td>DGP 5</td></tr><tr><td></td><td>70%</td><td>DGP 6</td></tr><tr><td rowspan="6">Clayton copula<br>θ = 3</td><td>F<sub>C</sub> ~ logN(v<sub>C</sub> = 7.30, τ<sub>C</sub> = 0.25)</td><td>15%</td><td>DGP 7</td></tr><tr><td></td><td>30%</td><td>DGP 8</td></tr><tr><td></td><td>70%</td><td>DGP 9</td></tr><tr><td>F<sub>C</sub> ~ N(μ<sub>C</sub> = 1500, σ<sub>C</sub> = 400)</td><td>15%</td><td>DGP 10</td></tr><tr><td></td><td>30%</td><td>DGP 11</td></tr><tr><td></td><td>70%</td><td>DGP 12</td></tr><tr><td rowspan="6">Gumbel copula<br>θ ≈ 0.809</td><td>F<sub>C</sub> ~ Γ shape<sub>C</sub> = 12, scale<sub>C</sub> = 125)</td><td>15%</td><td>DGP 13</td></tr><tr><td></td><td>30%</td><td>DGP 14</td></tr><tr><td></td><td>70%</td><td>DGP 15</td></tr><tr><td>F<sub>C</sub> ~ logN(v<sub>C</sub> = 7.30, τ<sub>C</sub> = 0.25)</td><td>15%</td><td>DGP 16</td></tr><tr><td></td><td>30%</td><td>DGP 17</td></tr><tr><td></td><td>70%</td><td>DGP 18</td></tr><tr><td rowspan="6">Gumbel copula<br>θ ≈ 0.809</td><td>F<sub>C</sub> ~ N(μ<sub>C</sub> = 1500, σ<sub>C</sub> = 400)</td><td>15%</td><td>DGP 19</td></tr><tr><td></td><td>30%</td><td>DGP 20</td></tr><tr><td></td><td>70%</td><td>DGP 21</td></tr><tr><td>F<sub>C</sub> ~ Γ shape<sub>C</sub> = 12, scale<sub>C</sub> = 125)</td><td>15%</td><td>DGP 22</td></tr><tr><td></td><td>30%</td><td>DGP 23</td></tr><tr><td></td><td>70%</td><td>DGP 24</td></tr><tr><td rowspan="6">Gumbel copula<br>θ ≈ 0.809</td><td>F<sub>C</sub> ~ logN(v<sub>C</sub> = 7.30, τ<sub>C</sub> = 0.25)</td><td>15%</td><td>DGP 25</td></tr><tr><td></td><td>30%</td><td>DGP 26</td></tr><tr><td></td><td>70%</td><td>DGP 27</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Inference on Kendall's tau according to the simulations. Dispersion of the estimated Kendall's tau according to the censoring level around the theoretical value of tau used to generate data: 0.6</figcaption></figure>
<table><caption>**Table 2** Information on the estimation of the Kendall's tau for each censoring level</caption><thead><tr><th>Censoring level</th><th>Mean (τ̂<sub>K</sub>)</th><th>Var (τ̂<sub>K</sub>)</th><th>Min (τ̂<sub>K</sub>)</th><th>Max (τ̂<sub>K</sub>)</th></tr></thead><tbody><tr><td>Censoring =0%</td><td>0.6002</td><td>0.00019</td><td>0.5488</td><td>0.6476</td></tr><tr><td>Censoring =15%</td><td>0.6011</td><td>0.00024</td><td>0.5416</td><td>0.6539</td></tr><tr><td>Censoring =30%</td><td>0.6030</td><td>0.00035</td><td>0.5319</td><td>0.6648</td></tr><tr><td>Censoring =70%</td><td>0.6089</td><td>0.00146</td><td>0.4624</td><td>0.7257</td></tr></tbody></table>

The outputs of the 9 simulations with a censoring of 0% are joint together in the information on the first line and the same for simulations censored at 15% on the second line, 30% at the third line and 70% at the last line

as is the dispersion on $\hat{\tau}_K$, since $\tau_K$ ranges in $[-1, 1]$ while, for example, for Clayton copula, $\theta$ ranges in $[-1, \infty) \setminus \{0\}$.

### Inference on the marginal distribution for the costs

To select the right marginal distribution for costs on each of the 500 simulations and for each DGP, we used the proposed criteria based on the deviance. Thus, on Fig. 2, one can see the performance of this criterion. One may remark that, even with a 70% rate of censoring, the chosen parametric distribution is almost always correctly estimated.

### Inference on the copula models

There exists a bunch of parametric copula families, but for the purpose of these simulations, we limited our selection to the most well-known ones: Gaussian, Student, Clayton, Gumbel, Frank and Joe copulas. However, with a real dataset, the reader may test any consistent

parametric copula using the proposed way. During simulations, for each iteration, we collected the information about the selected family using the information criterion. At Tables 3 and 4, one can see these results for each DGP. We indicated, in bold, the copula that was selected the most in the 500 iterations and claimed that it was the copula to be retained for the mentioned DGP. When the selection of the copula was only in-between those used to simulate the data (as on Table 3), it was obvious that the chosen copula family is the one used for the generation process. Otherwise, when intermediate copulas (for which dependence adequation is close to the one provoked by the generating copula) are introduced in the selection process, the results may differ as shown on Table 4. For the DGPs where the costs were simulated following a normal distribution (1,2,3,10,11,12,19,20 and 21), the choice of the copula was influenced by the dependence caused by the type of copula retained to generate the data. In facts, when the costs margin follows $N(\mu_C = 1500, \sigma_C = 400)$ and the QALY margin follows $N(\mu_{T_{adj}} = 4, \sigma_{T_{adj}} = 0.75)$, using $\tau_K = 0.60$, if the dependency distribution is mostly normal (i.e. comes from a Gaussian or a Gumbel copula), there will be a moderate tail on both tails of the distribution with an uniform cloud along the correlation path. For these reasons which are attributes of the Frank copula, it was the chosen copula for DGPs 1,2,3,19,20 and 21. However, when, from the generating Clayton copula, a strict dependence was imposed to the left-tail while a mostly total independence was imposed to the right-tail, only a Clayton copula appeared appropriate to modeling

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Frequency of selection of parametric marginal distributions for costs from the deviance criteria for each data generating process (DGP). The black bar represents the selection of the Normal distribution, the dark-gray bar represents the selection of the Gamma distribution and the light-gray stands for a logNormal distribution</figcaption></figure>
<table><caption>**Table 3** Frequency of choice of a copula for each DGP given 500 iterations for the three main copulas</caption><thead><tr><th>DGP</th><th>Gaussian copula</th><th>Clayton copula</th><th>Gumbel copula</th></tr></thead><tbody><tr><td>DPG 1</td><td><strong>460</strong></td><td>0</td><td>40</td></tr><tr><td>DPG 2</td><td><strong>454</strong></td><td>0</td><td>46</td></tr><tr><td>DPG 3</td><td>25</td><td>94</td><td><strong>381</strong></td></tr><tr><td>DPG 4</td><td><strong>415</strong></td><td>0</td><td>85</td></tr><tr><td>DPG 5</td><td><strong>429</strong></td><td>0</td><td>71</td></tr><tr><td>DPG 6</td><td>140</td><td>21</td><td><strong>339</strong></td></tr><tr><td>DPG 7</td><td><strong>473</strong></td><td>0</td><td>27</td></tr><tr><td>DPG 8</td><td><strong>490</strong></td><td>0</td><td>10</td></tr><tr><td>DPG 9</td><td><strong>261</strong></td><td>10</td><td>229</td></tr><tr><td>DPG 10</td><td>0</td><td><strong>500</strong></td><td>0</td></tr><tr><td>DPG 11</td><td>0</td><td><strong>500</strong></td><td>0</td></tr><tr><td>DPG 12</td><td>0</td><td><strong>499</strong></td><td>1</td></tr><tr><td>DPG 13</td><td>18</td><td><strong>482</strong></td><td>0</td></tr><tr><td>DPG 14</td><td>11</td><td><strong>489</strong></td><td>0</td></tr><tr><td>DPG 15</td><td>68</td><td><strong>371</strong></td><td>61</td></tr><tr><td>DPG 16</td><td>2</td><td><strong>498</strong></td><td>0</td></tr><tr><td>DPG 17</td><td>6</td><td><strong>493</strong></td><td>1</td></tr><tr><td>DPG 18</td><td>177</td><td><strong>531</strong></td><td>62</td></tr><tr><td>DPG 19</td><td>2</td><td>0</td><td><strong>498</strong></td></tr><tr><td>DPG 20</td><td>9</td><td>0</td><td><strong>491</strong></td></tr><tr><td>DPG 21</td><td>1</td><td>9</td><td><strong>490</strong></td></tr><tr><td>DPG 22</td><td>123</td><td>0</td><td><strong>377</strong></td></tr><tr><td>DPG 23</td><td>109</td><td>0</td><td><strong>391</strong></td></tr><tr><td>DPG 24</td><td>36</td><td>0</td><td><strong>467</strong></td></tr><tr><td>DPG 25</td><td>80</td><td>0</td><td><strong>420</strong></td></tr><tr><td>DPG 26</td><td>84</td><td>0</td><td><strong>416</strong></td></tr><tr><td>DPG 27</td><td>121</td><td>0</td><td><strong>379</strong></td></tr></tbody></table>

The chosen copula is in bold font

these data. That's the reason why it was the chosen copula for DGPs 10,11 and 12. In the situation where QALY follows $N(\mu_E = 4, \sigma_E = 0.75)$ while costs follows a skewed marginal distribution (gamma or lognormal), in any case, there was not a right-tailed distribution, but mostly a left-oriented data cloud with a fat left-tail. That is the reason why, even when the data were generated from a Clayton copula, the Student (t) copula was the chosen one. Therefore, the simulations showed that the Bayesian criterion of selection of the copula was in accordance with the theoretical properties of the parametric copulas.

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

<table><caption>**Table 4** Frequency of choice of copula for each DGP on 500 iterations beyond the three main copulas and intermediate copulas</caption><thead><tr><th>DGP</th><th>Gaussian copula</th><th>Student copula</th><th>Clayton copula</th><th>Gumbel copula</th><th>Frank copula</th><th>Joe copula</th></tr></thead><tbody><tr><td>DPG 1</td><td>7</td><td>45</td><td>6</td><td>25</td><td><strong>415</strong></td><td>2</td></tr><tr><td>DPG 2</td><td>2</td><td>44</td><td>4</td><td>25</td><td><strong>425</strong></td><td>1</td></tr><tr><td>DPG 3</td><td>4</td><td>54</td><td>4</td><td>17</td><td><strong>421</strong></td><td>0</td></tr><tr><td>DPG 4</td><td>32</td><td><strong>297</strong></td><td>0</td><td>56</td><td>112</td><td>3</td></tr><tr><td>DPG 5</td><td>26</td><td><strong>291</strong></td><td>0</td><td>59</td><td>120</td><td>4</td></tr><tr><td>DPG 6</td><td>23</td><td><strong>284</strong></td><td>0</td><td>79</td><td>113</td><td>1</td></tr><tr><td>DPG 7</td><td>47</td><td><strong>343</strong></td><td>1</td><td>45</td><td>38</td><td>26</td></tr><tr><td>DPG 8</td><td>55</td><td><strong>344</strong></td><td>0</td><td>43</td><td>33</td><td>25</td></tr><tr><td>DPG 9</td><td>55</td><td><strong>328</strong></td><td>0</td><td>56</td><td>37</td><td>24</td></tr><tr><td>DPG 10</td><td>0</td><td>12</td><td><strong>364</strong></td><td>0</td><td>124</td><td>0</td></tr><tr><td>DPG 11</td><td>0</td><td>17</td><td><strong>372</strong></td><td>0</td><td>111</td><td>0</td></tr><tr><td>DPG 12</td><td>0</td><td>5</td><td><strong>385</strong></td><td>0</td><td>110</td><td>0</td></tr><tr><td>DPG 13</td><td>0</td><td><strong>427</strong></td><td>56</td><td>0</td><td>17</td><td>0</td></tr><tr><td>DPG 14</td><td>2</td><td><strong>433</strong></td><td>55</td><td>0</td><td>10</td><td>0</td></tr><tr><td>DPG 15</td><td>1</td><td><strong>431</strong></td><td>58</td><td>0</td><td>10</td><td>0</td></tr><tr><td>DPG 16</td><td>1</td><td><strong>477</strong></td><td>19</td><td>2</td><td>1</td><td>0</td></tr><tr><td>DPG 17</td><td>2</td><td><strong>476</strong></td><td>18</td><td>0</td><td>4</td><td>0</td></tr><tr><td>DPG 18</td><td>4</td><td><strong>472</strong></td><td>19</td><td>0</td><td>4</td><td>1</td></tr><tr><td>DPG 19</td><td>0</td><td>18</td><td>0</td><td>110</td><td><strong>196</strong></td><td>176</td></tr><tr><td>DPG 20</td><td>0</td><td>20</td><td>0</td><td>107</td><td><strong>214</strong></td><td>159</td></tr><tr><td>DPG 21</td><td>0</td><td>17</td><td>0</td><td>109</td><td><strong>202</strong></td><td>172</td></tr><tr><td>DPG 22</td><td>13</td><td><strong>236</strong></td><td>0</td><td>153</td><td>67</td><td>31</td></tr><tr><td>DPG 23</td><td>19</td><td><strong>231</strong></td><td>0</td><td>157</td><td>58</td><td>35</td></tr><tr><td>DPG 24</td><td>23</td><td>161</td><td>0</td><td><strong>204</strong></td><td>75</td><td>37</td></tr><tr><td>DPG 25</td><td>32</td><td><strong>198</strong></td><td>0</td><td>177</td><td>13</td><td>80</td></tr><tr><td>DPG 26</td><td>39</td><td><strong>213</strong></td><td>0</td><td>148</td><td>13</td><td>87</td></tr><tr><td>DPG 27</td><td>37</td><td><strong>210</strong></td><td>0</td><td>167</td><td>17</td><td>69</td></tr></tbody></table>

The chosen copula is in bold font

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 left-skewed, 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
<table><caption>**Table 5** Information gained in the analysis process for costs and QALY in both arms</caption><thead><tr><th>Modelisation process</th><th>Control arm</th><th>Acupuncture arm</th></tr></thead><tbody><tr><td>Kendall's tau (τ<sub>K</sub>)</td><td>-0.1065</td><td>-0.1232</td></tr><tr><td>QALY distribution</td><td>T<sub>adj</sub> ~ N(μ<sub>T<sub>adj</sub></sub> = 0.7083, σ<sub>T<sub>adj</sub></sub> = 0.1118)</td><td>T<sub>adj</sub> ~ N(μ<sub>T<sub>adj</sub></sub> = 0.7268, σ<sub>T<sub>adj</sub></sub> = 0.1190)</td></tr><tr><td>Costs statistics</td><td>μ<sub>C</sub> = 217.20<br>σ<sub>C</sub> = 486.00</td><td>μ<sub>C</sub> = 403.40<br>σ<sub>C</sub> = 356.59</td></tr><tr><td>Costs distribution</td><td>C ~ logN(ν<sub>C</sub> = 4.4844, τ<sub>C</sub> = 1.3390)</td><td>C ~ logN(ν<sub>C</sub> = 5.7111, τ<sub>C</sub> = 0.7600)</td></tr><tr><td>Selected copula family</td><td>Gaussian</td><td>Student (t)</td></tr><tr><td>Copula parameter (θ̂)</td><td>-0.1664</td><td>-0.1923</td></tr></tbody></table>

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.

The data collection focuses on the measure of effectiveness in terms of QALY gained and the cumulative cost associated in UK pounds (£). Patients themselves reported unit costs associated with non-prescription drugs and private healthcare visits. The cost of the study intervention was estimated from the standard cost for a NHS professional multiplied by the contact time with the patient. Thus, patients in the treated arm had a mean time of 4.2 hours with study acupuncturist. No imputation for missing data was done if the three questionnaires were not complete and consequently for which QALY cannot be measured. Therefore, in the acupuncture arm, there was 136 participants and in the control arm, 119. The modeling process of both joint distributions function for QALY and costs in the two clinical arms is presented on Table 5. We remark that in both cases, dependence is weak since that Kendall's tau ranges between -0.10 and -0.15. Here, for the distribution of $T_{adj}$, we compared two possible choices: a Gamma distribution and a gaussian one. The normal distribution has the smallest deviance. For the choice of distributions for costs, in the two arms, the lognormal distribution was considered as the one with the smallest deviance. For the copula family selection using deviance information criteria, we compared, for each arm, Gaussian, Clayton, Student, Frank, Joe and Gumbel copulas. In the treated arm, Student copula was the considered one while in the control arm, it was the Gaussian copula. Therefore, the joint distribution function of the acupuncture arm was estimated by:

$$ \hat{F}(C_{j=1}, T_{adjj=1}) = C_{\hat{\theta}=-0.1923}^{(\text{Student})} \times (F_C \sim \log N(\hat{v}_C = 5.7111, \hat{\tau}_C = 0.7600), F_{T_{adj}} \sim N(\hat{\mu}_{T_{adj}} = 0.7268, \hat{\sigma}_{T_{adj}} = 0.1190)) $$

while, for the control arm, the estimation is:

$$ \hat{F}(C_{j=0}, T_{adjj=0}) = C_{\hat{\theta}=-0.1664}^{(\text{Gaussian})} \times (F_C \sim \log N(\hat{v}_C = 4.4844, \hat{\tau}_C = 1.3390), F_{T_{adj}} \sim N(\hat{\mu}_{T_{adj}} = 0.7083, \hat{\sigma}_{T_{adj}} = 0.1118)) $$

Using the approach given from the copula densities, the ICER is estimated such that $ICER = 10082.68\$/$ unit of effectiveness, with a confidence interval which is:

$$ ICER = \frac{1 + 12.44z_{1-\alpha}^2 \pm z_{1-\alpha/2} \sqrt{365336.81 - 9744.36 \times z_{1-\alpha/2}^2}}{1 - z_{1-0.0266\alpha/2}^2} $$

where $z_{1-\alpha}$ is the 100(1 - $\alpha$/2) percentile of the standard gaussian distribution. It means that it costs approximately 10082.68 per year to get an additional unit of effectiveness using acupuncture for headache.

The plot of the estimated INB with his 90 percent confidence limits versus lambda is presented on Fig. 3.

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Plot of INB versus $\lambda$ for acupuncture for headaches in primary care example. $\lambda$ stands for the willingness-to-pay for a unit of effectiveness</figcaption></figure>
<table><caption>**Table 6** Information gained in the analysis process QALY in both arms, where an artificial censoring around 30 percents has been created</caption><thead><tr><th>Modelisation process</th><th>Control arm</th><th>Acupuncture arm</th></tr></thead><tbody><tr><td>Kendall's tau (τ<sub>k</sub>)</td><td>-0.1388</td><td>-0.1011</td></tr><tr><td>QALY distribution</td><td>T<sub>adj</sub> ~ N(μ̂<sub>T<sub>adj</sub></sub> = 0.7133, σ̂<sub>T<sub>adj</sub></sub> = 0.1026)</td><td>T<sub>adj</sub> ~ N(μ̂<sub>T<sub>adj</sub></sub> = 0.7304, σ̂<sub>T<sub>adj</sub></sub> = 0.1005)</td></tr><tr><td>Selected copula family</td><td>Gaussian</td><td>Gaussian</td></tr><tr><td>Copula parameter (θ̂)</td><td>-0.2163</td><td>-0.1582</td></tr></tbody></table>

The vertical intercept shows the negative value of the variability of costs and its confidence interval while the horizontal intercept shows the estimated ICER. Since the number of covariates was limited, it has not been possible to determine the existence of subgroups. Whether available, the search for subgroups is possible using our approach.

### Impact of an artificially created censoring on the acupuncture example

To challenge our methodology via a censored dataset, we decided to create an artificial censoring on QALY variable

(since that costs are assumed to always be observed). This censoring variable followed an exponential distribution where $\lambda = 0.30$. Hence, the QALY variable has been censored on approximately 30% of data. Using the same methodology than for the acupuncture original dataset, we obtained the information shown on Table 6. One notes that we did not report the costs information on this table since that it does not change from the information on Table 5 (costs is not a censored variable). Using these information, one obtains an estimated ICER such that $IC̃ER = 10879.89£/ unit of effectiveness, with a confidence interval which is:

<figure><img src="image_6.png" /></figure>
$$ICER = \frac{\left(1 + 11.33z_{1-\alpha}^2 \pm z_{1-\alpha/2} \sqrt{360533.67 - 7443.98 \times z_{1-\alpha/2}^2}\right)}{1 - z_{1-0.0206\alpha/2}^2}$$

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 SF-36 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 ICER equals to 9180 £with a mean health gain for acupuncture treatment of 0.021 QALY.

Using the copula-based 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 copula-based ICER.

## Conclusion

One motivation for this work was generated by the limitations of the at standard regression models applied to the cost-effectiveness analysis where the dependence structure between costs and utility along with time was not taken into account. We provided a simple step-by-step procedure to find INB and ICER and their confidence intervals, even in case of censoring.

On Fig. 4, one sees the schematized method from the observational data produced by both clinical arms to the complete cost-effectiveness analysis. In a parallel way, one accomplishes steps 1 and 2 which are the measure of the dependence between QALY and cumulative costs in each arm via Kendall's tau and the determination of the marginal distributions of both random variables. Then, at step 3, one generates copulas from information gained in steps 1 and 2 and, using the information criterion, one selects the closest copula to the right joint distribution function at step 4. Finally, at step 5, one determines the INB and the ICER using joint cdf. In case of subgroups cohorts analysis, one reiterate the procedure from step 1 to 5 to get two supplementary copulas and then, the joint

cdf of costs and QALY for the crossed-arms covariance terms.

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 $T_{adjji}$ and $\eta_{ji}$ are dependent.

*Proof* To simplify the notation, let assume the following being for an individual $i$ on a therapeutic $j$. Also, let $T_{adj}(\omega) = \int_{0}^{T(\omega)} Q(t)dt$. Given the observed times $E(\omega) = \inf(\eta(\omega), T(\omega))$. If $\eta(\omega) \le T(\omega)$, one gets

$$ \begin{aligned} T_{adj}(\omega) &= \int_{0}^{\eta(\omega)} Q(t)dt + \int_{\eta(\omega)}^{T(\omega)} Q(t)dt \\ &= \eta_{adj}(\omega) + \int_{\eta(\omega)}^{T(\omega)} Q(t)dt \\ &= \eta_{adj}(\omega) + f(\eta(\omega)) \end{aligned} $$

where $f$ is a function of $\eta(\omega)$. Thus, $T_{adj}(\omega)$ is dependent of $\eta(\omega)$. $\square$

**Proposition 2** The random variables $T_{adjji}$ and $\eta_{adjji}$ are independent.

*Proof* Let

$$ \begin{aligned} T_{adj}(\omega) &= \int_{0}^{T(\omega)} Q(t)dt \\ &= H[T(\omega)] \end{aligned} $$

and

$$ \begin{aligned} \eta_{adj}(\omega) &= \int_{0}^{\eta(\omega)} Q(t)dt \\ &= H[\eta(\omega)] \end{aligned} $$

where $H$ is an invertible Borel function (hence monotone). Since that $T(\omega)$ and $\eta(\omega)$ are independent, consequently $H[T(\omega)]$ and $H[\eta(\omega)]$ are independent. Thus, $T_{adjji}$ and $\eta_{adjji}$ are independent. $\square$

**Proposition 3** The random variables $C_{ji}$ and $\eta_{ji}$ are dependent.

*Proof* To simplify the notation, let assume the following be for an individual $i$ on a therapeutic arm $j$. One has:

$$ C(\omega) = \begin{cases} \int_{0}^{T(\omega)} C_k(t)dt & \text{if } T(\omega) \le \eta(\omega); \\ \int_{0}^{\eta(\omega)} C_k(t)dt & \text{if } \eta(\omega) \le T(\omega). \end{cases} $$
Then, with $E(\omega) = \inf(\eta(\omega), T(\omega))$,

$$C_k(E(\omega)) = 1_{[T(\omega) \le \eta(\omega)]} C(T(\omega)) + 1_{[\eta(\omega) \le T(\omega)]} C(\eta(\omega)).$$

Thus, assuming that $T(\omega)$ is independent from $\eta(\omega)$, one sees that $\eta(\omega)$ and $C(\omega)$ are dependent. ☐

### Abbreviations

CEA: Cost-effectiveness analysis; DGP: Data generating process; DIC: Deviance information criterion; ICER: Incremental cost-effectiveness ratio; INB: Incremental net benefit; IPCW: Inverse probability censoring weighted; QALY: Quality adjusted life years

### 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 not-for-profit 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/1745-6215-7-15-S1.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 de-identified and authors got the authorization from Pr. A.J. Vickers for a secondary analysis of the data.

Received: 20 June 2016 Accepted: 2 February 2017

Published online: 15 February 2017

### References

1. Willan A, Lin D. Incremental net benefit in randomized clinical trial. Stat Med. 2001;20:1563-74.

2. Kaplan E, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53:457-81.

3. 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.

4. Weinstein MC, Stason WB. Foundations of cost-effectiveness analysis for health and medical practices. New England J Med. 1977;296(13):716-21.

5. Willan A, Chen E, Cook R, Lin D. Incremental net benefit in clinical trials with quality-adjusted survival. Stat Med. 2003;22:353-62.

6. Willan A, Lin D, Manca A. Regression methods for cost-effectiveness analysis with censored data. Stat Med. 2005;24:131-45.

7. Pullenayegum E, Willan A. Semi-parametric regression models for cost-effectiveness analysis: Improving the efficiency of estimation from censored data. Stat Med. 2007;26:3274-99.

8. Sklar M. Fonctions de répartition à n dimensions et leurs marges. Publ Inst Statist Univ Paris. 1959;8:229-31.

9. 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.

10. Frees EW, Valdez EA. Understanding relationships using copulas. North Am Actuar J. 1998;2(1):1-25.

11. Hougaard P. Analysis of multivariate survival data: Springer Science & Business Media; 2012.

12. Nelsen RB. An introduction to copulas. Springer Series in Statistics. New York: Springer; 2006.

13. Thompson S, Nixon R. How sensitive are cost-effectiveness analyses to choice of parametric distributions? Med Decision Making. 2007;4:416-23.

14. Stamey J, Beavers D, Faries D, Price K, Seaman J. Bayesian modeling of cost-effectiveness studies with unmeasured confounding: a simulation study. Pharm Stat. 2014;13:94-100.

15. Lin D. Linear regression analysis of censored medical data. Biostatistics. 2000;1:35-47.

16. Buckley J, James I. Linear regression with censored data. Biometrika. 1979;66:429-36.

17. Miller R, Halpern J. Regression with censored data. Biometrika. 1982;69-3:521-31.

18. Genest C, Rivest LP. Statistical inference procedures for bivariate archimedean copulas. J Am Stat Assoc. 1993;88:1034-43.

19. Genest C, Ghoudi K, Rivest LP. A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika. 1995;82:543-52.

20. Brahimi B, Necir A. A semiparametric estimation of copula models based on the method of moments. Stat Method. 2012;9(4):467-77.

21. Kendall M. A new measure of rank correlation. Biometrika. 1938;30:81-9.

22. Oakes D. A concordance test for independence in the presence of censoring. Biometrics. 1982;38:451-5.

23. Oakes D. On consistency of kendall's tau under censoring. Biometrika. 2008;95:4-997-1001.

24. Lakhal-Chaieb L. Copula inference under censoring. Biometrika. 2010;97:2:505-12.

25. Dos Santos Silva R, Freitas Lopes H. Copula, marginal distributions and model selection: a bayesian note. Stat Comput. 2008;18:313-20.

26. 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.

27. El Maache H, Lepage Y. Spearman's rho and kendall's tau for multivariate data sets. Math Stat Appl, Lect Notes-Monograph Series. 2003;42:113-30.

28. Fieller EC. Some problems in interval estimation. J Stat Royal Soc Ser B. 1954;16:175-85.

29. Willan A, O'Brien B. Confidence intervals for cost-effectiveness ratios: an application of fieller's theorem. Health Econ. 1996;5:297-305.

30. Chaudhary N, Stearns S. Confidence intervals for cost-effectiveness ratios: an example from randomized trial. Stat Med. 1996;15:1447-58.

31. Siani C, Moatti JP. Quelles méthodes de calcul des régions de confiance du ratio coût-efficacité incrémental choisir?: Universites d'Aix-Marseille II et III; 2002.

32. Bang H, Zhao H. Median-based incremental cost-effectiveness ratio (icer). J Stat Theory Prac. 2012;6(3):428-42.

33. Nixon R, Thompson S. Methods for incorporating covariate adjustment, subgroup analysis and between-centre differences into cost-effectiveness evaluations. Health Econ. 2005;14:1217-29.

34. Tsai K-T, Peace KE, et al. Analysis of subgroup data in clinical trials. Journal of Causal Inference. 2013;1(2):193-207. De Gruyter.

35. 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.

36. 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.

37. Vickers AJ. Whose data set is it anyway? sharing raw data from randomized trials. Trials. 2006;7(1):15.

38. 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.

39. Brechmann E, Schepsmeier U. Modeling dependence with c-and d-vine copulas: The r-package cdvine. J Stat Softw. 2013;52(3):1-27.