Survival Copula model

Note: Reference to be added

This vignette focuses on the copula-based approach by (Stoye) for the meta-analysis of diagnostic test accuracy studies with multiple thresholds. The goal of this approach is to model diagnostic accuracy by directly incorporating the distribution of test values in diseased and non-diseased populations, while flexibly capturing the dependence structure between them.

Unlike standard approaches that model sensitivity and specificity directly, this method is based on an underlying representation of diagnostic test results as interval-censored data. The model combines parametric assumptions for the distributions of test values with a copula-based formulation to describe the joint behavior of the two populations. This allows for a flexible modeling of both marginal distributions and their dependence, while naturally accounting for multiple thresholds within studies.

Firstly, we introduce the model itself, explaining how diagnostic test results are represented as interval-censored observations and how the likelihood is constructed using a copula-based framework. We then describe the choice of marginal distributions, including the use of mixture formulations to account for varying numbers of diseased and non-diseased individuals across studies. These components are combined to obtain a full likelihood-based model for meta-analysis.

Finally, we discuss how the model can be fitted in practice and how data can be simulated from the fitted model.

Model specification

The proposed approach is based on the representation of diagnostic test accuracy data as interval-censored observations. Consider a meta-analysis of \(I\) studies, where each study \(i = 1, \dots, I\) reports results for \(K_i\) diagnostic thresholds. These thresholds partition the observed data into \(K_i + 1\) disjoint regions.

For each study, individuals are divided into non-diseased and diseased populations. Let \(Y_h\) and \(Y_d\) denote the (unobserved) diagnostic test values for non-diseased and diseased individuals, respectively. Rather than observing exact values, the available data only provide information on how many individuals fall above or below certain thresholds. This induces censoring.

Each threshold defines bounds on the possible values of \((y_h, y_d)\). For a given censoring region \(j = 1, \dots, K_i + 1\) in study \(i\), we define lower and upper bounds

\[B_{hlij}, \; B_{huij}, \quad B_{dlij}, \; B_{duij},\] with \[B_{hlij} < B_{huij} ;B_{dlij} < B_{duij}.\]

Depending on the position of the region relative to the thresholds, three types of censoring occur:

Left-censoring: \(X_h \le B_{huij}\) and \(X_d \le B_{duij}\)
Right-censoring: \(X_h \ge B_{hlij}\) and \(X_d \ge B_{dlij}\)
Interval-censoring: \(B_{hlij} < X_h < B_{huij}\) and \(B_{dlij} < X_d < B_{duij}\)

Thus, each study contributes a set of bivariate censoring regions, rather than exact observations. These regions form the basis for constructing the likelihood of the model.

Model formulation

The underlying test values for non-diseased and diseased individuals are assumed to follow parametric distributions with CDF \(F_h(y)\) and \(F_d(y)\), and corresponding survival functions \(S_h(y) = 1 - F_h(y)\) and \(S_d(y) = 1 - F_d(y)\).

Commonly the Weibull distribution, parameterized as

\[ F_W(y; b, \lambda) = 1 - \exp\left( - (y \exp(-b))^{1/\lambda} \right), \] is used, where \(b\) and \(\lambda\) denote transformed scale and shape parameters. Separate parameters \((b_h, \lambda_h)\) and \((b_d, \lambda_d)\) are used for the non-diseased and diseased populations.

Since the observed data consist of counts of individuals above or below thresholds, we model these counts directly. Let \(H_i\) and \(D_i\) denote the number of non-diseased and diseased individuals in study \(i\).For a given threshold, the number of individuals exceeding the threshold is modeled (omitting the study index for simplicity) as

\[ X_h \sim \text{Bin}(H, S_h(y_h)), \] \[ X_d \sim \text{Bin}(D, S_d(y_d)).\]

These binomial distributions depend on the survival functions of the underlying test value distributions and thus link the individual-level model to the observed aggregated data. To obtain continuous marginals suitable for copula modeling, these distributions can be approximated by normal distributions:

