Meta-Regression

Subgroup Analysis

This vignette introduces the regression functionality of the dtametaTMB package for meta-analysis of diagnostic test accuracy (DTA) studies.

Validation

The subgroup HSROC and Reitsma implementations reproduce the Cochrane Handbook RF and Anti-CCP examples closely, including the baseline parameters (HSROC: accuracy, threshold, shape; Reitsma: logit sensitivity and specificity), between-study variance estimates, and subgroup effects on accuracy and threshold. The Schuetz CT/MRI analyses yielded results broadly consistent with the published subgroup parameter estimates and variance components. Likelihood-ratio tests led to identical substantive conclusions, although the LR statistics differed somewhat from the published values.

Dummy-coding vs. cell-means parameterization

For subgroup analyses, we distinguish between two equivalent parameterizations of the same model. In the reference-group parameterization, one subgroup serves as the baseline and the remaining subgroup effects are expressed as deviations from this reference. In the group-specific (cell-means) parameterization, each subgroup has its own parameter directly. The former is convenient for testing subgroup differences, whereas the latter is convenient for reporting subgroup-specific estimates and plotting subgroup-specific HSROC curves. The Reitsma subgroup model is fitted twice with both parameterizations. The subgroup HSROC model is estimated using a reference-group parameterization, and subgroup-specific parameters are recovered by evaluating the fitted linear predictors at subgroup-specific design points via Zpred.

How do we include covariates in the Reitsma model?

For sensitivity in study i, we have the number of diseased individuals testing positive: yAi ∼ ℬ(nAi, πAi).

Similarly for specificity, we have the number of non-diseased individuals testing negative: yBi ∼ ℬ(nBi, πBi).

Now we introduce a p-dimensional design-vector zi including study level covariates. Consequently, at the study level, we have

$$ \begin{pmatrix} \boldsymbol{z}_i^{\top}\boldsymbol{\mu}_{Ai} \\ \boldsymbol{z}_i^{\top}\boldsymbol{\mu}_{Bi} \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \boldsymbol{z}_i^{\top}\boldsymbol{\mu}_A \\ \boldsymbol{z}_i^{\top}\boldsymbol{\mu}_B \end{pmatrix}, \; \Sigma \right), \quad \text{with} \quad \Sigma = \begin{pmatrix} \sigma_A^2 & \sigma_{AB} \\ \sigma_{AB} & \sigma_B^2 \end{pmatrix}.$$

Let’s assume that zi includes a single binary covariate representing two subgroups. Then using dummy coding with subgroup 1 being the reference, we have

$$\boldsymbol{z}_i^{\top} \boldsymbol{\mu}_{Ai} = \begin{cases} \mu_{Ai} & \text{for subgroup 1}, \\ \mu_{Ai} + \nu_{A2} & \text{for subgroup 2}, \end{cases} \quad \quad \boldsymbol{z}_i^{\top} \boldsymbol{\mu}_{Bi} = \begin{cases} \mu_{Bi} & \text{for subgroup 1}, \\ \mu_{Bi} + \nu_{B2} & \text{for subgroup 2}, \end{cases} $$

$$ \boldsymbol{z}_i^{\top} \boldsymbol{\mu}_{A} = \begin{cases} \mu_{A} & \text{for subgroup 1}, \\ \mu_{A} + \nu_{A2} & \text{for subgroup 2}, \end{cases} \quad \quad \boldsymbol{z}_i^{\top} \boldsymbol{\mu}_B = \begin{cases} \mu_{B} & \text{for subgroup 1}, \\ \mu_{B} + \nu_{B2} & \text{for subgroup 2}. \end{cases} $$

data("anticcp")
reitsmasub <- fitReitsmaSubgroup(data=anticcp,
                                 TP=TP,
                                 FP=FP,
                                 FN=FN,
                                 TN=TN,
                                 study=study,
                                 subgroup=generation)
reitsmasub
#> 
#> Reitsma Subgroup Model
#> ----------------------
#> 
#> Number of studies   : 37 
#> Number of subgroups : 2 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0004844247 
#> -2 log likelihood   : 533.37 ( df = 7 )
#> AIC                 : 547.37 
#> BIC                 : 558.646 
#> 
#> Use summary() for parameter estimates.
summary(reitsmasub)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                               Estimate   Std_Error
#> mu_A.CCP1                 -0.096538884 0.220320590
#> mu_B.CCP1                  3.446719200 0.298243690
#> mu_A.CCP2                  0.866038480 0.120876820
#> mu_B.CCP2                  3.016490600 0.162237540
#> sigma2_A.sens.log_sigma_A  0.359828710 0.102176500
#> sigma2_B.spec.log_sigma_B  0.539909310 0.180157230
#> sigma_AB                  -0.196844990 0.098353424
#> nu_A.CCP2                  0.962571630 0.251342040
#> nu_B.CCP2                 -0.430210700 0.337711820
#> 
#> Sensitivity / specificity
#> -------------------------
#>           type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.CCP1 sens 0.47588401      0.95 0.37089964 0.58304390
#> mu_B.CCP1 spec 0.96913315      0.95 0.94594449 0.98255779
#> mu_A.CCP2 sens 0.70392072      0.95 0.65229087 0.75081294
#> mu_B.CCP2 spec 0.95331358      0.95 0.93693872 0.96559260
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>         Lambda       Theta       beta sigma2_alpha sigma2_theta
#> CCP1 3.0073766 -1.61053440 0.20288653   0.48784248   0.31880561
#> CCP2 3.6840001 -0.88349711 0.20288653   0.48784248   0.31880561
#> 
#> Subgroups
#> ---------
#> [1] "CCP1" "CCP2"

