Table of contents
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 Z pred
How do we include covariates in the Reitsma model?
For sensitivity in study i y A i n A i π A i
Similarly for specificity, we have the number of non-diseased individuals testing negative: y B i n B i π B i
Now we introduce a p z i
$$ \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 z i
$$\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" )<- fitReitsmaSubgroup (data= anticcp,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= generation)#> #> Reitsma Subgroup Model #> ---------------------- #> #> Number of studies : 37 #> Number of subgroups : 2 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.0004858134 #> -2 log likelihood : 533.37 ( df = 7 ) #> AIC : 547.37 #> BIC : 558.646 #> #> Use summary() for parameter estimates. summary (reitsmasub)#> $estimates #> Estimate Std_Error #> mu_A.CCP1 -0.09653883 0.22032058 #> mu_B.CCP1 3.44671920 0.29824368 #> mu_A.CCP2 0.86603855 0.12087681 #> mu_B.CCP2 3.01649066 0.16223753 #> sigma2_A.sens 0.35982866 0.10217649 #> sigma2_B.spec 0.53990927 0.18015721 #> sigma_AB -0.19684497 0.09835341 #> nu_A.CCP2 0.96257151 0.25134200 #> nu_B.CCP2 -0.43020982 0.33771185 #> #> $sensspec #> type Estimate conflevel CI_Lower CI_Upper #> mu_A.CCP1 sens 0.4758840 0.95 0.3708997 0.5830439 #> mu_B.CCP1 spec 0.9691331 0.95 0.9459445 0.9825578 #> mu_A.CCP2 sens 0.7039207 0.95 0.6522909 0.7508130 #> mu_B.CCP2 spec 0.9533136 0.95 0.9369387 0.9655926 #> #> $RutterGatsonis_recovered #> Lambda Theta beta sigma2_alpha sigma2_theta #> CCP1 3.007377 -1.6105343 0.2028866 0.4878424 0.3188056 #> CCP2 3.684000 -0.8834971 0.2028866 0.4878424 0.3188056 #> #> $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,
<- fitReitsmaSubgroup (data= anticcp,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= generation,constrain= "sigma2_A" )summary (constrainA)$ estimates#> Estimate Std_Error #> mu_A.CCP1 -5.439405e-02 0.05498530 #> mu_B.CCP1 3.446625e+00 0.29875482 #> mu_A.CCP2 9.179794e-01 0.03139172 #> mu_B.CCP2 2.996859e+00 0.16216233 #> sigma2_A.sens 4.930381e-32 0.00000000 #> sigma2_B.spec 5.396331e-01 0.17997640 #> sigma_AB 0.000000e+00 0.00000000 #> nu_A.CCP2 9.723735e-01 0.06331527 #> nu_B.CCP2 -4.497661e-01 0.33795382
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,
<- fitReitsmaSubgroup (data= anticcp,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= generation,sensspec_constrain= "sens" )summary (constrainsens)$ estimates#> Estimate Std_Error #> mu_A.CCP1 0.65330520 0.1274076 #> mu_B.CCP1 3.08122030 0.3256178 #> mu_A.CCP2 0.65330520 0.0000000 #> mu_B.CCP2 3.11800628 0.1745128 #> sigma2_A.sens 0.54192926 0.1462735 #> sigma2_B.spec 0.57765267 0.1996838 #> sigma_AB -0.27715832 0.1398143 #> nu_A.CCP2 0.00000000 0.0000000 #> nu_B.CCP2 0.03678907 0.3857733
assumes equal (logit) sensitivities in all subgroups. Note that you can also set subgroup_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?
<- 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 we include covariates in the Rutter and Gatsonis model?
The number of diseased individuals from study i y i 1n i 1π i 1
Similarly, the number of non-diseased individuals who test positive is y i 2n i 2π i 2
Now we introduce a p z i
logit (π i j z i ⊤ θ i z i ⊤ α i x i j z i ⊤ β x i j
z i ⊤ α i z i ⊤ Λ , σ α 2 ), z i ⊤ θ i z i ⊤ Θ , σ θ 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 z i
$$\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.
data ("RF" )<- RF[RF$ method %in% c ("LA" ,"ELISA" ,"Nephelometry" ),]$ method <- factor (RF2$ method,levels= c ("LA" ,"ELISA" ,"Nephelometry" ))<- fitRutterGatsonisSubgroup (data= RF2,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= method,constrain= "shape" ) # assumes equal # shapes in subgroups
#> #> Rutter & Gatsonis Subgroup Model #> -------------------------------- #> #> Number of studies : 47 #> Number of subgroups : 3 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.0001105815 #> -2 log likelihood : 753.722 ( df = 9 ) #> AIC : 771.722 #> BIC : 788.374 #> #> Use summary() for parameter estimates. summary (ruttergatsonissub)#> $estimates #> Estimate Std. Error #> Lambda_LA 2.45630477 0.32357688 #> xi_ELISA 0.24853315 0.43954107 #> xi_Nephelometry 0.33142788 0.44260056 #> Theta_LA -0.55005213 0.21334028 #> gamma_ELISA -0.19625486 0.26096194 #> gamma_Nephelometry 0.49585945 0.26223043 #> beta_LA 0.19768573 0.16983959 #> delta_ELISA 0.00000000 0.00000000 #> delta_Nephelometry 0.00000000 0.00000000 #> sigma2_alpha 1.27791458 0.30852474 #> sigma2_theta 0.47684025 0.11344515 #> Lambda_LA 2.45630477 0.32357688 #> Lambda_ELISA 2.70483792 0.32695410 #> Lambda_Nephelometry 2.78773265 0.30577909 #> Theta_LA -0.55005213 0.21334028 #> Theta_ELISA -0.74630699 0.20991299 #> Theta_Nephelometry -0.05419268 0.21213124 #> beta_LA 0.19768573 0.16983959 #> beta_ELISA 0.19768573 0.16983959 #> beta_Nephelometry 0.19768573 0.16983959 #> logitsens 0.82616623 0.28629916 #> logitsens 1.05130869 0.28606874 #> logitsens 1.12640178 0.27689464 #> sens 0.69554369 0.06062747 #> sens 0.74102612 0.05489842 #> sens 0.75517425 0.05119397 #> #> $sensspec #> subgroup spec conflevel logitsens Std_Error CI_Lower CI_Upper #> 1 LA 0.8461538 0.95 0.8261662 0.2862992 0.2650302 1.387302 #> 2 ELISA 0.8461538 0.95 1.0513087 0.2860687 0.4906243 1.611993 #> 3 Nephelometry 0.8461538 0.95 1.1264018 0.2768946 0.5836983 1.669105 #> Sens SensCI_Lower SensCI_Upper #> 1 0.6955437 0.5658724 0.8001612 #> 2 0.7410261 0.6202535 0.8336879 #> 3 0.7551742 0.6419179 0.8414565 #> #> $Reitsma_recovered #> mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec sigma_AB #> LA 0.6142809 1.962947 0.6534814 0.9703778 -0.1573616 #> ELISA 0.5490678 2.316769 0.6534814 0.9703778 -0.1573616 #> Nephelometry 1.2135903 1.598502 0.6534814 0.9703778 -0.1573616 #> #> $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.,
constrain="sigma2_alpha" sets σ α 2 constrain="sigma2_theta" sets σ θ 2 constrain="accuracy" assumes equal accuracy parameters across subgroups, i.e. ξ 2 = ξ 3 = … = 0constrain="threshold" assumes equal threshold parameters across subgroups, i.e. γ 2 = γ 3 = … = 0constrain="shape" assumes equal shape parameters across subgroups, i.e. δ 2 = δ 3 = … = 0constrain="shape_zero" fixes all shape parameters at zero.
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 plot() nor forest() methods for fitRutterGatsonisReg().
Let’s reproduce the subgroup-analysis from before.
# Specify design matrix Z <- model.matrix (~ method,data= RF2)# For study level-covariates, we need to two identical consecutive # rows per study (diseased and non-diseased). <- Z[rep (seq_len (nrow (Z)), each = 2 ), , drop = FALSE ]# Specify prediction matrix Z_pred <- matrix (c (1 ,0 ,0 ,1 ,1 ,0 ,1 ,0 ,1 ),ncol= 3 ,nrow= 3 ,byrow= T)<- list (shape_coef= factor (c (1 , rep (NA , ncol (Z2) - 1 ))))<- fitRutterGatsonisReg (data= RF2,TP= TP,FP= FP,FN= FN,TN= TN,study= study,Z= Z2,Z_pred= Z_pred,map= constrain)#> #> Rutter & Gatsonis Regression Model #> ---------------------------------- #> #> Number of studies : 47 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.0005885859 #> -2 log likelihood : 753.722 ( df = 9 ) #> AIC : 771.722 #> BIC : 788.374 #> #> Use summary() for parameter estimates. summary (ruttergatsonisreg)#> $estimates #> Estimate Std. Error #> accuracy_coef 2.45631156 0.32357673 #> accuracy_coef 0.24852208 0.43954066 #> accuracy_coef 0.33141584 0.44260014 #> threshold_coef -0.55005084 0.21334149 #> threshold_coef -0.19625009 0.26096393 #> threshold_coef 0.49585695 0.26223232 #> shape_coef 0.19768191 0.16983950 #> shape_coef 0.00000000 0.00000000 #> shape_coef 0.00000000 0.00000000 #> sigma2_alpha 1.27791200 0.30852400 #> sigma2_theta 0.47684805 0.11344753 #> Lambda_Pred 2.45631156 0.32357673 #> Lambda_Pred 2.70483365 0.32695391 #> Lambda_Pred 2.78772741 0.30577879 #> Theta_Pred -0.55005084 0.21334149 #> Theta_Pred -0.74630093 0.20991415 #> Theta_Pred -0.05419388 0.21213237 #> beta_Pred 0.19768191 0.16983950 #> beta_Pred 0.19768191 0.16983950 #> beta_Pred 0.19768191 0.16983950 #> logitsens 0.82617129 0.28629952 #> logitsens 1.05130416 0.28606887 #> logitsens 1.12639652 0.27689490 #> sens 0.69554476 0.06062743 #> sens 0.74102525 0.05489857 #> sens 0.75517328 0.05119416 #> #> $sensspec #> spec conflevel logitsens Std_Error CI_Lower CI_Upper Sens #> 1 0.8461538 0.95 0.8261713 0.2862995 0.2650345 1.387308 0.6955448 #> 2 0.8461538 0.95 1.0513042 0.2860689 0.4906195 1.611989 0.7410253 #> 3 0.8461538 0.95 1.1263965 0.2768949 0.5836925 1.669101 0.7551733 #> SensCI_Lower SensCI_Upper #> 1 0.5658735 0.8001621 #> 2 0.6202524 0.8336873 #> 3 0.6419166 0.8414559
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.
<- fitRutterGatsonisSubgroup (data= RF2,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= method,constrain= NULL )#> #> Rutter & Gatsonis Subgroup Model #> -------------------------------- #> #> Number of studies : 47 #> Number of subgroups : 3 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.000337596 #> -2 log likelihood : 753.553 ( df = 11 ) #> AIC : 775.553 #> BIC : 795.905 #> #> Use summary() for parameter estimates. summary (ruttergatsonissubfull)#> $estimates #> Estimate Std. Error #> Lambda_LA 2.4243669 0.33049225 #> xi_ELISA 0.2974583 0.51315392 #> xi_Nephelometry 0.3716007 0.45086693 #> Theta_LA -0.5029487 0.24451662 #> gamma_ELISA -0.2578242 0.37024676 #> gamma_Nephelometry 0.3970656 0.35944299 #> beta_LA 0.2777491 0.26642343 #> delta_ELISA -0.1028887 0.42661221 #> delta_Nephelometry -0.1603831 0.39753788 #> sigma2_alpha 1.2730677 0.30798947 #> sigma2_theta 0.4762396 0.11343569 #> Lambda_LA 2.4243669 0.33049225 #> Lambda_ELISA 2.7218252 0.39299597 #> Lambda_Nephelometry 2.7959676 0.30708175 #> Theta_LA -0.5029487 0.24451662 #> Theta_ELISA -0.7607729 0.27817185 #> Theta_Nephelometry -0.1058831 0.26322843 #> beta_LA 0.2777491 0.26642343 #> beta_ELISA 0.1748604 0.33321541 #> beta_Nephelometry 0.1173661 0.29500322 #> logitsens 0.8186925 0.27519154 #> logitsens 1.0626991 0.32405869 #> logitsens 1.1206495 0.28834985 #> sens 0.6939587 0.05844518 #> sens 0.7432060 0.06184687 #> sens 0.7541092 0.05346829 #> #> $sensspec #> subgroup spec conflevel logitsens Std_Error CI_Lower CI_Upper #> 1 LA 0.8461538 0.95 0.8186925 0.2751915 0.2793270 1.358058 #> 2 ELISA 0.8461538 0.95 1.0626991 0.3240587 0.4275558 1.697843 #> 3 Nephelometry 0.8461538 0.95 1.1206495 0.2883499 0.5554942 1.685805 #> Sens SensCI_Lower SensCI_Upper #> 1 0.6939587 0.5693812 0.7954439 #> 2 0.7432060 0.6052899 0.8452527 #> 3 0.7541092 0.6354094 0.8436717 #> #> $Reitsma_recovered #> mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec sigma_AB #> LA 0.6172734 1.970652 0.6018282 1.0488717 -0.1579727 #> ELISA 0.5498977 2.315536 0.6670472 0.9463208 -0.1579727 #> Nephelometry 1.2184583 1.594759 0.7065226 0.8934472 -0.1579727 #> #> $subgroups #> [1] "LA" "ELISA" "Nephelometry" <- fitRutterGatsonisSubgroup (data= RF2,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= method,constrain= c ("shape" ,"accuracy" ))#> #> Rutter & Gatsonis Subgroup Model #> -------------------------------- #> #> Number of studies : 47 #> Number of subgroups : 3 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.00030144 #> -2 log likelihood : 754.332 ( df = 7 ) #> AIC : 768.332 #> BIC : 781.283 #> #> Use summary() for parameter estimates. summary (ruttergatsonissubacc)#> $estimates #> Estimate Std. Error #> Lambda_LA 2.65809857 0.19181608 #> xi_ELISA 0.00000000 0.00000000 #> xi_Nephelometry 0.00000000 0.00000000 #> Theta_LA -0.55745914 0.21293277 #> gamma_ELISA -0.19385254 0.26093654 #> gamma_Nephelometry 0.49702816 0.26213094 #> beta_LA 0.18949665 0.16472998 #> delta_ELISA 0.00000000 0.00000000 #> delta_Nephelometry 0.00000000 0.00000000 #> sigma2_alpha 1.29064650 0.31233907 #> sigma2_theta 0.47671277 0.11344490 #> Lambda_LA 2.65809857 0.19181608 #> Lambda_ELISA 2.65809857 0.19181608 #> Lambda_Nephelometry 2.65809857 0.19181608 #> Theta_LA -0.55745914 0.21293277 #> Theta_ELISA -0.75131168 0.20809954 #> Theta_Nephelometry -0.06043098 0.20958332 #> beta_LA 0.18949665 0.16472998 #> beta_ELISA 0.18949665 0.16472998 #> beta_Nephelometry 0.18949665 0.16472998 #> logitsens 1.00734453 0.16526474 #> logitsens 1.00734453 0.16526474 #> logitsens 1.00734453 0.16526474 #> sens 0.73250015 0.03238258 #> sens 0.73250015 0.03238258 #> sens 0.73250015 0.03238258 #> #> $sensspec #> subgroup spec conflevel logitsens Std_Error CI_Lower CI_Upper #> 1 LA 0.8461538 0.95 1.007345 0.1652647 0.6834316 1.331257 #> 2 ELISA 0.8461538 0.95 1.007345 0.1652647 0.6834316 1.331257 #> 3 Nephelometry 0.8461538 0.95 1.007345 0.1652647 0.6834316 1.331257 #> Sens SensCI_Lower SensCI_Upper #> 1 0.7325001 0.6645042 0.7910486 #> 2 0.7325001 0.6645042 0.7910486 #> 3 0.7325001 0.6645042 0.7910486 #> #> $Reitsma_recovered #> mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec sigma_AB #> LA 0.7018398 2.073994 0.6613828 0.9661567 -0.1540511 #> ELISA 0.5255112 2.287112 0.6613828 0.9661567 -0.1540511 #> Nephelometry 1.1539375 1.527570 0.6613828 0.9661567 -0.1540511 #> #> $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
The tests below formally investigate whether the HSROC curve shapes and/or accuracy parameters are equal in all subgroups.
anova (ruttergatsonissub,#> Df logLik Df.diff Chisq Pr(>Chisq) #> Model 1 9 -376.86 #> Model 2 11 -376.78 2 0.1692 0.9189 anova (ruttergatsonissubacc,#> Df logLik Df.diff Chisq Pr(>Chisq) #> Model 1 7 -377.17 #> Model 2 9 -376.86 2 0.6095 0.7373 anova (ruttergatsonissubacc,#> 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 $ test <- factor (schuetz$ test,levels= c ("MRI" ,"CT" ))<- fitReitsmaSubgroup (data= schuetz,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= test,variances= "common" )#> #> Reitsma Subgroup Model #> ---------------------- #> #> Number of studies : 108 #> Number of subgroups : 2 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.0006484896 #> -2 log likelihood : 951.843 ( df = 7 ) #> AIC : 965.843 #> BIC : 984.618 #> #> Use summary() for parameter estimates. summary (schuetzreitsma)#> $estimates #> Estimate Std_Error #> mu_A.MRI 2.0934567 0.2639212 #> mu_B.MRI 0.8604695 0.2572264 #> mu_A.CT 3.4827280 0.1569859 #> mu_B.CT 1.9290688 0.1206664 #> sigma2_A.sens 0.8500392 0.2195458 #> sigma2_B.spec 0.8526826 0.1683378 #> sigma_AB 0.1857463 0.1342067 #> nu_A.CT 1.3893006 0.3015220 #> nu_B.CT 1.0686026 0.2834029 #> #> $sensspec #> type Estimate conflevel CI_Lower CI_Upper #> mu_A.MRI sens 0.8902656 0.95 0.8286629 0.9315491 #> mu_B.MRI spec 0.7027587 0.95 0.5881481 0.7965102 #> mu_A.CT sens 0.9701923 0.95 0.9598842 0.9779126 #> mu_B.CT spec 0.8731463 0.95 0.8445615 0.8971148 #> #> $RutterGatsonis_recovered #> Lambda Theta beta sigma2_alpha sigma2_theta #> MRI 2.954884 0.6176402 0.001552424 2.074212 0.3328068 #> CT 5.413004 0.7789302 0.001552424 2.074212 0.3328068 #> #> $subgroups #> [1] "MRI" "CT" plot (schuetzreitsma,nudge_legend= - 0.2 ,size= "se" ,col= c ("red" ,"black" ))<- fitReitsmaSubgroup (data= schuetz,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= test,variances= "unequal" )#> #> Reitsma Subgroup Model #> ---------------------- #> #> Number of studies : 108 #> Number of subgroups : 2 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 0.0002796148 #> -2 log likelihood : 943.305 ( df = 10 ) #> AIC : 963.305 #> BIC : 990.127 #> #> Use summary() for parameter estimates. summary (schuetzreitsma1)#> $estimates #> Estimate Std_Error #> mu_A.MRI 1.9670363 0.1624331 #> mu_B.MRI 0.8410193 0.2406372 #> sigma2_A.MRI 0.1126405 0.1600126 #> sigma2_B.MRI 0.7121292 0.3424753 #> sigma_AB.MRI -0.1462045 0.1370573 #> mu_A.CT 3.5594295 0.1739508 #> mu_B.CT 1.9316450 0.1224654 #> sigma2_A.CT 1.1069322 0.3099921 #> sigma2_B.CT 0.8800412 0.1903877 #> sigma_AB.CT 0.3048890 0.1746537 #> nu_A.CT 1.5923912 0.2379991 #> nu_B.CT 1.0906240 0.2700068 #> #> $sensspec #> type Estimate conflevel CI_Lower CI_Upper #> mu_A.MRI sens 0.8772924 0.95 0.8387117 0.9076606 #> mu_B.MRI spec 0.6986799 0.95 0.5913089 0.7879579 #> mu_A.CT sens 0.9723322 0.95 0.9615243 0.9801668 #> mu_B.CT spec 0.8734314 0.95 0.8444367 0.8976767 #> #> $RutterGatsonis_recovered #> Lambda Theta beta sigma2_alpha sigma2_theta #> MRI 3.649477 1.2943550 0.9220289 0.2740348 0.2147132 #> CT 5.406709 0.6577024 -0.1146895 2.5837546 0.3410496 #> #> $subgroups #> [1] "MRI" "CT" plot (schuetzreitsma1,nudge_legend= - 0.2 ,size= "se" ,col= c ("red" ,"black" ))<- fitReitsmaSubgroup (data= schuetz,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= test,sensspec_constrain= "sens" )<- 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?
<- subset (schuetz,indirect== 0 )<- fitReitsmaSubgroup (data= schuetz2,TP= TP,FP= FP,FN= FN,TN= TN,study= study,subgroup= test,constrain= "sigma2_A" )#> #> Reitsma Subgroup Model #> ---------------------- #> #> Number of studies : 10 #> Number of subgroups : 2 #> Optimizer : Converged #> Hessian : Positive definite #> Max |grad| : 1.86672e-05 #> -2 log likelihood : 88.734 ( df = 5 ) #> AIC : 98.734 #> BIC : 100.247 #> #> Use summary() for parameter estimates. round (summary (schuetzreitsma4)$ estimates,5 )#> Estimate Std_Error #> mu_A.MRI 1.80829 0.24124 #> mu_B.MRI 1.06318 0.41604 #> mu_A.CT 2.81341 0.34319 #> mu_B.CT 1.80140 0.43629 #> sigma2_A.sens 0.00000 0.00000 #> sigma2_B.spec 0.58153 0.43287 #> sigma_AB 0.00000 0.00000 #> nu_A.CT 1.00512 0.41949 #> nu_B.CT 0.73822 0.60560 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.