\[X_h \sim \mathcal{N}(H S_h(y_h), \; H S_h(y_h) F_h(y_h)),\] \[X_d \sim \mathcal{N}(D S_d(y_d), \; D S_d(y_d) F_d(y_d)).\]

The joint distribution of \((y_h, y_d)\) is then constructed using a copula. By Sklar’s theorem, the joint distribution can be expressed as

\[ F(y_h, y_d) = C(F_h(y_h), F_d(y_d)), \]

where \(C(\cdot, \cdot)\) is a bivariate copula function. In this framework, a survival copula \(\hat{C}\) is used to model the joint survival function

\[ S(y_h, y_d) = \hat{C}(S_h(y_h), S_d(y_d)).\]

As an example, a Clayton copula can be used, which introduces a dependence parameter \(\theta \ge 0\) controlling the strength of association between the two populations. The copula allows flexible modeling of the dependence structure independently of the marginal distributions.

By combining the marginal models for the counts with the copula-based dependence structure, the model provides a flexible framework for jointly describing the distribution of diagnostic test results in diseased and non-diseased individuals across multiple studies and thresholds.

Likelihood

The full log-likelihood of the model is obtained by combining contributions from all censoring regions across studies. Let \(c_{ij} \in \mathbb{R}^3\) denote the likelihood contribution of region \(j\) in study \(i\), corresponding to the three possible censoring types (left, interval, right):

\[ c_{ij} = \begin{pmatrix} \max(C(F_{ijhu}, F_{ijdu}), \epsilon) \\ \max\big(\hat{C}(S_{ijhl}, S_{ijdl}) - \hat{C}(S_{ijhu}, S_{ijdl}) - \hat{C}(S_{ijhl}, S_{ijdu}) + \hat{C}(S_{ijhu}, S_{ijdu}), \epsilon \big) \\ \max(\hat{C}(S_{ijhl}, S_{ijdl}), \epsilon) \end{pmatrix}, \]

where \(\epsilon > 0\) is a small constant to prevent numerical underflow. The first, second, and third entries correspond to left-, interval-, and right-censoring, respectively. The indicator vector \(z_{ij} \in \{0,1\}^3\) selects the appropriate component depending on the censoring type.

The full log-likelihood is then:

\[\ell(x_1,x_2|\theta, b, \lambda, B) = \sum_{i=1}^I \sum_{j=1}^{K_i + 1} \log(z_{ij}^\top c_{ij}),\]

where \(B\) represents all lower and upper bounds of the censoring intervals, \(b = (b_h, b_d)\) and \(\lambda = (\lambda_h, \lambda_d)\) are the marginal distribution parameters, and \(\theta\) is the copula dependence parameter.

Optimization

Model parameters are estimated using numerical optimization. Two main strategies can be used:

One-step full optimization: All parameters \((b_h, \lambda_h, b_d, \lambda_d, \theta)\) are optimized simultaneously from suitable starting values.
Iterative three-step approach:
- Step 1: Estimate marginal parameters assuming independence.
- Step 2: Estimate copula parameter(s) with fixed marginals.
- Step 3: Re-optimize marginal parameters using the estimated copula dependence.

ROC construction

From the estimated marginal distributions and copula parameter, summary sensitivity and specificity at a threshold \(y\) are computed as:

\[\hat{\text{sens}}(y) = S_d(y \mid \hat{b}_d, \hat{\lambda}_d),\] \[\hat{\text{spec}}(y) = F_h(y \mid \hat{b}_h, \hat{\lambda}_h),\] where \(S_d\) is the survival function of the diseased population and \(F_h\) is the CDF of the non-diseased population, evaluated using the estimated distributional parameters \((\hat{b}_d, \hat{\lambda}_d)\) and \((\hat{b}_h, \hat{\lambda}_h)\), respectively.

These estimates form the summary ROC curve by plotting \(1 - \hat{\text{spec}}(y)\) on the x-axis and \(\hat{\text{sens}}(y)\) on the y-axis. The area under the curve (AUC) is then computed numerically using the trapezoidal rule on the estimated points, providing a single summary measure of overall diagnostic accuracy.

Model application in MetaROC