How do I get a summary plot of the Reitsma model?

plot(reitsmasub,
     scale=0.01,
     nudge_legend=-0.2,
     size="se",
     col=c("black","red"))

How do I get a coupled forest plot?

Note: Rendering forest plots may take longer in an interactive R session due to the underlying grid graphics. Performance is typically faster when knitting the vignette (e.g. via RMarkdown or Quarto).

forest(reitsmasub,subgroup_label="Generation")

How do I constrain parameters in the Reitsma model?

In sparse data one may wish to fix parameters of the random effects (variance-covariance matrix) at zero. This can be done via the constrain argument.

For example,

constrainA <- fitReitsmaSubgroup(data=anticcp,
                                 TP=TP,
                                 FP=FP,
                                 FN=FN,
                                 TN=TN,
                                 study=study,
                                 subgroup=generation,
                                 constrain="sigma2_A")
summary(constrainA)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                                 Estimate   Std_Error
#> mu_A.CCP1                 -5.4394046e-02 0.054985303
#> mu_B.CCP1                  3.4466250e+00 0.298754820
#> mu_A.CCP2                  9.1797945e-01 0.031391722
#> mu_B.CCP2                  2.9968587e+00 0.162162330
#> sigma2_A.sens.log_sigma_A  4.9303807e-32 0.000000000
#> sigma2_B.spec.log_sigma_B  5.3963306e-01 0.179976400
#> sigma_AB                   0.0000000e+00 0.000000000
#> nu_A.CCP2                  9.7237352e-01 0.063315273
#> nu_B.CCP2                 -4.4976608e-01 0.337953820
#> 
#> Sensitivity / specificity
#> -------------------------
#>           type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.CCP1 sens 0.48640484      0.95 0.45954779 0.51334062
#> mu_B.CCP1 spec 0.96913033      0.95 0.94588841 0.98257334
#> mu_A.CCP2 sens 0.71463023      0.95 0.70191901 0.72701020
#> mu_B.CCP2 spec 0.95243201      0.95 0.93577764 0.96492938
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>          Lambda      Theta     beta sigma2_alpha  sigma2_theta
#> CCP1 -3128642.1 -1564321.1 35.73522 3.262267e-16 8.1556674e-17
#> CCP2 52800433.0 26400216.0 35.73522 3.262267e-16 8.1556674e-17
#> 
#> Subgroups
#> ---------
#> [1] "CCP1" "CCP2"

fixes the logit sensitivity variance to zero. This also implies that the random effects covariance is zero. Note that you can also set constrain to "sigma_AB", "sigma2_B", or "all".

Constraining fixed effects is controlled by the sensspec_constrain argument. For example,

constrainsens <- fitReitsmaSubgroup(data=anticcp,
                                    TP=TP,
                                    FP=FP,
                                    FN=FN,
                                    TN=TN,
                                    study=study,
                                    subgroup=generation,
                                    sensspec_constrain="sens")
summary(constrainsens)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                               Estimate  Std_Error
#> mu_A.CCP1                  0.653305180 0.12740877
#> mu_B.CCP1                  3.081222100 0.32561613
#> mu_A.CCP2                  0.653305180 0.00000000
#> mu_B.CCP2                  3.118005800 0.17451227
#> sigma2_A.sens.log_sigma_A  0.541939800 0.14627725
#> sigma2_B.spec.log_sigma_B  0.577648680 0.19968290
#> sigma_AB                  -0.277162040 0.13981590
#> nu_A.CCP2                  0.000000000 0.00000000
#> nu_B.CCP2                  0.036831361 0.38577185
#> 
#> Sensitivity / specificity
#> -------------------------
#>           type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.CCP1 sens 0.65775489      0.95 0.59954955 0.71157008
#> mu_B.CCP1 spec 0.95611150      0.95 0.92004997 0.97632597
#> mu_A.CCP2 sens 0.65775489      0.95 0.59954955 0.71157008
#> mu_B.CCP2 spec 0.95762939      0.95 0.94136328 0.96952925
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>         Lambda      Theta        beta sigma2_alpha sigma2_theta
#> CCP1 3.6962689 -1.1843237 0.031905471   0.56469478   0.41833574
#> CCP2 3.7324705 -1.2024245 0.031905471   0.56469478   0.41833574
#> 
#> Subgroups
#> ---------
#> [1] "CCP1" "CCP2"

assumes equal (logit) sensitivities in all subgroups. Note that you can also set sensspec_constrain to "spec" or c("sens","spec").

How can I compare constrainsens with the full model?

We can perform a likelihood ratio test via anova, essentially testing whether there exist subgroup differences in sensitivity.

anova(constrainsens,reitsmasub)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1  6 -272.77                              
#> Model 2  7 -266.69       1 12.181  0.0004827 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How do I allow for different subgroup-specific random-effects (co-)variances in the Reitsma model?

heteroskedastic <- fitReitsmaSubgroup(data=anticcp,
                                      TP=TP,
                                      FP=FP,
                                      FN=FN,
                                      TN=TN,
                                      study=study,
                                      subgroup=generation,
                                      variances="unequal")
plot(heteroskedastic,
     scale=0.01,
     nudge_legend=-0.2,
     size="se",
     col=c("black","red"))

