Latent Class Meta-Analysis

Reitsma and Rutter-Gatsonis Models with an Imperfect Reference Standard

What are latent class models?

In conventional meta-analysis of diagnostic test accuracy studies, the reference standard is assumed to classify disease status without error. In practice, however, the reference test may itself be imperfect. Latent class models relax the gold-standard assumption by treating the true disease status as an unobserved (latent) variable and jointly modelling the results of the index and reference tests.

Which latent class models does dtametaTMB support?

The latent class models in dtametaTMB are based on the reduced latent class formulations proposed by Liu et al. (2015).

Liu, Y., Chen, Y., & Chu, H. (2015). A unification of models for meta-analysis of diagnostic accuracy studies without a gold standard. Biometrics, 71(2), 538–547.

The latent class implementations were validated by reproducing the SAS PROC NLMIXED examples provided by Liu et al. (2015).

What is conditional independence?

The latent class models in dtametaTMB assume conditional independence between the index test and the reference test given the true disease status. Under this assumption, any association between the index and reference tests is fully explained by the latent disease status. Consequently, the latent-class probabilities can be expressed as functions of disease prevalence and the sensitivities and specificities of the two tests without introducing additional covariance parameters.

To account for residual dependence between the index and reference tests, more complex latent class models do include covariance terms among diseased and/or non-diseased individuals. Although potentially more flexible, these models are also more difficult to identify and estimate reliably, particularly in sparse data.

How do latent class models work?

Consider study i. Let

Let

ni = y11i + y10i + y01i + y00i

denote the total sample size in study i.

Furthermore, let

p11i,  p10i,  p01i,  p00i

denote the corresponding latent-class cell probabilities.

Conditional on the model parameters, the observed counts are assumed to follow a multinomial distribution,

(Y11i, Y10i, Y01i, Y00i) ∼ Multinomial(ni, (p11i, p10i, p01i, p00i)).

Ignoring the multinomial coefficient, which does not depend on the model parameters, the study-specific likelihood is proportional to

Li ∝ p11iy11ip10iy10ip01iy01ip00iy00i.

Equivalently, the study-specific log-likelihood can be written as

i = y11ilog (p11i) + y10ilog (p10i) + y01ilog (p01i) + y00ilog (p00i),

which corresponds directly to the likelihood contribution implemented in the TMB template.

How do we model disease prevalence?

Let πi denote the prevalence of disease in study i. We model the study-specific prevalence by

logit(πi) = μprev + uprev, i,

where

uprev, i ∼ 𝒩(0, σprev2).

The parameter μprev represents the average logit prevalence across studies, while σprev2 describes between-study heterogeneity in disease prevalence.

How do we model the index test? (Reitsma version)

Let Sei and Spi denote the sensitivity and specificity of the index test in study i.

The study-specific logit sensitivities and specificities are modeled as

logit(Sei) = μA, index + uAi,

and

logit(Spi) = μB, index + uBi,

where

$$ \begin{pmatrix} u_{Ai}\\ u_{Bi} \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} 0\\ 0 \end{pmatrix}, \begin{pmatrix} \sigma^2_{A,\mathrm{index}} & \sigma_{AB,\mathrm{index}}\\[4pt] \sigma_{AB,\mathrm{index}} & \sigma^2_{B,\mathrm{index}} \end{pmatrix} \right). $$

The parameters

μA, index  and  μB, index

represent the average logit sensitivity and specificity of the index test, respectively.

The variance components

σA, index2  and  σB, index2

describe heterogeneity in sensitivity and specificity, whereas

σAB, index

captures the association between them.

How do we model the index test? (Rutter and Gatsonis version)

Let xij indicate disease status, where

$$ x_{ij} = \begin{cases} +0.5, & \text{diseased},\\ -0.5, & \text{non-diseased}. \end{cases} $$

The Rutter and Gatsonis model is given by

ηij = [Θ + θi + (Λ + αi)xij]exp (−βxij),

where

θi ∼ 𝒩(0, σθ2),  and  αi ∼ 𝒩(0, σα2). The parameter Λ controls diagnostic accuracy, whereas Θ controls the diagnostic threshold. The shape parameter β determines the asymmetry of the summary ROC curve.

