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