anova(reitsmasub,heteroskedastic)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)  
#> Model 1  7 -266.69                            
#> Model 2 10 -262.42       3 8.5306    0.03623 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How do I fit the Reitsma subgroup model using the custom TMB backend?

By default, fitReitsmaSubgroup() fits the Reitsma model using glmmTMB. The package also provides fitReitsmaSubgroupTMB(), which fits the same binomial-normal model using a custom Template Model Builder (TMB) template.

The custom TMB backend may offer improved numerical robustness in difficult or near-boundary cases.

data("anticcp")
reitsmasubTMB <- fitReitsmaSubgroupTMB(data=anticcp,
                                       TP=TP,
                                       FP=FP,
                                       FN=FN,
                                       TN=TN,
                                       study=study,
                                       variances="unequal",
                                       subgroup=generation)
reitsmasubTMB
#> 
#> Reitsma Subgroup TMB Model
#> --------------------------
#> 
#> Number of studies   : 37 
#> Number of subgroups : 2 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 7.201756e-05 
#> -2 log likelihood   : 524.839 ( df = 10 )
#> AIC                 : 544.839 
#> BIC                 : 560.948 
#> 
#> Use summary() for parameter estimates.
summary(reitsmasubTMB)
#> 
#> Summary: Reitsma Subgroup TMB Model
#> -----------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                   Estimate  Std. Error
#> mu_A.CCP1     -0.077734976 0.093745405
#> mu_A.CCP2      0.872414330 0.137541560
#> mu_B.CCP1      3.458614000 0.290964660
#> mu_B.CCP2      3.014784100 0.164919640
#> sigma2_A.CCP1  0.042771854 0.033050876
#> sigma2_A.CCP2  0.481452150 0.152370450
#> sigma2_B.CCP1  0.479322870 0.320124020
#> sigma2_B.CCP2  0.560479850 0.215712700
#> sigma_AB.CCP1 -0.062584493 0.079018858
#> sigma_AB.CCP2 -0.233390560 0.129708840
#> rho_AB.CCP1   -0.437092770 0.445469900
#> rho_AB.CCP2   -0.449289960 0.200921540
#> nu_A.CCP2      0.950149300 0.166450840
#> nu_B.CCP2     -0.443829930 0.334453170
#> 
#> Sensitivity / specificity
#> -------------------------
#>           type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.CCP1 sens 0.48057604      0.95 0.43500175 0.52647587
#> mu_A.CCP2 sens 0.70524782      0.95 0.64630529 0.75804497
#> mu_B.CCP1 spec 0.96948699      0.95 0.94726671 0.98251710
#> mu_B.CCP2 spec 0.95323757      0.95 0.93652604 0.96571036
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>         Lambda       Theta        beta sigma2_alpha sigma2_theta
#> CCP1 1.7480891 -1.01627220 1.208247100   0.16119811   0.10288402
#> CCP2 3.8085824 -0.99809026 0.075993224   0.57214955   0.37642794
#> 
#> Subgroups
#> ---------
#> [1] "CCP1" "CCP2"
plot(reitsmasubTMB,
     scale=0.01,
     nudge_legend=-0.2,
     size="se",
     col=c("black","red"))

How do we include covariates in the Rutter and Gatsonis model?

The number of diseased individuals from study i who test positive is denoted by yi1 ∼ ℬ(ni1, πi1).

Similarly, the number of non-diseased individuals who test positive is yi2 ∼ ℬ(ni2, πi2).

Now we introduce a p-dimensional design-vector zi including study level covariates. Consequently, at the study level, we have

logit (πij) = (zi⊤θi + zi⊤αixij)exp (−zi⊤βxij),

zi⊤αi ∼ 𝒩(zi⊤Λ, σα2),    zi⊤θi ∼ 𝒩(zi⊤Θ, σθ2),

where $$x_{ij} = \begin{cases} -0.5 & \text{for non-diseased individuals}, \\ \phantom{-}0.5 & \text{for diseased individuals}. \end{cases}$$

Let’s assume that zi includes a single categorical covariate representing three subgroups. Then using dummy coding with subgroup 1 being the reference, we have

$$\boldsymbol{z}_i^{\top} \boldsymbol{\theta}_i = \begin{cases} \theta_i & \text{for subgroup 1}, \\ \theta_i + \gamma_2 & \text{for subgroup 2},\\ \theta_i + \gamma_3 & \text{for subgroup 3}, \end{cases} \quad \quad \boldsymbol{z}_i^{\top} \boldsymbol{\alpha}_i = \begin{cases} \alpha_i & \text{for subgroup 1}, \\ \alpha_i + \xi_2 & \text{for subgroup 2},\\ \alpha_i + \xi_3 & \text{for subgroup 3}, \end{cases} $$ $$\boldsymbol{z}_i^{\top} \boldsymbol{\beta} = \begin{cases} \beta & \text{for subgroup 1}, \\ \beta + \delta_2 & \text{for subgroup 2},\\ \beta + \delta_3 & \text{for subgroup 3}, \end{cases}$$ $$\boldsymbol{z}_i^{\top} \boldsymbol{\Lambda} = \begin{cases} \Lambda & \text{for subgroup 1}, \\ \Lambda + \xi_2 & \text{for subgroup 2},\\ \Lambda + \xi_3 & \text{for subgroup 3}, \end{cases} \quad \quad \boldsymbol{z}_i^{\top} \boldsymbol{\Theta}_i = \begin{cases} \Theta & \text{for subgroup 1}, \\ \Theta + \gamma_2 & \text{for subgroup 2},\\ \Theta +\gamma_3 & \text{for subgroup 3}. \end{cases}$$