Using the Rutter and Gatsonis parameterization, the study-specific sensitivity and specificity of the index test are

$$ \text{logit}(Se_{i}) = \left( \Theta + \theta_i + \frac{\Lambda+\alpha_i}{2} \right) \exp\!\left(-\frac{\beta}{2}\right) $$

and

$$ \text{logit}(Sp_{i}) = - \left( \Theta + \theta_i - \frac{\Lambda+\alpha_i}{2} \right) \exp\!\left(\frac{\beta}{2}\right) . $$

How do we model the reference test?

The reduced latent class models of Liu et al. (2015) assume that the reference test has common sensitivity and specificity across studies.

SeR = logit−1(μA, ref),

and

SpR = logit−1(μB, ref).

The parameters

μA, ref  and  μB, ref

represent the logit sensitivity and specificity of the reference test, respectively.

How do we model the latent-class probabilities?

Under conditional independence the latent class probabilities for study i are then given by

p11i = πiSeiSeR + (1 − πi)(1 − Spi)(1 − SpR),

p10i = πiSei(1 − SeR) + (1 − πi)(1 − Spi)SpR,

p01i = πi(1 − Sei)SeR + (1 − πi)Spi(1 − SpR),

p00i = πi(1 − Sei)(1 − SeR) + (1 − πi)SpiSpR.

These probabilities correspond to the four observable cells

(y11i, y10i, y01i, y00i),

while the true disease status remains unobserved.

Model parameters (Reitsma version)

The reduced Reitsma latent class model therefore estimates

μprev,  μA, index,  μB, index,  μA, ref,  μB, ref,

together with

σprev2,  σA, index2,  σB, index2,  σAB, index.

Model parameters (Rutter and Gatsonis version)

The reduced HSROC latent class model estimates

μprev,  Λ,  Θ,  β,  μA, ref,  μB, ref,

together with

σprev2,  σα2,  σθ2.

For some datasets, the threshold heterogeneity parameter σθ2 may be estimated close to zero, suggesting that a reduced model with fixed threshold heterogeneity may provide an adequate description of the data. This phenomenon is discussed by Liu et al. (2015) and can also arise in practical applications when most between-study variation is attributable to changes in diagnostic accuracy rather than diagnostic threshold.

How does the model provide study-specific sensitivities/specificities?

Because the true disease status is unobserved, study-specific sensitivities and specificities are not directly available from the observed data. Instead, they are obtained as empirical Bayes estimates derived from the fitted latent class model. Conditional empirical Bayes variances are used to quantify uncertainty in study-level summaries and graphical displays.

How do we obtain starting values?

Starting values are obtained from simpler nested models. For example, the latent class models are initialized using parameter estimates from the corresponding non-latent-class Reitsma or HSROC models. Subgroup models use the fitted non-subgroup latent class model to obtain initial values for common parameters, while subgroup contrasts are initialized at zero.

How do I fit the Reitsma latent class model?

data("pap")
papLCA <- fitReitsmaLCA(data=pap,
                        y11=y11,
                        y10=y10,
                        y01=y01,
                        y00=y00,
                        study=id)
papLCA
#> 
#> Reitsma LCA Model
#> -----------------
#> 
#> Number of studies : 59 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.002932841 
#> -2 log likelihood : 39752.25 ( df = 9 )
#> AIC               : 39770.25 
#> BIC               : 39788.94 
#> 
#> mu_prev           : 0.556 
#> mu_A.index (sens) : 0.64 
#> mu_B.index (spec) : 1.62 
#> 
#> Use summary() for parameter estimates.
summary(papLCA)
#> $estimates
#>                  Estimate Std. Error
#> mu_prev         0.5560178  0.2081842
#> mu_A.index      0.6402414  0.1843872
#> mu_B.index      1.6195488  0.2294594
#> sigma2_prev     2.1528711  0.4809724
#> sigma2_A.index  1.6679041  0.3510982
#> sigma2_B.index  1.6093096  0.4167185
#> sigma_AB.index -0.8845144  0.3075786
#> rho_AB.index   -0.5398829  0.1357354
#> mu_A.ref        2.2247864  0.1471293
#> mu_B.ref        4.4675482  1.2790747
#> 
#> $sensspec
#>            type  Estimate conflevel  CI_Lower  CI_Upper
#> mu_prev    Prev 0.6355306      0.95 0.5369287 0.7239322
#> mu_A.index Sens 0.6548080      0.95 0.5692640 0.7313797
#> mu_B.index Spec 0.8347329      0.95 0.7631119 0.8878814
#> mu_A.ref   Sens 0.9024534      0.95 0.8739581 0.9250589
#> mu_B.ref   Spec 0.9886548      0.95 0.8765987 0.9990654
#> 
#> $RutterGatsonis_recovered
#>                        Lambda      Theta        beta sigma2_alpha sigma2_theta
#> Estimate (recovered) 2.268636 -0.4997754 -0.01788129     1.507661      1.26143

