Bivariate meta-analysis of sensitivity and specificity

## Compare INLA manual, section 3.4 for the model formulation

##

## Call using the bivariate joint prior as described in

## Paul et al; Stat Med 2010; 29:1325-1339

library(INLA)

data(BivMetaAnalysis)

# Twice the number of diagnostic tests

# (meaning number of sensitivity AND specificity values)

N2 <- dim(BivMetaAnalysis)[1]

# the first two elements of param specify the shape and rate parameters for

# log gamma prior for the first log precision, the second two for the second

# log precision and the third two for the normal prior of the Fisher-transformed

# correlation parameter.

formula <- Y~f(diid, model="2diid", param=c(0.25, 0.025, 0.25, 0.025, 0, 0.4), n=N2) +

lag.tp + lag.tn + ct.tp + ct.tn + mr.tp + mr.tn -1

model = inla(formula, family="binomial", data=BivMetaAnalysis, Ntrials=N, verbose=T)

## get the summary estimates for sensitivity and specificity

## by applying a backtransformation using the invlogit

invlogit <- function(x){return(1/(1+exp(-x)))}

# LAG

round(invlogit(model$summary.fixed["lag.tp",c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(model$summary.fixed["lag.tn",c("0.5quant","0.025quant", "0.975quant")]),2)

# CT

round(invlogit(model$summary.fixed["ct.tp", c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(model$summary.fixed["ct.tn", c("0.5quant","0.025quant", "0.975quant")]),2)

# MR

round(invlogit(model$summary.fixed["mr.tp", c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(model$summary.fixed["mr.tn", c("0.5quant","0.025quant", "0.975quant")]),2)

########################################################

########################################################

## Call using the Wishart prior

## (Note: The original model type "2diidwishart" is deprecated,

## so that we use the model "iid2d")

# Twice the number of diagnostic tests

# (meaning number of sensitivity AND specificity values)

N2 <- dim(BivMetaAnalysis)[1]

# Rearrange the data from an alternating structure to a block structure.

# The first block contains the values for TP (lag.tp, ct.tp, mr.tp)

# and the second block the values for TN (lag.tn, ct.tn, mr.tn).

BivMetaAnalysis_Wishart <- rbind(BivMetaAnalysis[seq(1,N2,by=2),], BivMetaAnalysis[seq(2,N2,by=2),])

# use a new index

BivMetaAnalysis_Wishart$diid <- 1:N2

formulaWish=Y~f(diid, model="iid2d", n=N2, hyper=list(prec1=list(prior="wishart2d", param=c(4,1,2, 0.1))))+ lag.tp + lag.tn + ct.tp + ct.tn + mr.tp + mr.tn -1

modelWish=inla(formulaWish, family="binomial", data=BivMetaAnalysis_Wishart, Ntrials=N, verbose=TRUE)

## get the summary estimates for sensitivity and specificity

## by applying a backtransformation using the invlogit

invlogit <- function(x){return(1/(1+exp(-x)))}

# LAG

round(invlogit(modelWish$summary.fixed["lag.tp",c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(modelWish$summary.fixed["lag.tn",c("0.5quant","0.025quant", "0.975quant")]),2)

# CT

round(invlogit(modelWish$summary.fixed["ct.tp", c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(modelWish$summary.fixed["ct.tn", c("0.5quant","0.025quant", "0.975quant")]),2)

# MR

round(invlogit(modelWish$summary.fixed["mr.tp", c("0.5quant","0.025quant", "0.975quant")]),2)

round(invlogit(modelWish$summary.fixed["mr.tn", c("0.5quant","0.025quant", "0.975quant")]),2)