Note: For prediction fitRutterGatsonisSubgroup() uses the prediction matrix $$\boldsymbol{Z}_{\mathrm{pred}}=\begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix}$$ for the case above and therefore recovers the threshold, accuracy, and shape parameters for each subgroup as if it were the reference group.

How do I perform a Rutter and Gatsonis meta-regression with a single categorical study-level covariate?

data("RF")
RF2        <- RF[RF$method %in% c("LA","ELISA","Nephelometry"),]
RF2$method <- factor(RF2$method,levels=c("LA","ELISA","Nephelometry"))
ruttergatsonissub <- fitRutterGatsonisSubgroup(data=RF2,
                                               TP=TP,
                                               FP=FP,
                                               FN=FN,
                                               TN=TN,
                                               study=study,
                                               subgroup=method,
                                               constrain="shape") # assumes equal 
                                                                  # shapes in subgroups
ruttergatsonissub
#> 
#> Rutter & Gatsonis Subgroup Model
#> --------------------------------
#> 
#> Number of studies   : 47 
#> Number of subgroups : 3 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0002867077 
#> -2 log likelihood   : 753.722 ( df = 9 )
#> AIC                 : 771.722 
#> BIC                 : 788.374 
#> 
#> Use summary() for parameter estimates.
summary(ruttergatsonissub)
#> 
#> Summary: Rutter & Gatsonis Subgroup Model
#> -----------------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                         Estimate  Std. Error
#> Lambda_LA            2.456309800 0.323577140
#> xi_ELISA             0.248533990 0.439541130
#> xi_Nephelometry      0.331428920 0.442600640
#> Theta_LA            -0.550060260 0.213340320
#> gamma_ELISA         -0.196255960 0.260961790
#> gamma_Nephelometry   0.495863020 0.262230320
#> beta_LA              0.197677370 0.169839730
#> delta_ELISA          0.000000000 0.000000000
#> delta_Nephelometry   0.000000000 0.000000000
#> sigma2_alpha         1.277914900 0.308524820
#> sigma2_theta         0.476839640 0.113444950
#> Lambda_LA            2.456309800 0.323577140
#> Lambda_ELISA         2.704843800 0.326954560
#> Lambda_Nephelometry  2.787738800 0.305779190
#> Theta_LA            -0.550060260 0.213340320
#> Theta_ELISA         -0.746316220 0.209913080
#> Theta_Nephelometry  -0.054197236 0.212131240
#> beta_LA              0.197677370 0.169839730
#> beta_ELISA           0.197677370 0.169839730
#> beta_Nephelometry    0.197677370 0.169839730
#> logitsens            0.826168420 0.286300410
#> logitsens            1.051312600 0.286070100
#> logitsens            1.126406200 0.276895840
#> sens                 0.695544150 0.060627687
#> sens                 0.741026870 0.054898582
#> sens                 0.755175060 0.051194078
#> 
#> Sensitivity / specificity
#> -------------------------
#>       subgroup       spec conflevel  logitsens  Std_Error   CI_Lower  CI_Upper
#> 1           LA 0.84615385      0.95 0.82616842 0.28630041 0.26502992 1.3873069
#> 2        ELISA 0.84615385      0.95 1.05131260 0.28607010 0.49062549 1.6119997
#> 3 Nephelometry 0.84615385      0.95 1.12640620 0.27689584 0.58370030 1.6691121
#>         Sens SensCI_Lower SensCI_Upper
#> 1 0.69554415   0.56587235   0.80016196
#> 2 0.74102687   0.62025377   0.83368883
#> 3 0.75517506   0.64191840   0.84145740
#> 
#> Recovered bivariate parameters (Reitsma)
#> ----------------------------------------
#>               mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec    sigma_AB
#> LA           0.61427837 1.9629508    0.65348637    0.97036902 -0.15736092
#> ELISA        0.54906437 2.3167723    0.65348637    0.97036902 -0.15736092
#> Nephelometry 1.21359400 1.5985037    0.65348637    0.97036902 -0.15736092
#> 
#> Subgroups
#> ---------
#> [1] "LA"           "ELISA"        "Nephelometry"

How do I get a summary plot of the Rutter and Gatsonis subgroup model?

plot(ruttergatsonissub,
     specrange=c(0.3,0.995),
     size="se",
     col=c("red","black","green"),
     scale=0.015)

How do I get a coupled forest plot?

Note: Rendering forest plots may take longer in an interactive R session due to the underlying grid graphics. Performance is typically faster when knitting the vignette (e.g. via RMarkdown or Quarto).

forest(ruttergatsonissub,subgroup_label = "Method")

How do I constrain parameters in the Rutter and Gatsonis subgroup model?

In the Rutter and Gatsonis model, all parameter constraints are controlled by constrain, i.e.,

Constraints can also be combined, for example constrain=c("shape","sigma2_theta").

How do I use the general Rutter and Gatsonis regression function?

This method is for advanced users who feel comfortable specifying their own design and prediction matrices. For study-level covariates, the design matrix Z needs two identical consecutive rows per study, one for the diseased and one for the non-diseased. Of note, there are neither plot() nor forest() methods for fitRutterGatsonisReg().

Let’s reproduce the subgroup-analysis from before.