How do I get a summary plot?

plot(papLCA)

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(papLCA)

How do I fit the Rutter and Gatsonis latent class model?

papLCA2 <- fitRutterGatsonisLCA(data=pap,
                                y11=y11,
                                y10=y10,
                                y01=y01,
                                y00=y00,
                                study=id)
papLCA2
#> 
#> Rutter & Gatsonis LCA Model
#> ---------------------------
#> 
#> Number of studies : 59 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.006392151 
#> -2 log likelihood : 39752.25 ( df = 9 )
#> AIC               : 39770.25 
#> BIC               : 39788.94 
#> 
#> mu_prev           : 0.556 
#> Lambda            : 2.269 
#> Theta             : -0.5 
#> beta              : -0.018 
#> 
#> Use summary() for parameter estimates.
summary(papLCA2)
#> $estimates
#>                 Estimate Std. Error
#> mu_prev       0.55603885 0.20818876
#> Lambda        2.26862858 0.24787026
#> Theta        -0.49976968 0.17614575
#> beta         -0.01790735 0.15271733
#> sigma2_prev   2.15293499 0.48100048
#> sigma2_alpha  1.50784637 0.45147296
#> sigma2_theta  1.26121398 0.27922251
#> mu_A.ref      2.22479234 0.14713722
#> mu_B.ref      4.46769888 1.27954190
#> logitsens     0.87768953 0.24730606
#> sens          0.70634321 0.05129684
#> 
#> $sensspec
#>   spec conflevel logitsens Std_Error  CI_Lower CI_Upper      Sens SensCI_Lower
#> 1  0.8      0.95 0.8776895 0.2473061 0.3929786 1.362401 0.7063432    0.5969995
#>   SensCI_Upper
#> 1    0.7961496
#> 
#> $Reitsma_recovered
#>                      mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec   sigma_AB
#> Estimate (recovered) 0.6402516  1.619518      1.667775      1.609101 -0.8842524

How do I get a summary plot?

plot(papLCA2,specrange=c(0.3,0.995))

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(papLCA2)

How do I perform subgroup analyses?

In the subgroup models, the reference test’s logit mean sensitivity and specificity are assumed common across all subgroups, whereas the logit mean prevalence and its between-study random effect variance can differ across subgroups g = 1, …, G, i.e.,

μprev, g,  σprev, g2.

However, by default σprev, 12 = … = σprev, 22 = … = σprev, G2 are assumed to be common. Separate prevalence heterogeneity parameters are estimated only when prev_variances="unequal".

Reitsma subgroup latent class model

For the index test, the Reitsma subgroup latent class model allows for different logit mean sensitivities and specificities across subgroups, i.e.,

μA, index, g,  μB, index, g.

By default (variances="common"), a common between-study variance-covariance structure is assumed across all subgroups, i.e.,

Σ1 = Σ2 = … = ΣG. Alternatively, subgroup-specific variance-covariance matrices can be estimated by specifying variances="unequal". In this case, each subgroup has its own between-study variance-covariance matrix

$$\Sigma_{g} = \begin{pmatrix} \sigma^2_{A,\text{index},g} & \sigma_{AB,\text{index},g} \\ \sigma_{AB,\text{index},g} & \sigma^2_{B,\text{index},g} \end{pmatrix}.$$ For convenience, subgroup differences relative to the reference subgroup are also reported on the logit scale through the parameters νprev, νA, index, νB, index.