Inside the package, the model can be fitted using either the fit_metaROC() function or the metaROC.metaROC() method by setting action = "estimate". To fit the model, set model = "stoye2024copula" and specify the type of copula with copula_kind, which can be "Joe" or "Clayton".

Users can also choose the marginal distributions of the latent parameters via marginal_kind, selecting from "weibull", "gompertz", "loglogistic", or "lognormal". Additional options include using a normal approximation for the likelihood (norm_approx) or fixing the copula parameters during estimation (fix_cop_par).

We will demonstrate how to fit both types of copulas. Since the Joe copula exhibits convergence issues on this dataset when using the same hyperparameters as the Clayton copula, we will avoid using a normal approximation and will not use a lognormal marginal for fitting:

library(metaROC)
set.seed(7)
data(hba1c)
fit_clay <- fit_metaROC(hba1c, model = "stoye2024copula", copula_kind = "clayton", norm_approx = TRUE,
                            marginal_kind = "weibull")

## Hello and welcome to metaROC!
## Requested model: stoye2024copula 
## This is a copula model for multiple thresholds.
##  See XX for more details.

fit_joe <- fit_metaROC(hba1c, model = "stoye2024copula", copula_kind = "joe", norm_approx = FALSE,
                            marginal_kind = "lognormal")

## Hello and welcome to metaROC!
## Requested model: stoye2024copula 
## This is a copula model for multiple thresholds.
##  See XX for more details.

Inside the fit_metaROC_stoye2024copula() function, the copula model is fitted as follows. The main parameter estimation is performed with stats::optim(), which numerically maximizes the full log-likelihood defined by the copula and marginal distributions. Depending on the settings, the copula parameters can be fixed (fix_cop_par = TRUE) or estimated jointly with the marginals, and a normal approximation (norm_approx) can optionally be applied for the mixture marginals. After optimization, standard errors are computed from the Hessian matrix using numDeriv::hessian() and HelpersMG::SEfromHessian(), providing uncertainty estimates for the fitted parameters.

The function returns an S3 object containing the fitted model, convergence status, estimated parameters with standard errors, the thresholds included in the study data, and metadata on the copula and marginal specification. These components can be used directly to construct the summary ROC curve and derive predicted sensitivity and specificity values across thresholds.

summary(fit_clay, ci_type = "wald")

## 
## *** Results of Multiple threshold method (MTM) ***
## 
## Model: Copula model for multiple thresholds 
## Copula kind: clayton 
## Marignal distribution: weibull 
## Normal approximation: TRUE 
## Fixed copula parameters during estimation: FALSE 
## 
## Total number of studies: 38 
## Total number of thresholds: 124 
## Number of different thresholds: 26 
## 
## Results with Wald confidence intervals: 
## 
## Youden index (sensitivity weight = 0.5): 0.5873
## Optimal threshold value: 5.8549
## 
## Estimated Sensitivity and Specificity [95% CI]:
##  Sens: 0.6939 [0.6838; 0.7037]
##  Spec: 0.8934 [0.8870; 0.8997]
## 
## AUC: 0.8559

Firstly, the model name is returned, along with a reminder that the copula model is a multiple threshold method. Additionally, we gain the specific information about the copula we fitted. Afterwards, the summary provides a basic overview of the meta-analysis. The output then provides the largest Youden index together with the threshold at which it is achieved.

Following the estimation of the model, the summary reports the estimated sensitivity and specificity for the diagnostic test together with their respective 95% confidence intervals. Finally, the AUC is displayed, providing a summary measure of the overall diagnostic accuracy.

summary(fit_joe, ci_type = "wald")

## 
## *** Results of Multiple threshold method (MTM) ***
## 
## Model: Copula model for multiple thresholds 
## Copula kind: joe 
## Marignal distribution: lognormal 
## Normal approximation: FALSE 
## Fixed copula parameters during estimation: FALSE 
## 
## Total number of studies: 38 
## Total number of thresholds: 124 
## Number of different thresholds: 26 
## 
## Results with Wald confidence intervals: 
## 
## Youden index (sensitivity weight = 0.5): 0.6035
## Optimal threshold value: 5.7578
## 
## Estimated Sensitivity and Specificity [95% CI]:
##  Sens: 0.7398 [0.0000; 0.9883]
##  Spec: 0.8637 [0.0098; 1.0000]
## 
## AUC: 0.8738