# Specify design matrix Z
Z  <- model.matrix(~method,data=RF2)
# For study level-covariates, we need to two identical consecutive 
# rows per study (diseased and non-diseased).
Z2 <- Z[rep(seq_len(nrow(Z)), each = 2), , drop = FALSE]
# Specify prediction matrix Z_pred
Z_pred <- matrix(c(1,0,0,1,1,0,1,0,1),ncol=3,nrow=3,byrow=T)
constrain <- list(shape_coef=factor(c(1, rep(NA, ncol(Z2) - 1))))

ruttergatsonisreg <- fitRutterGatsonisReg(data=RF2,
                                          TP=TP,
                                          FP=FP,
                                          FN=FN,
                                          TN=TN,
                                          study=study,
                                          Z=Z2,
                                          Z_pred=Z_pred,
                                          map=constrain)

ruttergatsonisreg
#> 
#> Rutter & Gatsonis Regression Model
#> ----------------------------------
#> 
#> Number of studies : 47 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.0003981125 
#> -2 log likelihood : 753.722 ( df = 9 )
#> AIC               : 771.722 
#> BIC               : 788.374 
#> 
#> Use summary() for parameter estimates.
summary(ruttergatsonisreg)
#> 
#> Summary: Rutter & Gatsonis Regression Model
#> -------------------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                    Estimate  Std. Error
#> accuracy_coef   2.456306200 0.323578110
#> accuracy_coef   0.248536610 0.439542730
#> accuracy_coef   0.331419970 0.442602170
#> threshold_coef -0.550050810 0.213340320
#> threshold_coef -0.196254760 0.260962000
#> threshold_coef  0.495866140 0.262230550
#> shape_coef      0.197682730 0.169839700
#> shape_coef      0.000000000 0.000000000
#> shape_coef      0.000000000 0.000000000
#> sigma2_alpha    1.277925600 0.308528070
#> sigma2_theta    0.476840430 0.113445210
#> Lambda_Pred     2.456306200 0.323578110
#> Lambda_Pred     2.704842800 0.326955320
#> Lambda_Pred     2.787726100 0.305780160
#> Theta_Pred     -0.550050810 0.213340320
#> Theta_Pred     -0.746305570 0.209913120
#> Theta_Pred     -0.054184665 0.212131320
#> beta_Pred       0.197682730 0.169839700
#> beta_Pred       0.197682730 0.169839700
#> beta_Pred       0.197682730 0.169839700
#> logitsens       0.826166620 0.286300670
#> logitsens       1.051312600 0.286070290
#> logitsens       1.126395500 0.276896050
#> sens            0.695543770 0.060627783
#> sens            0.741026870 0.054898619
#> sens            0.755173080 0.051194397
#> 
#> Sensitivity / specificity
#> -------------------------
#>         spec conflevel  logitsens  Std_Error   CI_Lower  CI_Upper       Sens
#> 1 0.84615385      0.95 0.82616662 0.28630067 0.26502762 1.3873056 0.69554377
#> 2 0.84615385      0.95 1.05131260 0.28607029 0.49062509 1.6120000 0.74102687
#> 3 0.84615385      0.95 1.12639550 0.27689605 0.58368917 1.6691018 0.75517308
#>   SensCI_Lower SensCI_Upper
#> 1   0.56587179   0.80016175
#> 2   0.62025368   0.83368888
#> 3   0.64191584   0.84145602

How do I compare models?

In the previous section we fitted the RF data set using the Rutter and Gatsonis subgroup model while keeping the shape parameter equal across all subgroups constrain="shape". Now let’s fit the full model allowing for different shape parameters across subgroups and check whether the data lend support to this approach. In addition, let’s also consider a model that constrains both shape and accuracy.

ruttergatsonissubfull <- fitRutterGatsonisSubgroup(data=RF2,
                                                   TP=TP,
                                                   FP=FP,
                                                   FN=FN,
                                                   TN=TN,
                                                   study=study,
                                                   subgroup=method,
                                                   constrain=NULL)