Rutter and Gatsonis subgroup latent class model

For the index test, the Rutter and Gatsonis subgroup latent class model allows for different accuracy, threshold and shape effects per subgroups, i.e.,

Λg,  Θg,  βg, while the random effects variances σα2 and σθ2 are assumed common across subgroups.

For convenience, subgroup differences relative to the reference subgroup are also reported through the parameters νg,  ξg,  γg,  δg. ##

How do I fit the Reitsma subgroup latent class model?

pap$type <- factor(pap$type,levels=c("SC","FU"))
papLCA3 <- fitReitsmaSubgroupLCA(data=pap,
                                 y11=y11,
                                 y10=y10,
                                 y01=y01,
                                 y00=y00,
                                 subgroup=type,
                                 study=id)
papLCA3
#> 
#> Reitsma Subgroup LCA Model
#> --------------------------
#> 
#> Number of studies : 59 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.004538098 
#> -2 log likelihood : 39749.29 ( df = 12 )
#> AIC               : 39773.29 
#> BIC               : 39798.22 
#> 
#> Use summary() for parameter estimates.
summary(papLCA3)
#> $estimates
#>                      Estimate Std. Error
#> mu_prev.SC         0.36347507  0.2909301
#> mu_prev.FU         0.72398280  0.2754397
#> mu_A.index.SC      0.39605108  0.2638889
#> mu_A.index.FU      0.82651596  0.2407450
#> mu_B.index.SC      1.56097544  0.2989015
#> mu_B.index.FU      1.62872239  0.2948260
#> sigma2_prev.SC     2.08585705  0.4690105
#> sigma2_prev.FU     2.08585705  0.4690105
#> sigma2_A.index.SC  1.61280022  0.3391763
#> sigma2_A.index.FU  1.61280022  0.3391763
#> sigma2_B.index.SC  1.58707540  0.4074292
#> sigma2_B.index.FU  1.58707540  0.4074292
#> sigma_AB.index.SC -0.89163884  0.2992091
#> sigma_AB.index.FU -0.89163884  0.2992091
#> rho_AB.index.SC   -0.55731395  0.1322830
#> rho_AB.index.FU   -0.55731395  0.1322830
#> mu_A.ref           2.24793832  0.1511945
#> mu_B.ref           4.78612624  1.8863456
#> nu_prev.FU         0.36050773  0.3895942
#> nu_A.index.FU      0.43046488  0.3485269
#> nu_B.index.FU      0.06774695  0.3852836
#> 
#> $sensspec
#>               type  Estimate conflevel  CI_Lower  CI_Upper
#> mu_prev.SC    Prev 0.5898814      0.95 0.4484990 0.7178228
#> mu_prev.FU    Prev 0.6734835      0.95 0.5459031 0.7796855
#> mu_A.index.SC Sens 0.5977385      0.95 0.4697466 0.7136676
#> mu_A.index.FU Sens 0.6956177      0.95 0.5877482 0.7855601
#> mu_B.index.SC Spec 0.8264933      0.95 0.7261427 0.8953705
#> mu_B.index.FU Spec 0.8359945      0.95 0.7409427 0.9008380
#> mu_A.ref      Sens 0.9044726      0.95 0.8756213 0.9271879
#> mu_B.ref      Spec 0.9917243      0.95 0.7481852 0.9997931
#> 
#> $RutterGatsonis_recovered
#>      Lambda      Theta         beta sigma2_alpha sigma2_theta
#> SC 1.961725 -0.5864003 -0.008039493     1.416495     1.245762
#> FU 2.458483 -0.4060412 -0.008039493     1.416495     1.245762

How do I get a summary plot?

plot(papLCA3,
     nudge_legend = -0.2)

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(papLCA3)

How do I fit the Rutter and Gatsonis latent class model?

pap$type <- factor(pap$type,levels=c("SC","FU"))
papLCA4 <- fitRutterGatsonisSubgroupLCA(data=pap,
                                        y11=y11,
                                        y10=y10,
                                        y01=y01,
                                        y00=y00,
                                        subgroup=type,
                                        constrain="shape",
                                        study=id)