Looking at the summary for the Joe copula, we can see that a different copula type as well as an alternative marginal distribution is returned. While the optimal threshold and Youden index are quite similar, the estimated sensitivity and specificity differ substantially, and the confidence intervals are much wider, almost spanning the entire range from 0 to 1, indicating a high degree of uncertainty in the estimates.

Furthermore, the Joe copula yields a slightly higher AUC, but the estimated sensitivity and specificity are highly uncertain. In contrast, the Clayton copula provides more stable and precise estimates, making it the more reliable choice for this dataset despite the marginally lower AUC.

SROC <- SROC(fit_clay)
head(SROC$sroc_df, 10)

##    threshold   sens_lo sensitivity   sens_up     spec_lo specificity
## 1   3.510000 0.9984219   0.9984591 0.9984952 0.006948698 0.006838290
## 2   3.514855 0.9983983   0.9984362 0.9984729 0.007057000 0.006945670
## 3   3.519710 0.9983745   0.9984130 0.9984503 0.007166830 0.007054579
## 4   3.524565 0.9983503   0.9983894 0.9984273 0.007278208 0.007165035
## 5   3.529419 0.9983258   0.9983656 0.9984041 0.007391154 0.007277059
## 6   3.534274 0.9983010   0.9983414 0.9983806 0.007505687 0.007390670
## 7   3.539129 0.9982758   0.9983169 0.9983567 0.007621827 0.007505888
## 8   3.543984 0.9982503   0.9982921 0.9983326 0.007739595 0.007622734
## 9   3.548839 0.9982245   0.9982670 0.9983081 0.007859010 0.007741229
## 10  3.553694 0.9981983   0.9982415 0.9982833 0.007980093 0.007861392
##        spec_up youden_index
## 1  0.006729679  0.005297377
## 2  0.006836146  0.005381845
## 3  0.006944142  0.005467534
## 4  0.007053685  0.005554458
## 5  0.007164796  0.005642635
## 6  0.007277494  0.005732079
## 7  0.007391800  0.005822808
## 8  0.007507734  0.005914837
## 9  0.007625317  0.006008184
## 10 0.007744570  0.006102865

The summary ROC curve is computed as described in the model specification section. For each threshold, the estimated sensitivity and specificity are obtained from the fitted model parameters using the marginal survival and cumulative distribution functions. Confidence intervals for both sensitivity and specificity are derived by propagating the standard errors of the estimated parameters through the same functions. These estimates are combined into a single data frame, and the Youden index is calculated to assess the optimal balance between sensitivity and specificity. Finally, the area under the SROC curve (AUC) is computed numerically from the summary ROC points, providing an overall measure of diagnostic accuracy across all thresholds.

Simulation

Simulation of datasets from the copula models can be performed by supplying the model-specific parameters in a list. For example, model_pars can contain beta00 and lambda0 specifying the marginal distribution parameters for the non-diseased population, beta01 and lambda1 for the diseased population, copula_kind_sim defining the copula type ("clayton" or "joe"), cop_params providing the copula dependence parameters (a single value for Clayton or a vector of three for Joe), norm_approx_sim indicating whether a normal approximation is used for the mixture marginals, and marginal_kind_sim specifying the chosen marginal distribution.

The main difference between the copula kinds when simulating is that Clayton requires only one copula parameter, whereas Joe requires three. We will demonstrate simulation for both copula kinds in the following examples.

control_sim_clay <- list(n_studies_min = 5, n_studies_max = 10,
                        n_subj_min = 50, n_subj_max = 200,
                        n_thresholds_min = 2, n_thresholds_max = 5,
                        thresholds_lower_bound = 0.1, thresholds_upper_bound = 0.9,
                        true_preval_min = 0.1,  true_preval_max = 0.5,
                        model = "stoye2024copula",
                        model_pars = list(beta00 = 0.5, lambda0 = 1.2,
                                          beta01 = -0.3, lambda1 = 0.8,
                                          copula_kind_sim = "clayton",
                                          cop_params = 2,
                                          norm_approx_sim = TRUE,
                                          marginal_kind_sim = "lognormal"
  )
)