ruttergatsonissubfull
#> 
#> Rutter & Gatsonis Subgroup Model
#> --------------------------------
#> 
#> Number of studies   : 47 
#> Number of subgroups : 3 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0001936155 
#> -2 log likelihood   : 753.553 ( df = 11 )
#> AIC                 : 775.553 
#> BIC                 : 795.905 
#> 
#> Use summary() for parameter estimates.
summary(ruttergatsonissubfull)
#> 
#> Summary: Rutter & Gatsonis Subgroup Model
#> -----------------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                        Estimate  Std. Error
#> Lambda_LA            2.42435650 0.330490060
#> xi_ELISA             0.29743443 0.513150050
#> xi_Nephelometry      0.37161351 0.450864220
#> Theta_LA            -0.50294636 0.244516610
#> gamma_ELISA         -0.25783318 0.370244910
#> gamma_Nephelometry   0.39706424 0.359443140
#> beta_LA              0.27775702 0.266422660
#> delta_ELISA         -0.10290405 0.426607590
#> delta_Nephelometry  -0.16039260 0.397537120
#> sigma2_alpha         1.27305060 0.307984420
#> sigma2_theta         0.47624252 0.113436600
#> Lambda_LA            2.42435650 0.330490060
#> Lambda_ELISA         2.72179090 0.392992730
#> Lambda_Nephelometry  2.79597000 0.307080030
#> Theta_LA            -0.50294636 0.244516610
#> Theta_ELISA         -0.76077954 0.278169440
#> Theta_Nephelometry  -0.10588212 0.263228630
#> beta_LA              0.27775702 0.266422660
#> beta_ELISA           0.17485297 0.333210110
#> beta_Nephelometry    0.11736442 0.295002910
#> logitsens            0.81868527 0.275188790
#> logitsens            1.06266640 0.324053890
#> logitsens            1.12065150 0.288348430
#> sens                 0.69395719 0.058444763
#> sens                 0.74319976 0.061846941
#> sens                 0.75410954 0.053467971
#> 
#> Sensitivity / specificity
#> -------------------------
#>       subgroup       spec conflevel  logitsens  Std_Error   CI_Lower  CI_Upper
#> 1           LA 0.84615385      0.95 0.81868527 0.27518879 0.27932515 1.3580454
#> 2        ELISA 0.84615385      0.95 1.06266640 0.32405389 0.42753242 1.6978003
#> 3 Nephelometry 0.84615385      0.95 1.12065150 0.28834843 0.55549896 1.6858041
#>         Sens SensCI_Lower SensCI_Upper
#> 1 0.69395719   0.56938077   0.79544184
#> 2 0.74319976   0.60528428   0.84524723
#> 3 0.75410954   0.63541045   0.84367155
#> 
#> Recovered bivariate parameters (Reitsma)
#> ----------------------------------------
#>               mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec    sigma_AB
#> LA           0.61726846 1.9706516    0.60182246    1.04887820 -0.15797987
#> ELISA        0.54987796 2.3155155    0.66705102    0.94631215 -0.15797987
#> Nephelometry 1.21846130 1.5947581    0.70652253    0.89344424 -0.15797987
#> 
#> Subgroups
#> ---------
#> [1] "LA"           "ELISA"        "Nephelometry"
ruttergatsonissubacc <- fitRutterGatsonisSubgroup(data=RF2,
                                                  TP=TP,
                                                  FP=FP,
                                                  FN=FN,
                                                  TN=TN,
                                                  study=study,
                                                  subgroup=method,
                                                  constrain=c("shape","accuracy"))
ruttergatsonissubacc
#> 
#> Rutter & Gatsonis Subgroup Model
#> --------------------------------
#> 
#> Number of studies   : 47 
#> Number of subgroups : 3 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0005592484 
#> -2 log likelihood   : 754.332 ( df = 7 )
#> AIC                 : 768.332 
#> BIC                 : 781.283 
#> 
#> Use summary() for parameter estimates.
summary(ruttergatsonissubacc)
#> 
#> Summary: Rutter & Gatsonis Subgroup Model
#> -----------------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                        Estimate  Std. Error
#> Lambda_LA            2.65810190 0.191815720
#> xi_ELISA             0.00000000 0.000000000
#> xi_Nephelometry      0.00000000 0.000000000
#> Theta_LA            -0.55746273 0.212932400
#> gamma_ELISA         -0.19386257 0.260935930
#> gamma_Nephelometry   0.49703130 0.262130440
#> beta_LA              0.18948853 0.164729720
#> delta_ELISA          0.00000000 0.000000000
#> delta_Nephelometry   0.00000000 0.000000000
#> sigma2_alpha         1.29063570 0.312335810
#> sigma2_theta         0.47671067 0.113444250
#> Lambda_LA            2.65810190 0.191815720
#> Lambda_ELISA         2.65810190 0.191815720
#> Lambda_Nephelometry  2.65810190 0.191815720
#> Theta_LA            -0.55746273 0.212932400
#> Theta_ELISA         -0.75132530 0.208099240
#> Theta_Nephelometry  -0.06043143 0.209582940
#> beta_LA              0.18948853 0.164729720
#> beta_ELISA           0.18948853 0.164729720
#> beta_Nephelometry    0.18948853 0.164729720
#> logitsens            1.00734590 0.165264820
#> logitsens            1.00734590 0.165264820
#> logitsens            1.00734590 0.165264820
#> sens                 0.73250042 0.032382577
#> sens                 0.73250042 0.032382577
#> sens                 0.73250042 0.032382577
#> 
#> Sensitivity / specificity
#> -------------------------
#>       subgroup       spec conflevel logitsens  Std_Error   CI_Lower CI_Upper
#> 1           LA 0.84615385      0.95 1.0073459 0.16526482 0.68343284 1.331259
#> 2        ELISA 0.84615385      0.95 1.0073459 0.16526482 0.68343284 1.331259
#> 3 Nephelometry 0.84615385      0.95 1.0073459 0.16526482 0.68343284 1.331259
#>         Sens SensCI_Lower SensCI_Upper
#> 1 0.73250042   0.66450444   0.79104882
#> 2 0.73250042   0.66450444   0.79104882
#> 3 0.73250042   0.66450444   0.79104882
#> 
#> Recovered bivariate parameters (Reitsma)
#> ----------------------------------------
#>               mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec    sigma_AB
#> LA           0.70184091 2.0739910    0.66138418    0.96614309 -0.15405174
#> ELISA        0.52550245 2.2871192    0.66138418    0.96614309 -0.15405174
#> Nephelometry 1.15394330 1.5275659    0.66138418    0.96614309 -0.15405174
#> 
#> Subgroups
#> ---------
#> [1] "LA"           "ELISA"        "Nephelometry"

How do I get the log likelihood, the AIC, and the BIC of a model?

logLik(ruttergatsonissubfull)
#> 'log Lik.' -376.7765 (df=11)
AIC(ruttergatsonissubfull)
#> [1] 775.5531
BIC(ruttergatsonissubfull)
#> [1] 795.9047

logLik(ruttergatsonissub)
#> 'log Lik.' -376.8612 (df=9)
AIC(ruttergatsonissub)
#> [1] 771.7223
BIC(ruttergatsonissub)
#> [1] 788.3736