papLCA4
#> 
#> Rutter & Gatsonis Subgroup LCA Model
#> ------------------------------------
#> 
#> Number of studies : 59 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.0048649 
#> -2 log likelihood : 39749.29 ( df = 12 )
#> AIC               : 39773.29 
#> BIC               : 39798.22 
#> 
#> Use summary() for parameter estimates.
summary(papLCA4)
#> $estimates
#>                    Estimate Std. Error
#> mu_prev.SC      0.363495687 0.29094405
#> nu_prev.FU      0.360583606 0.38960774
#> Lambda_SC       1.961742518 0.32030038
#> xi_FU           0.496626279 0.38711138
#> Theta_SC       -0.586492738 0.23519964
#> gamma_FU        0.180728522 0.31327319
#> beta_SC        -0.007932961 0.15042333
#> delta_FU        0.000000000 0.00000000
#> sigma2_prev.SC  2.086003195 0.46911055
#> sigma2_prev.FU  2.086003195 0.46911055
#> sigma2_alpha    1.416302259 0.42901818
#> sigma2_theta    1.245762077 0.27344454
#> Lambda_SC       1.961742518 0.32030038
#> Lambda_FU       2.458368797 0.30223518
#> Theta_SC       -0.586492738 0.23519964
#> Theta_FU       -0.405764216 0.23297316
#> beta_SC        -0.007932961 0.15042333
#> beta_FU        -0.007932961 0.15042333
#> mu_prev.SC      0.363495687 0.29094405
#> mu_prev.FU      0.724079294 0.27545393
#> mu_A.ref        2.248031811 0.15124425
#> mu_B.ref        4.787467501 1.89052978
#> logitsens       0.572203667 0.31884101
#> logitsens       1.070803716 0.30173486
#> sens            0.639271504 0.07352584
#> sens            0.744749730 0.05735906
#> 
#> $sensspec
#>   subgroup spec conflevel logitsens Std_Error    CI_Lower CI_Upper      Sens
#> 1       SC  0.8      0.95 0.5722037 0.3188410 -0.05271323 1.197121 0.6392715
#> 2       FU  0.8      0.95 1.0708037 0.3017349  0.47941427 1.662193 0.7447497
#>   SensCI_Lower SensCI_Upper
#> 1    0.4868247    0.7680122
#> 2    0.6176096    0.8405322
#> 
#> $Reitsma_recovered
#>    mu_A.sens mu_B.spec sigma2_A.sens sigma2_B.spec   sigma_AB
#> SC 0.3959459  1.561159       1.61258      1.587196 -0.8916865
#> FU 0.8266927  1.628476       1.61258      1.587196 -0.8916865

How do I get a summary plot?

plot(papLCA4,
     specrange=c(0.3,0.995),
     nudge_legend = -0.2)

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(papLCA4)

How do I apply these methods to the schuetz data set?

data("schuetz")
schuetzLCA <- fitReitsmaLCA(data=schuetz,
                            y11=TP,
                            y10=FP,
                            y01=FN,
                            y00=TN,
                            study=study)
schuetzLCA
#> 
#> Reitsma LCA Model
#> -----------------
#> 
#> Number of studies : 108 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 3.950918e-05 
#> -2 log likelihood : 16121.38 ( df = 9 )
#> AIC               : 16139.38 
#> BIC               : 16163.52 
#> 
#> mu_prev           : 0.302 
#> mu_A.index (sens) : 3.269 
#> mu_B.index (spec) : 2.161 
#> 
#> Use summary() for parameter estimates.
summary(schuetzLCA)
#> $estimates
#>                  Estimate   Std. Error
#> mu_prev         0.3024396   0.08705366
#> mu_A.index      3.2688653   0.15041192
#> mu_B.index      2.1614724   0.24027626
#> sigma2_prev     0.6256940   0.10189887
#> sigma2_A.index  1.1494500   0.28118311
#> sigma2_B.index  1.4160394   0.36985002
#> sigma_AB.index  0.5530858   0.21827396
#> rho_AB.index    0.4335207   0.13074933
#> mu_A.ref        3.5220777   0.39936864
#> mu_B.ref       18.3550248 924.61451052
#> 
#> $sensspec
#>            type  Estimate conflevel  CI_Lower  CI_Upper
#> mu_prev    Prev 0.5750388      0.95 0.5329068 0.6161082
#> mu_A.index Sens 0.9633451      0.95 0.9513885 0.9724460
#> mu_B.index Spec 0.8967360      0.95 0.8442951 0.9329182
#> mu_A.ref   Sens 0.9713095      0.95 0.9393079 0.9866767
#> mu_B.ref   Spec 1.0000000      0.95 0.0000000 1.0000000
#> 
#> $RutterGatsonis_recovered
#>                        Lambda     Theta      beta sigma2_alpha sigma2_theta
#> Estimate (recovered) 5.495493 0.6960962 0.1042901     3.657772    0.3613573

