---
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.*

