Copula function
To introduce the theory beyond the new proposed estimator for costeffectiveness 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 ddimensional multivariate distribution \(H(x_{1},x_{2},\ldots,x_{d})=\mathbb {P}(X_{1}\leq x_{1}, X_{2}\leq \)
x
_{2},…,X
_{
d
}≤x
_{
d
}) from a random vector (X
_{1},X
_{2},…,X
_{
d
})^{T} with marginal distributions \(F_{i}(x)=\mathbb {P}(X_{i} \leq x)\) can be written as:
$$ H(x_{1},x_{2},\ldots,x_{d})=C(F_{1}(x_{1}),F_{2}(x_{2}),\ldots,F_{d}(x_{d})) $$
where C is the cumulative distribution functionof 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},\ldots,u_{d})=\mathbb {P}(F_{1}(X_{1})\leq u_{1},F_{2}(X_{2})\leq u_{2},\ldots,F_{d}(X_{d})\leq u_{d})\). Therefore, the copula density can be written as follows:
$$ c(u_{1},u_{2},\ldots,u_{d})=\frac{\partial^{d} C(u_{1},u_{2},\ldots,u_{d})}{\partial u_{1},u_{2},\ldots,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}(\dot)\) 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 practicioners have the following measures, \(\phantom {\dot {i}\!}E_{{adj}_{ji}}\), which is such that:
$$E_{{adj}_{ji}}={inf}_{i} \left[T_{{adj}_{ji}},\eta_{{adj}_{ji}} \right] $$
where \(\phantom {\dot {i}\!}\eta _{{adj}_{ji}}\) represents the censoring adjusted on quality of life and \(\phantom {\dot {i}\!}T_{{adj}_{ji}}\) 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.
\(\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
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 nonadministrative censoring exists as the main consideration of the estimation. Let \(\boldsymbol {Z_{ji}^{C}}\) be the dvector of covariates that affect costs for arm j, j=0,1, for the grouped population, and \(\boldsymbol {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}+\ldots+\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}}+\ldots+\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 α
_{
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}^{\star}}{\hat{G}\left(X_{jki}^{\star}\right)}Z_{j}^{C} \left(Z^{C}_{j}\right)^{t} \right)^{1} \sum_{i=1}^{n} \frac{\delta_{jki}^{\star} C_{jki}}{\hat{G}\left(X_{jki}^{\star}\right)} Z_{j}^{C} $$
where \(X_{jki}^{\star }=\min (X_{ji}, a_{k+1}), \hat {G}(\bullet)\) is the KaplanMeier estimator of and, as described in “Background” section, X
_{
ji
} is the minimum between time from randomization to death and time from randomization to censoring, and . The same approach is used to find β, 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}^{2}_{C_{j}}=\frac{1}{\sum_{l=1}^{n}\delta_{l} 2}\sum_{i=1}^{n} \delta_{i}\left(\hat{e}^{0}_{i}\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}\boldsymbol {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 μ
_{
C
} be the mean and \(\sigma ^{2}_{C}\) 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 \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.
In the case of QALY, the choice may be any symmetric or skewsymmetric 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}_{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.
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 \(\left (C^{\{1\}},T_{{adj}}^{\{1\}}\right)\) and \(\phantom {\dot {i}\!}\left (C^{\{2\}},T_{{adj}}^{\{2\}}\right)\) two independent joint observations of the couple (C,T
_{
adj
}). Then, the pair is said concordant if \(\left (C^{\{1\}}C^{\{2\}}\right)\left (T_{{adj}}^{\{1\}}T_{{adj}}^{\{2\}}\right)>0\) and discordant elsewhere. The Kendall’s tau, which is in fact a concordance measure, is defined in Kendall [21] by
where \(\mathbb {E}\) is the expectation, is the indicator function, and . In a more general frame, one has the couples \(\left (C^{\{1\}},T_{adj}^{\{1\}}\right),\left (C^{\{2\}}, T_{adj}^{\{2\}}\right), \ldots,\left (C^{\{n\}}, T_{adj}^{\{n\}}\right)\) where all the values of \(C^{\{r\}}, T_{adj}^{\{r\}}, r=1,\ldots,n\) are unique. Thus, one can write and where r and s are the index of the independent replications. In absence of censoring, the estimation of τ is given by its sample version:
$$\hat{\tau}_{K}=\binom {n}{2}^{1}\sum_{1\leq r < s \leq 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 to that equation such that:
$$\tilde{\tau}_{K}=\binom {n}{2}^{1}\sum_{1\leq r < s \leq 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 highdependent 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\leq r < s \leq n\}} L_{rs}a_{rs}b_{rs}}{\sum_{\{1\leq r < s \leq 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 goodnessoffit, LakhalChaieb [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
}(χ) the cumulative distribution functions for QALY and costs for a given randomization arm. Thus, one gets:
$${} \begin{aligned} f(y,\chi\Phi)&= c(F_{T_{adj}}(y\Phi_{T_{adj}}),F_{C}(\chi\Phi_{C})\Phi_{\theta})\\ &\quad\times f_{T_{adj}}(y\Phi_{T_{adj}})f_{C}(\chi\Phi_{C}) \end{aligned} $$
where Φ
_{
C
} stands for a parameter vector comprising parameters of the distribution of costs, \(\Phi _{T_{adj}}\) is the parameter vector for QALY distribution, Φ
_{
θ
} 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 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:
$${{\begin{aligned} &\!\! L\left(\boldsymbol{x}\Phi,\mathcal{M}_{k}\right)=\prod_{j=1}^{n} \left[ c\left(F_{T_{adj}}\left(y_{j}\Phi_{T_{adj}},\mathcal{M}_{k}\right),F_{C}\right.\right.\\ &\times\left.\left.\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\!\!\Phi_{\theta},\mathcal{M}_{k}{\vphantom{F_{T_{adj}}(y_{j}\Phi_{T_{adj}},\mathcal{M}_{k}),F_{C}}}\right) f_{T_{adj}}\!\left(y_{j}\Phi_{T_{adj}},\mathcal{M}_{k}\!\right)\!f_{C}\!\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\right]^{\!\delta_{j}} \\ & \times \left[1\,\,C\left(F_{T_{adj}}\left(y_{j}\Phi_{T_{adj}},\mathcal{M}_{k}\right),F_{C}\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\!\!\Phi_{\theta},\mathcal{M}_{k}\right)\!\right]^{1\delta_{j}} \\ & \times \left[ c\left(F_{\eta_{adj}}\left(y_{j}\Phi_{\eta_{adj}},\mathcal{M}_{k}\right),F_{C}\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\Phi_{\theta},\mathcal{M}_{k}\right)\right.\\ &\times\left. f_{\eta_{adj}}\left(y_{j}\Phi_{\eta_{adj}},\mathcal{M}_{k}\right)f_{C}\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\right]^{1\delta_{j}} \\ & \times \left[1C\left(F_{\eta_{adj}}\left(y_{j}\Phi_{\eta_{adj}},\mathcal{M}_{k}\right),F_{C}\left(\chi_{j}\Phi_{C},\mathcal{M}_{k}\right)\Phi_{\theta},\mathcal{M}_{k}\right)\right]^{\delta_{j}}\!, \end{aligned}}} $$
where δ
_{
j
} indicates if an individual is censored or not. Then, using the deviance function which is \(D(\Phi _{k})=2ll(\boldsymbol {x}\Phi,\mathcal {M}_{k})\) where ll stands for the loglikelihood function, Dos Santos Silva et al. [25] proposed to use the deviance information criterion (DIC) which is:
$$\begin{array}{@{}rcl@{}} DIC(\mathcal{M}_{k})=2\mathbb{E}[\!D(\Phi_{k})\boldsymbol{x},\mathcal{M}_{k}]D(\mathbb{E}[\!\Phi_{k}\boldsymbol{x},\mathcal{M}_{k}]). \end{array} $$
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
From the estimated joint densities f(y
_{
j
},χ
_{
j
}),j∈{0,1}, one writes, for costs:
$$\begin{array}{@{}rcl@{}} \mathbb{E}[\!C_{j}] &=& \int_{\mathbb{D}_{C_{j}}} \int_{\mathbb{D}_{T_{{adj}_{j}}}} \chi_{j} f\left(\chi_{j},y_{j}\right) {dy}_{j} d\chi_{j} \\ &\approx & \int_{\mathbb{D}_{C_{j}}} \int_{\mathbb{D}_{T_{{adj}_{j}}}} \chi_{j} c_{\hat{\theta}}^{(i)}\left(\tilde{F}_{T_{adj}}\left(y\hat{\Phi}_{T_{adj}}\right),\tilde{F}_{C}\left(\chi\hat{\Phi}_{C}\right)\right)\\ &&\times\tilde{f}_{C}\left(y\hat{\Phi}_{C}\right) \tilde{f}_{T_{adj}}\left(\chi\hat{\Phi}_{T_{adj}}\right) {dy}_{j} d\chi_{j} \end{array} $$
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
} for arm j. For QALY, one has the following:
$$\begin{array}{@{}rcl@{}} \mathbb{E}[\!E_{j}]&=&\int_{\mathbb{D}_{T_{{adj}_{j}}}} \int_{\mathbb{D}_{C_{j}}} y_{j} f\left(\chi_{j},y_{j}\right) d\chi_{j} {dy}_{j} \\ &\approx & \int_{\mathbb{D}_{T_{{adj}_{j}}}} \int_{\mathbb{D}_{C_{j}}} y_{j} c_{\hat{\theta}}^{(i)}\left(\tilde{F}_{T_{adj}}\left(y\hat{\Phi}_{T_{adj}}\right),\tilde{F}_{C}\left(\chi\hat{\Phi}_{C}\right)\right)\\ &&\times\tilde{f}_{C}\left(y\hat{\Phi}_{c}\right) \tilde{f}_{T_{adj}}\left(\chi\hat{\Phi}_{T_{adj}}\right) d\chi_{j} {dy}_{j}. \end{array} $$
Thus, given the expected costs and survival time adjusted to quality of life, the adjusted ICER is estimated by
$$\widehat{ICER}=\frac{\widehat{\mathbb{E}(C_{j=1})}\widehat{\mathbb{E}(C_{j=0})}} {\widehat{\mathbb{E}(T_{{adj}_{j=1}})}\widehat{\mathbb{E}(E_{{adj}_{j=0}})}} $$
and using Fieller’s theorem (Fieller [28], Willan and O’Brien [29], Chaudhary and Stearns [30]), one gets the 100(1−α)% confidence interval such that
$$\widehat{ICER}\left(\frac{\left(1z^{2}_{1\alpha /2} \hat{\sigma}_{\Delta_{C}\Delta_{T_{adj}}} \right) \pm z_{1\alpha /2} \sqrt{\hat{\sigma}^{2}_{\Delta_{T_{adj}}}+\hat{\sigma}^{2}_{\Delta_{C}}2\hat{\sigma}^{2}_{\Delta_{C} \Delta_{T_{adj}}} z^{2}_{1\alpha /2}\big(\hat{\sigma}^{2}_{\Delta_{T_{adj}}}\hat{\sigma}^{2}_{\Delta_{C}}\hat{\sigma}^{2}_{\Delta_{C} \Delta_{T_{adj}}}\big) }}{1z^{2}_{1\alpha /2}\hat{\sigma}^{2}_{\Delta_{T_{adj}}}}\right). $$
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.
Let \(\phantom {\dot {i}\!}T_{{adj}_{j=1,k=1,i}}\) be the effectiveness for smokers individuals i on the treated arm, \(\phantom {\dot {i}\!}T_{{adj}_{j=1,k=0,i}}\) be the effectiveness for nonsmokers individuals i on the treated arm; and the same for the control arm. One does a similar writing for costs. Let \(\mathbb {E}(\overline {T}_{{adj}_{j}})=\mathbb {E}(T_{{adj}_{j,k=1}})\mathbb {E}(T_{{adj}_{j,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 \(\overline {INB}(\lambda)=\lambda (\mathbb {E}(\overline {T}_{{adj}_{j=1}})\mathbb {E}(\overline {T}_{{adj}_{j=0}}))(\mathbb {E}(\overline {C}_{j=1})\mathbb {E}(\overline {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{array}{@{}rcl@{}} \mathbb{V}ar(\overline{INB}(\lambda)) &=& \lambda^{2}\sigma^{2}_{\Delta_{\overline{T_{adj}}}}+\sigma^{2}_{\Delta_{\overline{C}}}2\lambda\sigma_{\Delta_{\overline{T_{adj}}} \Delta_{\overline{C}}} \\ &=& \lambda^{2} \mathbb{V}ar\left(\overline{T}_{{adj}_{j=1}}\overline{T}_{{adj}_{j=0}}\right)\\ &&+\mathbb{V}ar\left(\overline{C}_{j=1}\overline{C}_{j=0}\right) \\ && 2\lambda cov\left(\overline{T}_{{adj}_{j=1}}\overline{T}_{{adj}_{j=0}},\overline{C}_{j=1}\overline{C}_{j=0}\right) \end{array} $$
where variances are computed in the standard path. For the covariance term, \(\sigma _{\Delta _{\overline {T_{adj}}} \Delta _{\overline {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} &cov\left(\overline{T}_{{adj}_{j=1}}\overline{T}_{{adj}_{j=0}},\overline{C}_{j=1}\overline{C}_{j=0}\right)\\ &\qquad\quad= cov\left(\overline{T}_{{adj}_{j=1}}, \overline{C}_{j=1}\right) + cov\left(\overline{T}_{{adj}_{j=0}}, \overline{C}_{j=0}\right) \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{array}{@{}rcl@{}} cov(\overline{T}_{{adj}_{j=1}}&&\overline{T}_{{adj}_{j=0}},\overline{C}_{j=1}\overline{C}_{j=0})\\ =&& cov\left(\overline{T}_{{adj}_{j=1}}, \overline{C}_{j=1}\right)+ cov\left(\overline{T}_{{adj}_{j=0}}, \overline{C}_{j=0}\right)\\ && cov\left(\overline{T}_{{adj}_{j=1}}, \overline{C}_{j=0}\right)cov\left(\overline{T}_{{adj}_{j=0}}, \overline{C}_{j=1}\right). \end{array} $$
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.
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 α if
$$\frac{\left\overline{INB}(\lambda)\right}{\sqrt{\mathbb{V}ar(\overline{INB}(\lambda))}} $$
is greater than the z
_{1−α/2} percentile of a standard gaussian distribution.