schuetzLCA2 <- fitReitsmaSubgroupLCA(data=schuetz,
                                     y11=TP,
                                     y10=FP,
                                     y01=FN,
                                     y00=TN,
                                     subgroup=test,
                                     study=study)
schuetzLCA2
#> 
#> Reitsma Subgroup LCA Model
#> --------------------------
#> 
#> Number of studies : 108 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.007977732 
#> -2 log likelihood : 16092.65 ( df = 12 )
#> AIC               : 16116.65 
#> BIC               : 16148.83 
#> 
#> Use summary() for parameter estimates.
summary(schuetzLCA2)
#> $estimates
#>                      Estimate   Std. Error
#> mu_prev.CT          0.2798869 9.437801e-02
#> mu_prev.MRI         0.4288622 2.043445e-01
#> mu_A.index.CT       3.4806545 1.562806e-01
#> mu_A.index.MRI      2.0912725 2.629100e-01
#> mu_B.index.CT       2.3839558 2.639987e-01
#> mu_B.index.MRI      1.0343138 3.122716e-01
#> sigma2_prev.CT      0.6260525 1.017571e-01
#> sigma2_prev.MRI     0.6260525 1.017571e-01
#> sigma2_A.index.CT   0.8418377 2.174135e-01
#> sigma2_A.index.MRI  0.8418377 2.174135e-01
#> sigma2_B.index.CT   1.1632812 3.014621e-01
#> sigma2_B.index.MRI  1.1632812 3.014621e-01
#> sigma_AB.index.CT   0.2680248 1.683303e-01
#> sigma_AB.index.MRI  0.2680248 1.683303e-01
#> rho_AB.index.CT     0.2708435 1.541600e-01
#> rho_AB.index.MRI    0.2708435 1.541600e-01
#> mu_A.ref            3.4844030 3.979323e-01
#> mu_B.ref           18.3550258 1.698887e+03
#> nu_prev.MRI         0.1489753 2.205362e-01
#> nu_A.index.MRI     -1.3893820 3.002888e-01
#> nu_B.index.MRI     -1.3496420 3.630471e-01
#> 
#> $sensspec
#>                type  Estimate conflevel  CI_Lower  CI_Upper
#> mu_prev.CT     Prev 0.5695185      0.95 0.5237096 0.6141675
#> mu_prev.MRI    Prev 0.6056019      0.95 0.5070881 0.6962217
#> mu_A.index.CT  Sens 0.9701323      0.95 0.9598575 0.9778379
#> mu_A.index.MRI Sens 0.8900520      0.95 0.8286342 0.9312830
#> mu_B.index.CT  Spec 0.9155956      0.95 0.8660560 0.9479148
#> mu_B.index.MRI Spec 0.7377514      0.95 0.6040270 0.8383978
#> mu_A.ref       Sens 0.9702407      0.95 0.9372899 0.9861338
#> mu_B.ref       Spec 1.0000000      0.95 0.0000000 1.0000000
#> 
#> $RutterGatsonis_recovered
#>       Lambda     Theta      beta sigma2_alpha sigma2_theta
#> CT  5.972558 0.7874873 0.1617064     2.515235     0.360784
#> MRI 3.221360 0.6567022 0.1617064     2.515235     0.360784
plot(schuetzLCA2,
     nudge_legend=-0.2,
     size="eb",
     col=c("red","black"))