logLik(ruttergatsonissubacc)
#> 'log Lik.' -377.1659 (df=7)
AIC(ruttergatsonissubacc)
#> [1] 768.3318
BIC(ruttergatsonissubacc)
#> [1] 781.2829

How do I perform a likelihood ratio test of the two models?

The tests below formally investigate whether the HSROC curve shapes and/or accuracy parameters are equal in all subgroups.

anova(ruttergatsonissub,
      ruttergatsonissubfull)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)
#> Model 1  9 -376.86                          
#> Model 2 11 -376.78       2 0.1692     0.9189
anova(ruttergatsonissubacc,
      ruttergatsonissub)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)
#> Model 1  7 -377.17                          
#> Model 2  9 -376.86       2 0.6095     0.7373
anova(ruttergatsonissubacc,
      ruttergatsonissubfull)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)
#> Model 1  7 -377.17                          
#> Model 2 11 -376.78       4 0.7787     0.9413

How do I replicate results from the Cochrane Handbook with respect to the schuetz data set?

What are the results for the full data set?

data(schuetz)
head(schuetz)
#>   test          study TP FP FN TN indirect
#> 1   CT Achenbach 2005 25  4  0 19        1
#> 2   CT   Alkadhi 2008 57 12  2 79        1
#> 3   CT  Andreini 2007 17  0  0 44        1
#> 4   CT    Bayrak 2008 64  4  0 32        1
#> 5  MRI    Bedaux 2002  7  1  0  1        1
#> 6  MRI   Bogaert 2003 12  3  3  1        1
schuetz$test   <- factor(schuetz$test,levels=c("MRI","CT"))
schuetzreitsma <- fitReitsmaSubgroup(data=schuetz,
                                     TP=TP,FP=FP,FN=FN,TN=TN,
                                     study=study,
                                     subgroup=test,
                                     variances="common")
schuetzreitsma
#> 
#> Reitsma Subgroup Model
#> ----------------------
#> 
#> Number of studies   : 108 
#> Number of subgroups : 2 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0009488887 
#> -2 log likelihood   : 951.843 ( df = 7 )
#> AIC                 : 965.843 
#> BIC                 : 984.618 
#> 
#> Use summary() for parameter estimates.
summary(schuetzreitsma)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                             Estimate  Std_Error
#> mu_A.MRI                  2.09345490 0.26392178
#> mu_B.MRI                  0.86047288 0.25722689
#> mu_A.CT                   3.48272650 0.15698592
#> mu_B.CT                   1.92906950 0.12066660
#> sigma2_A.sens.log_sigma_A 0.85004519 0.21954718
#> sigma2_B.spec.log_sigma_B 0.85268704 0.16833879
#> sigma_AB                  0.18575027 0.13420795
#> nu_A.CT                   1.38928530 0.30152230
#> nu_B.CT                   1.06858890 0.28340397
#> 
#> Sensitivity / specificity
#> -------------------------
#>          type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.MRI sens 0.89026540      0.95 0.82866250 0.93154909
#> mu_B.MRI spec 0.70275944      0.95 0.58814870 0.79651089
#> mu_A.CT  sens 0.97019227      0.95 0.95988410 0.97791261
#> mu_B.CT  spec 0.87314639      0.95 0.84456148 0.89711495
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>        Lambda      Theta        beta sigma2_alpha sigma2_theta
#> MRI 2.9548851 0.61763696 0.001551535    2.0742307   0.33280741
#> CT  5.4130029 0.77892791 0.001551535    2.0742307   0.33280741
#> 
#> Subgroups
#> ---------
#> [1] "MRI" "CT"
plot(schuetzreitsma,
     nudge_legend=-0.2,
     size="se",
     col=c("red","black"))


schuetzreitsma1 <- fitReitsmaSubgroup(data=schuetz,
                                      TP=TP,FP=FP,FN=FN,TN=TN,
                                      study=study,
                                      subgroup=test,
                                      variances="unequal")
schuetzreitsma1
#> 
#> Reitsma Subgroup Model
#> ----------------------
#> 
#> Number of studies   : 108 
#> Number of subgroups : 2 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 0.0001117501 
#> -2 log likelihood   : 943.305 ( df = 10 )
#> AIC                 : 963.305 
#> BIC                 : 990.127 
#> 
#> Use summary() for parameter estimates.
summary(schuetzreitsma1)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                 Estimate  Std_Error
#> mu_A.MRI      1.96704120 0.16243413
#> mu_B.MRI      0.84102032 0.24063701
#> sigma2_A.MRI  0.11264202 0.16001562
#> sigma2_B.MRI  0.71212692 0.34247434
#> sigma_AB.MRI -0.14620818 0.13705823
#> mu_A.CT       3.55943060 0.17395078
#> mu_B.CT       1.93164220 0.12246545
#> sigma2_A.CT   1.10692640 0.30999078
#> sigma2_B.CT   0.88004569 0.19038875
#> sigma_AB.CT   0.30489374 0.17465399
#> nu_A.CT       1.59239900 0.23799792
#> nu_B.CT       1.09062870 0.27000566
#> 
#> Sensitivity / specificity
#> -------------------------
#>          type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.MRI sens 0.87729295      0.95 0.83871205 0.90766116
#> mu_B.MRI spec 0.69868006      0.95 0.59130927 0.78795796
#> mu_A.CT  sens 0.97233226      0.95 0.96152431 0.98016685
#> mu_B.CT  spec 0.87343107      0.95 0.84443634 0.89767647
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>        Lambda      Theta        beta sigma2_alpha sigma2_theta
#> MRI 3.6494747 1.29435110  0.92202065   0.27403018   0.21471572
#> CT  5.4067108 0.65771141 -0.11468430   2.58376400   0.34104725
#> 
#> Subgroups
#> ---------
#> [1] "MRI" "CT"
plot(schuetzreitsma1,
     nudge_legend=-0.2,
     size="se",
     col=c("red","black"))


