---
title: "Meta-Regression"
format:
html:
toc: TRUE
# pdf
vignette: >
%\VignetteIndexEntry{Meta-Regression}
%\VignetteEncoding{UTF-8}
%\VignetteEngine{quarto::html}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(dtametaTMB)
```
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 $\boldsymbol{Z}_\mathrm{pred}$.
## How do we include covariates in the Reitsma model?
For sensitivity in study $i$, we have the number of diseased individuals testing positive: $y_{Ai} \sim \mathcal{B}(n_{Ai},\pi_{Ai})$.
Similarly for specificity, we have the number of non-diseased individuals testing negative: $y_{Bi} \sim \mathcal{B}(n_{Bi},\pi_{Bi})$.
Now we introduce a $p$-dimensional design-vector $\boldsymbol{z}_i$ 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 $\boldsymbol{z}_i$ 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} $$
```{r}
#| echo: TRUE
data("anticcp")
reitsmasub <- fitReitsmaSubgroup(data=anticcp,
TP=TP,
FP=FP,
FN=FN,
TN=TN,
study=study,
subgroup=generation)
reitsmasub
summary(reitsmasub)
```
## How do I get a summary plot of the Reitsma model?
```{r}
#| fig.height: 7
#| fig.width: 8
#| echo: TRUE
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).
```{r}
#| fig.height: 11
#| fig.width: 13
#| echo: TRUE
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,
```{r}
#| echo: TRUE
constrainA <- fitReitsmaSubgroup(data=anticcp,
TP=TP,
FP=FP,
FN=FN,
TN=TN,
study=study,
subgroup=generation,
constrain="sigma2_A")
summary(constrainA)$estimates
```
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,
```{r}
#| echo: TRUE
constrainsens <- fitReitsmaSubgroup(data=anticcp,
TP=TP,
FP=FP,
FN=FN,
TN=TN,
study=study,
subgroup=generation,
sensspec_constrain="sens")
summary(constrainsens)$estimates
```
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.
```{r}
anova(constrainsens,reitsmasub)
```
## How do I allow for different subgroup-specific random-effects (co-)variances in the Reitsma model?
```{r}
#| fig.height: 7
#| fig.width: 8
#| echo: TRUE
heteroskedastic <- fitReitsmaSubgroup(data=anticcp,
TP=TP,
FP=FP,
FN=FN,
TN=TN,
study=study,
subgroup=generation,
variances="unequal")
heteroskedastic
plot(heteroskedastic,
scale=0.01,
nudge_legend=-0.2,
size="se",
col=c("black","red"))
anova(reitsmasub,heteroskedastic)
```
## 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 $y_{i1} \sim \mathcal{B}(n_{i1},\pi_{i1})$.
Similarly, the number of non-diseased individuals who test positive is $y_{i2} \sim \mathcal{B}(n_{i2},\pi_{i2})$.
Now we introduce a $p$-dimensional design-vector $\boldsymbol{z}_i$ including study level covariates.
Consequently, at the study level, we have
$$\operatorname{logit}(\pi_{ij}) =
\left( \boldsymbol{z}_i^{\top} \boldsymbol{\theta}_i + \boldsymbol{z}_i^{\top} \boldsymbol{\alpha}_i x_{ij} \right)
\exp(-\boldsymbol{z}_i^{\top}\boldsymbol{\beta} x_{ij}), $$
$$ \boldsymbol{z}_i^{\top} \boldsymbol{\alpha}_i \sim \mathcal{N}(\boldsymbol{z}_i^{\top}\boldsymbol{\Lambda}, \sigma_\alpha^2), \quad \quad \boldsymbol{z}_i^{\top} \boldsymbol{\theta}_i \sim \mathcal{N}(\boldsymbol{z}_i^{\top}\boldsymbol{\Theta}, \sigma_\theta^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 $\boldsymbol{z}_i$ 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?
```{r}
#| echo: TRUE
#| results: 'hide'
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
```
```{r}
#| echo: TRUE
ruttergatsonissub
summary(ruttergatsonissub)
```
## How do I get a summary plot of the Rutter and Gatsonis subgroup model?
```{r}
#| fig.height: 6
#| fig.width: 8
#| echo: TRUE
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).
```{r}
#| fig.height: 13
#| fig.width: 13
#| echo: TRUE
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 $\sigma_\alpha^2$ to zero,
- `constrain="sigma2_theta"` sets $\sigma_\theta^2$ to zero,
- `constrain="accuracy"` assumes equal accuracy parameters across subgroups, i.e. $\xi_2=\xi_3=\dots=0$,
- `constrain="threshold"` assumes equal threshold parameters across subgroups, i.e. $\gamma_2=\gamma_3=\dots=0$,
- `constrain="shape"` assumes equal shape parameters across subgroups, i.e. $\delta_2=\delta_3=\dots=0$,
- `constrain="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 $\boldsymbol{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.
```{r}
#| echo: TRUE
# 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
summary(ruttergatsonisreg)
```
## 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.
```{r}
#| echo: TRUE
ruttergatsonissubfull <- fitRutterGatsonisSubgroup(data=RF2,
TP=TP,
FP=FP,
FN=FN,
TN=TN,
study=study,
subgroup=method,
constrain=NULL)
ruttergatsonissubfull
summary(ruttergatsonissubfull)
```
## How do I get the log likelihood, the AIC, and the BIC of a model?
```{r}
#| echo: TRUE
logLik(ruttergatsonissubfull)
AIC(ruttergatsonissubfull)
BIC(ruttergatsonissubfull)
logLik(ruttergatsonissub)
AIC(ruttergatsonissub)
BIC(ruttergatsonissub)
```
## How do I perform a likelihood ratio test of the two models?
The test below formally investigates whether the HSROC curve shapes are equal in all subgroups.
```{r}
#| echo: TRUE
anova(ruttergatsonissub,
ruttergatsonissubfull)
```
## 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?
```{r}
#| echo: TRUE
#| fig.height: 6
#| fig.width: 8
data(schuetz)
head(schuetz)
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)
schuetzreitsma
summary(schuetzreitsma)
plot(schuetzreitsma,
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(schuetzreitsma2,schuetzreitsma)
anova(schuetzreitsma3,schuetzreitsma)
```
## What are the results for the direct comparisons?
```{r}
#| echo: TRUE
#| fig.height: 6
#| fig.width: 8
schuetz2 <- subset(schuetz,indirect==0)
schuetzreitsma4 <- fitReitsmaSubgroup(data=schuetz2,
TP=TP,FP=FP,FN=FN,TN=TN,
study=study,
subgroup=test,
constrain="sigma2_A")
round(summary(schuetzreitsma4)$estimates,5)
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.*