schuetzLCA3 <- fitReitsmaSubgroupLCA(data=schuetz,
                                     y11=TP,
                                     y10=FP,
                                     y01=FN,
                                     y00=TN,
                                     subgroup=test,
                                     variances="unequal",
                                     prev_variances="unequal",
                                     study=study)
#> Warning in sqrt(diag(cov)): NaNs produced
schuetzLCA3
#> 
#> Reitsma Subgroup LCA Model
#> --------------------------
#> 
#> Number of studies : 108 
#> Optimizer         : Converged 
#> Hessian           : Not positive definite 
#> Max |grad|        : 0.001991942 
#> -2 log likelihood : 16083.64 ( df = 16 )
#> AIC               : 16115.64 
#> BIC               : 16158.55 
#> 
#> Use summary() for parameter estimates.
plot(schuetzLCA3,
     nudge_legend=-0.2,
     size="eb",
     col=c("red","black"))



schuetzLCA4 <- fitReitsmaSubgroupLCA(data=schuetz,
                                     y11=TP,
                                     y10=FP,
                                     y01=FN,
                                     y00=TN,
                                     subgroup=test,
                                     variances="unequal",
                                     prev_variances="common",
                                     study=study)
#> Warning in sqrt(diag(cov)): NaNs produced

anova(schuetzLCA ,schuetzLCA2)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1  9 -8060.7                              
#> Model 2 12 -8046.3       3 28.727  2.555e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(schuetzLCA2,schuetzLCA3)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)  
#> Model 1 12 -8046.3                            
#> Model 2 16 -8041.8       4 9.0112    0.06082 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(schuetzLCA2,schuetzLCA4)
#>         Df  logLik Df.diff Chisq Pr(>Chisq)  
#> Model 1 12 -8046.3                           
#> Model 2 15 -8042.0       3 8.733    0.03306 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

schuetzLCA5 <- fitReitsmaSubgroupLCA(data=schuetz,
                                     y11=TP,
                                     y10=FP,
                                     y01=FN,
                                     y00=TN,
                                     subgroup=test,
                                     variances="unequal",
                                     sensspec_constrain = "sens",
                                     prev_variances="common",
                                     study=study)
schuetzLCA5
#> 
#> Reitsma Subgroup LCA Model
#> --------------------------
#> 
#> Number of studies : 108 
#> Optimizer         : Converged 
#> Hessian           : Positive definite 
#> Max |grad|        : 0.002178951 
#> -2 log likelihood : 16108.43 ( df = 14 )
#> AIC               : 16136.43 
#> BIC               : 16173.99 
#> 
#> Use summary() for parameter estimates.
anova(schuetzLCA5,schuetzLCA4)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1 14 -8054.2                              
#> Model 2 15 -8042.0       1 24.519  7.358e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

schuetzLCA6 <- fitReitsmaSubgroupLCA(data=schuetz,
                                     y11=TP,
                                     y10=FP,
                                     y01=FN,
                                     y00=TN,
                                     subgroup=test,
                                     variances="unequal",
                                     sensspec_constrain = "spec",
                                     prev_variances="common",
                                     study=study)
#> Warning in sqrt(diag(cov)): NaNs produced
schuetzLCA6
#> 
#> Reitsma Subgroup LCA Model
#> --------------------------
#> 
#> Number of studies : 108 
#> Optimizer         : Converged 
#> Hessian           : Not positive definite 
#> Max |grad|        : 0.002580734 
#> -2 log likelihood : 16096.9 ( df = 14 )
#> AIC               : 16124.9 
#> BIC               : 16162.44 
#> 
#> Use summary() for parameter estimates.
anova(schuetzLCA6,schuetzLCA4)
#>         Df  logLik Df.diff  Chisq Pr(>Chisq)    
#> Model 1 14 -8048.4                              
#> Model 2 15 -8042.0       1 12.979   0.000315 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

References

Liu, Y., Chen, Y., & Chu, H. (2015). A unification of models for meta-analysis of diagnostic accuracy studies without a gold standard. Biometrics, 71(2), 538–547.