schuetzreitsma2 <- fitReitsmaSubgroup(data=schuetz,
                                      TP=TP,FP=FP,FN=FN,TN=TN,
                                      study=study,
                                      subgroup=test,
                                      sensspec_constrain="sens")
schuetzreitsma3 <- fitReitsmaSubgroup(data=schuetz,
                                      TP=TP,FP=FP,FN=FN,TN=TN,
                                      study=study,
                                      subgroup=test,
                                      sensspec_constrain="spec")
anova(schuetzreitsma,schuetzreitsma1)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)  
#> Model 1  7 -475.92                            
#> Model 2 10 -471.65       3 8.5381    0.03611 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(schuetzreitsma2,schuetzreitsma)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1  6 -485.79                              
#> Model 2  7 -475.92       1 19.745  8.848e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(schuetzreitsma3,schuetzreitsma)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1  6 -482.64                              
#> Model 2  7 -475.92       1 13.432  0.0002474 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What are the results for the direct comparisons?

schuetz2 <- subset(schuetz,indirect==0)
schuetzreitsma4 <- fitReitsmaSubgroup(data=schuetz2,
                                      TP=TP,FP=FP,FN=FN,TN=TN,
                                      study=study,
                                      subgroup=test,
                                      constrain="sigma2_A")
schuetzreitsma4
#> 
#> Reitsma Subgroup Model
#> ----------------------
#> 
#> Number of studies   : 10 
#> Number of subgroups : 2 
#> Optimizer           : Converged 
#> Hessian             : Positive definite 
#> Max |grad|          : 1.866797e-05 
#> -2 log likelihood   : 88.734 ( df = 5 )
#> AIC                 : 98.734 
#> BIC                 : 100.247 
#> 
#> Use summary() for parameter estimates.
summary(schuetzreitsma4)
#> 
#> Summary: Reitsma Subgroup Model
#> -------------------------------
#> 
#> Parameter estimates
#> -------------------
#>                                Estimate  Std_Error
#> mu_A.MRI                  1.8082893e+00 0.24124001
#> mu_B.MRI                  1.0631770e+00 0.41604406
#> mu_A.CT                   2.8134122e+00 0.34318787
#> mu_B.CT                   1.8013974e+00 0.43628527
#> sigma2_A.sens.log_sigma_A 4.9303807e-32 0.00000000
#> sigma2_B.spec.log_sigma_B 5.8153371e-01 0.43286772
#> sigma_AB                  0.0000000e+00 0.00000000
#> nu_A.CT                   1.0051219e+00 0.41949313
#> nu_B.CT                   7.3821676e-01 0.60559889
#> 
#> Sensitivity / specificity
#> -------------------------
#>          type   Estimate conflevel   CI_Lower   CI_Upper
#> mu_A.MRI sens 0.85915500      0.95 0.79174360 0.90730054
#> mu_B.MRI spec 0.74329721      0.95 0.56162156 0.86745121
#> mu_A.CT  sens 0.94339630      0.95 0.89480371 0.97028818
#> mu_B.CT  spec 0.85831895      0.95 0.72036925 0.93441053
#> 
#> Recovered HSROC parameters (Rutter-Gatsonis)
#> --------------------------------------------
#>        Lambda    Theta     beta  sigma2_alpha  sigma2_theta
#> MRI 105972100 52986049 35.77261 3.3865514e-16 8.4663784e-17
#> CT  164875820 82437911 35.77261 3.3865514e-16 8.4663784e-17
#> 
#> Subgroups
#> ---------
#> [1] "MRI" "CT"
plot(schuetzreitsma4,predlevel=0.000001,
     nudge_legend=-0.2,
     size="se",scale=0.0025,
     connectstudies = TRUE,
     col=c("red","black"))

References

Reitsma, J. B., et al. (2005). Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. Journal of Clinical Epidemiology, 58(10), 982–990.

Rutter, C. M., & Gatsonis, C. A. (2001). A hierarchical regression approach to meta-analysis of diagnostic test accuracy evaluations. Statistics in Medicine, 20(19), 2865–2884.

Harbord, R. M., Deeks, J. J., Egger, M., Whiting, P., & Sterne, J. A. C. (2007). A unification of models for meta-analysis of diagnostic accuracy studies. Biostatistics, 8(2), 239–251.

Riley, R. D., Ensor, J., Jackson, D., & Burke, D. L. (2018). Deriving percentage study weights in multi-parameter meta-analysis models. Statistical Methods in Medical Research, 27(10), 2885–2905.

Hoyer, A., Hirt, S., Kuss, O. (2018). Meta-analysis of full ROC curves using bivariate time-to-event models for interval-censored data. Research Synthesis Methods, 9(1), 62-72.

Deeks, J. J., Bossuyt, P. M., Leeflang, M. M., & Takwoingi, Y. (editors) (2023). Cochrane Handbook for Systematic Reviews of Diagnostic Test Accuracy. Version 2.0 (updated July 2023). Cochrane.