sim_cop_clay <- simulate_metaROC(type = "phaseII", n_sims = 5,
                 control_sim = control_sim_clay)
sim_cop_clay

## $data
## # A tibble: 118 × 7
##    replicate study threshold    TP    FN    FP    TN
##        <int> <int>     <dbl> <dbl> <dbl> <dbl> <dbl>
##  1         1     1       0.4    13     5   115    19
##  2         1     1       0.5    13     5   115    19
##  3         1     1       0.7     9     9   103    31
##  4         1     1       0.8     9     9    98    36
##  5         1     2       0.4    23     8    29     4
##  6         1     2       0.7    18    13    26     7
##  7         1     2       0.8    12    19    23    10
##  8         1     3       0.2    23     0   139     5
##  9         1     3       0.4    16     7   120    24
## 10         1     3       0.8    15     8   114    30
## # ℹ 108 more rows
## 
## $data_info
## $data_info$model_id
## [1] "stoye2024copula"
## 
## $data_info$true_sens
## [1] 0.9938467 0.9491636 0.8707547 0.7794578 0.6884402 0.6039297 0.5282389
## [8] 0.4617323 0.4038868
## 
## $data_info$true_spec
## [1] 0.009758982 0.039385247 0.077807140 0.118951934 0.160041018 0.199794725
## [7] 0.237645811 0.273380881 0.306966589
## 
## $data_info$true_AUC
## [1] 0.3568727
## 
## $data_info$model_pars
## $data_info$model_pars$beta00
## [1] 0.5
## 
## $data_info$model_pars$lambda0
## [1] 1.2
## 
## $data_info$model_pars$beta01
## [1] -0.3
## 
## $data_info$model_pars$lambda1
## [1] 0.8
## 
## $data_info$model_pars$copula_kind_sim
## [1] "clayton"
## 
## $data_info$model_pars$cop_params
## [1] 2
## 
## $data_info$model_pars$norm_approx_sim
## [1] TRUE
## 
## $data_info$model_pars$marginal_kind_sim
## [1] "lognormal"
## 
## $data_info$model_pars$all_available_thresholds
## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

control_sim_joe <- list(n_studies_min = 5, n_studies_max = 10,
                        n_subj_min = 50, n_subj_max = 200,
                        n_thresholds_min = 2, n_thresholds_max = 5,
                        thresholds_lower_bound = 0.1, thresholds_upper_bound = 0.9,
                        true_preval_min = 0.1,  true_preval_max = 0.5,
                        model = "stoye2024copula",
                        model_pars = list(beta00 = 0.5, lambda0 = 1.2,
                                          beta01 = -0.3, lambda1 = 0.8,
                                          copula_kind_sim = "joe",
                                          cop_params = c(1,1,2),
                                          norm_approx_sim = TRUE,
                                          marginal_kind_sim = "lognormal"
  )
)

sim_cop_joe <- simulate_metaROC(type = "phaseII", n_sims = 5,
                 control_sim = control_sim_joe)

Computationally Simulation is done as follows:

Simulation from the copula models proceeds by first generating pairs of uniform random variables from the specified copula. For the Clayton copula, this is done using copula::rCopula() with the copula parameter provided, whereas the Joe copula requires a more complex acceptance-rejection sampling procedure to account for its three copula parameters. These uniform samples are then transformed into counts of true positives and false positives at each threshold using the inverse of the chosen marginal survival distributions, either exactly via binomial quantiles or approximately via normal approximations for mixture marginals. A failsafe ensures that the counts remain consistent across thresholds, preserving the monotonicity of cumulative events.

To gain an overview of how to plot a fitted model and how to conduct simulation studies, particularly for evaluating models such as the copula (Stoye) model, please refer to the other vignettes included in this package, which provide more detailed guidance on these topics.

2026-04-28