## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(dtametaTMB)


## -----------------------------------------------------------------------------
#| echo: TRUE
data("anticcp")
reitsmasub <- fitReitsmaSubgroup(data=anticcp,
                                 TP=TP,
                                 FP=FP,
                                 FN=FN,
                                 TN=TN,
                                 study=study,
                                 subgroup=generation)
reitsmasub
summary(reitsmasub)


## -----------------------------------------------------------------------------
#| fig.height: 7
#| fig.width: 8
#| echo: TRUE
plot(reitsmasub,
     scale=0.01,
     nudge_legend=-0.2,
     size="se",
     col=c("black","red"))


## -----------------------------------------------------------------------------
#| fig.height: 11
#| fig.width: 13
#| echo: TRUE
forest(reitsmasub,subgroup_label="Generation")


## -----------------------------------------------------------------------------
#| echo: TRUE
constrainA <- fitReitsmaSubgroup(data=anticcp,
                                 TP=TP,
                                 FP=FP,
                                 FN=FN,
                                 TN=TN,
                                 study=study,
                                 subgroup=generation,
                                 constrain="sigma2_A")
summary(constrainA)$estimates


## -----------------------------------------------------------------------------
#| echo: TRUE
constrainsens <- fitReitsmaSubgroup(data=anticcp,
                                    TP=TP,
                                    FP=FP,
                                    FN=FN,
                                    TN=TN,
                                    study=study,
                                    subgroup=generation,
                                    sensspec_constrain="sens")
summary(constrainsens)$estimates


## -----------------------------------------------------------------------------
anova(constrainsens,reitsmasub)


## -----------------------------------------------------------------------------
#| 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")
plot(heteroskedastic,
     scale=0.01,
     nudge_legend=-0.2,
     size="se",
     col=c("black","red"))
anova(reitsmasub,heteroskedastic)


## -----------------------------------------------------------------------------
#| 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


## -----------------------------------------------------------------------------
#| echo: TRUE
ruttergatsonissub
summary(ruttergatsonissub)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| fig.height: 13
#| fig.width: 13
#| echo: TRUE
forest(ruttergatsonissub,subgroup_label = "Method")


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| echo: TRUE
ruttergatsonissubfull <- fitRutterGatsonisSubgroup(data=RF2,
                                                   TP=TP,
                                                   FP=FP,
                                                   FN=FN,
                                                   TN=TN,
                                                   study=study,
                                                   subgroup=method,
                                                   constrain=NULL)
ruttergatsonissubfull
summary(ruttergatsonissubfull)
ruttergatsonissubacc <- fitRutterGatsonisSubgroup(data=RF2,
                                                  TP=TP,
                                                  FP=FP,
                                                  FN=FN,
                                                  TN=TN,
                                                  study=study,
                                                  subgroup=method,
                                                  constrain=c("shape","accuracy"))
ruttergatsonissubacc
summary(ruttergatsonissubacc)



## -----------------------------------------------------------------------------
#| echo: TRUE
logLik(ruttergatsonissubfull)
AIC(ruttergatsonissubfull)
BIC(ruttergatsonissubfull)

logLik(ruttergatsonissub)
AIC(ruttergatsonissub)
BIC(ruttergatsonissub)

logLik(ruttergatsonissubacc)
AIC(ruttergatsonissubacc)
BIC(ruttergatsonissubacc)


## -----------------------------------------------------------------------------
#| echo: TRUE
anova(ruttergatsonissub,
      ruttergatsonissubfull)
anova(ruttergatsonissubacc,
      ruttergatsonissub)
anova(ruttergatsonissubacc,
      ruttergatsonissubfull)


## -----------------------------------------------------------------------------
#| 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,
                                     variances="common")
schuetzreitsma
summary(schuetzreitsma)
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
summary(schuetzreitsma1)
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)
anova(schuetzreitsma2,schuetzreitsma)
anova(schuetzreitsma3,schuetzreitsma)


## -----------------------------------------------------------------------------
#| 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")
schuetzreitsma4
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"))