MASTERfile.R
In esm-ispm-unibe-ch/rankingagreement:
rm(list=ls())
library(devtools)
library(nmadb)
library(readxl)
install_github("esm-ispm-unibe-ch/rankingagreement")
library(rankingagreement)
install_github("esm-ispm-unibe-ch/NMAJags")
library(NMAJags)
library(R2jags)
library(netmeta)
library(ircor)
install.packages("xlsx")
library(xlsx)
loadOrRun = function (filename, runfunction){
if(file.exists(filename)){
out = readRDS(filename)
}else{
out = runfunction()
}
return(out)
}
# access database or locally saved file
nmadb=loadOrRun("nmadb.RData",
function () {
ndb = getNMADB()
saveRDS(file = "./nmadb.RData", object = ndb)
return(ndb)}
)
#nmadb = getNMADB()
# nmadb_used = nmadb[nmadb$Verified=="True" & nmadb$Format!="iv" & (nmadb$Type.of.Outcome.=="Binary" | nmadb$Type.of.Outcome.=="Continuous"),]
# get IDs of networks separately for binary and continuous outcome; those with inverse-variance excluded
binaryIDs = nmadb[nmadb$Verified=="True" & nmadb$Type.of.Outcome.=="Binary" & nmadb$Format!="iv",]$Record.ID
continuousIDs = nmadb[nmadb$Verified=="True" & nmadb$Type.of.Outcome.=="Continuous" & nmadb$Format!="iv",]$Record.ID
# the following functions get the netmetas and the datas for all the networks with continuous outcomes
getContinuousNMAs = function() {
out = lapply(continuousIDs,
function(rid)
return(list(rid=rid,netobj=runnetmeta(rid)))
)
return(out)
}
continuousNetObjects = loadOrRun("continuousNetObjects.RData",
function () {
nmas = getContinuousNMAs()
saveRDS(file = "./continuousNetObjects.RData", object = nmas)
return(nmas)}
)
names(continuousNetObjects) <- continuousIDs
getContinuousDatasets = function() {
out = lapply(continuousIDs,
function(rid)
return(readByID(rid))
)
return(out)
}
continuousDatasets = loadOrRun("continuousDatasets.RData",
function () {
datas = getContinuousDatasets()
saveRDS(file = "./continuousDatasets.RData", object = datas)
return(datas)}
)
names(continuousDatasets) <- continuousIDs
# calculate ranking metrics for continuous outcome networks
getContinuous = function () {
out = lapply( 1:length(continuousNetObjects),
function(i) {
tryCatch({
nma = continuousNetObjects[[i]]$netobj
rid = continuousNetObjects[[i]]$rid
netd = continuousDatasets[[i]]
if (nmadb[nmadb$Record.ID==rid,]$Harmful.Beneficial.Outcome=="Beneficial"){
nmaranks = netmetaranks_B(nma,1000)
altnma = alternativenma(nma, small.values = "bad")
jagsranks = nmajagsranks_con_B(netd$data)
}
else {
nmaranks = netmetaranks_H(nma,1000)
altnma = alternativenma(nma)
jagsranks = nmajagsranks_con(netd$data)
}
return(list("no. treatments"=nma$n, "no. studies"=nma$k,"sample size"=sum(netd$data$n),
"ranking metrics"=cbind(nmaranks[order(as.numeric(rownames(nmaranks))),],
jagsranks, "Avg TE"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE,
"Avg TE ranks"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE_ranks,
"Avg Pscore"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$Pscoreaverage),
"Avg TE prec range"=(max(altnma$averages$seTE^2)-min(altnma$averages$seTE^2))/max(altnma$averages$seTE^2),
"Avg TE prec avg"=mean(altnma$averages$seTE^2)))
}, error=function(cond){
message(cond)
return(list(recid=rid,error=cond))
}
)
}
)
return(out)
}
pdf("traceplots_cont.pdf")
continuous_rm = loadOrRun("continuousRM.RData",
function () {
conts = getContinuous()
saveRDS(file = "./continuousRM.RData", object = conts)
return(conts)}
)
# continuous_rm = getContinuous()
dev.off()
names(continuous_rm) <- as.character(continuousIDs)
head(continuous_rm)
# the following functions get the netmetas and the datas for all the networks with continuous outcomes
getBinaryNMAs = function() {
out = lapply(binaryIDs,
function(rid)
return(list(rid=rid,netobj=runnetmeta(rid)))
)
return(out)
}
binaryNetObjects = loadOrRun("binaryNetObjects.RData",
function () {
nmas = getBinaryNMAs()
saveRDS(file = "./binaryNetObjects.RData", object = nmas)
return(nmas)}
)
names(binaryNetObjects) <- binaryIDs
getBinaryDatasets = function() {
out = lapply(binaryIDs,
function(rid)
return(readByID(rid))
)
return(out)
}
binaryDatasets = loadOrRun("binaryDatasets.RData",
function () {
datas = getBinaryDatasets()
saveRDS(file = "./binaryDatasets.RData", object = datas)
return(datas)}
)
names(binaryDatasets) <- binaryIDs
# calculate ranking metrics for binary outcome networks
getBinary = function () {
out = lapply( 1:length(binaryNetObjects),
function(i) {
tryCatch({
nma = binaryNetObjects[[i]]$netobj
rid = binaryNetObjects[[i]]$rid
netd = binaryDatasets[[i]]
if (nmadb[nmadb$Record.ID==rid,]$Harmful.Beneficial.Outcome=="Beneficial"){
nmaranks = netmetaranks_B(nma,1000)
altnma = alternativenma(nma, small.values = "bad")
jagsranks = nmajagsranks_bin_B(netd$data)
}
else {
nmaranks = netmetaranks_H(nma,1000)
altnma = alternativenma(nma)
jagsranks = nmajagsranks_bin(netd$data)
}
return(list("no. treatments"=nma$n, "no. studies"=nma$k,"sample size"=sum(netd$data$n),
"ranking metrics"=cbind(nmaranks[order(as.numeric(rownames(nmaranks))),],
jagsranks, "Avg TE"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE,
"Avg TE ranks"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE_ranks,
"Avg Pscore"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$Pscoreaverage),
"Avg TE prec range"=(max(altnma$averages$seTE^2)-min(altnma$averages$seTE^2))/max(altnma$averages$seTE^2),
"Avg TE prec avg"=mean(altnma$averages$seTE^2))) }, error=function(cond){
message(cond)
return(list(recid=rid,error=cond))
}
)
}
)
}
pdf("traceplots_bin.pdf")
binary_rm = loadOrRun("binaryRM.RData",
function () {
bin = getBinary()
saveRDS(file = "./binaryRM.RData", object = bin)
return(bin)}
)
# binary_rm = getBinary()
dev.off()
names(binary_rm) <- as.character(binaryIDs)
head(binary_rm)
# create lists with only ranks
con_ranks <- lapply(1:length(continuous_rm), function(i) continuous_rm[[i]][["ranking metrics"]][,grepl("rank",colnames(continuous_rm[[i]][["ranking metrics"]]))])
names(con_ranks) <- as.character(continuousIDs)
head(con_ranks)
bin_ranks <- lapply(1:length(binary_rm), function(i) binary_rm[[i]][["ranking metrics"]][,grepl("rank",colnames(binary_rm[[i]][["ranking metrics"]]))])
names(bin_ranks) <- as.character(binaryIDs)
head(bin_ranks)
# create lists with only ranking metrics values for spearman correlation
# con_values <- lapply(1:length(continuous_rm), function(i) continuous_rm[[i]][["ranking metrics"]][,c(1,3,5,7,8,9,11,13)])
# bin_values <- lapply(1:length(binary_rm), function(i) binary_rm[[i]][["ranking metrics"]][,c(1,3,5,7,8,9,11,13)])
# calculate kendall correlation
kendall_con <- lapply(con_ranks, cor, method="kendall")
names(kendall_con) <- as.character(continuousIDs)
head(kendall_con)
kendall_bin <- lapply(bin_ranks, cor, method="kendall")
names(kendall_bin) <- as.character(binaryIDs)
head(kendall_bin)
# calculate spearman correlation
spearman_con <- lapply(con_ranks, cor, method="spearman")
names(spearman_con) <- as.character(continuousIDs)
head(spearman_con)
spearman_bin <- lapply(bin_ranks, cor, method="spearman")
names(spearman_bin) <- as.character(binaryIDs)
head(spearman_bin)
# checks for P-score and SUCRA ranks
pSCOREvsSUCRA_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["Pscore_ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["Pscore_ranks","SUCRA_ranks"], digits = 3)))
names(pSCOREvsSUCRA_s) <- as.character(c(continuousIDs,binaryIDs))
res_pSCOREvsSUCRA_s <- paste0(summary(pSCOREvsSUCRA_s, digits = 3)["Median"], " (", summary(pSCOREvsSUCRA_s, digits = 3)["1st Qu."], ", ", summary(pSCOREvsSUCRA_s, digits = 3)["3rd Qu."], ")")
sum(pSCOREvsSUCRA_s>0.99)/length(pSCOREvsSUCRA_s) # % of networks with spearman correlation >0.99
print(pSCOREvsSUCRA_s[pSCOREvsSUCRA_s<0.99])
# prepare matrix to store results
results <- matrix(nrow = 4, ncol = 4,
dimnames = list(c("Spearman rho", "Kendall tau", "Yilmaz tauAP", "Average Overlap"),
c("pBV vs SUCRA", "SUCRA vs ATE", "pBV vs ATE", "SUCRA vs SUCRAjags")))
# save all pBV vs SUCRA in a vector separately for kendall, spearman and AP then store median and interquartile range
pBVvsSUCRA_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)))
names(pBVvsSUCRA_s) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)))
names(pBVvsSUCRA_k) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"pBV ranks"], con_ranks[[i]][,"SUCRA_ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"pBV ranks"], bin_ranks[[i]][,"SUCRA_ranks"], decreasing=F), digits = 3)))
names(pBVvsSUCRA_AP) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"pBV ranks"]),order(con_ranks[[i]][,"SUCRA_ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"pBV ranks"]),order(bin_ranks[[i]][,"SUCRA_ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(pBVvsSUCRA_AO) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AO <- Filter(Negate(anyNA), pBVvsSUCRA_AO) ## exclude any NAs
head(pBVvsSUCRA_s)
head(pBVvsSUCRA_k)
head(pBVvsSUCRA_AP)
head(pBVvsSUCRA_AO)
results["Spearman rho","pBV vs SUCRA"] <- paste0(summary(pBVvsSUCRA_s, digits = 2)["Median"], " (", summary(pBVvsSUCRA_s, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_k, digits = 2)["Median"], " (", summary(pBVvsSUCRA_k, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_AP, digits = 2)["Median"], " (", summary(pBVvsSUCRA_AP, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_AO, digits = 2)["Median"], " (", summary(pBVvsSUCRA_AO, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_AO, digits = 2)["3rd Qu."], ")")
sum(pBVvsSUCRA_s<.9)
sum(pBVvsSUCRA_s<.9)/length(pBVvsSUCRA_s)
# save all SUCRA vs Avg TE in a vector separately for kendall, spearman and AP, then store median and interquartile range
SUCRAvsAvgTE_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)))
names(SUCRAvsAvgTE_s) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)))
names(SUCRAvsAvgTE_k) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"SUCRA_ranks"], con_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"SUCRA_ranks"], bin_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)))
names(SUCRAvsAvgTE_AP) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"SUCRA_ranks"]),order(con_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"SUCRA_ranks"]),order(bin_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(SUCRAvsAvgTE_AO) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AO <- Filter(Negate(anyNA), SUCRAvsAvgTE_AO) ## exclude any NAs
head(SUCRAvsAvgTE_s)
head(SUCRAvsAvgTE_k)
head(SUCRAvsAvgTE_AP)
head(SUCRAvsAvgTE_AO)
results["Spearman rho","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_s, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_s, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_k, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_k, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_AP, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_AP, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","SUCRA vs ATE"] <-paste0(summary(SUCRAvsAvgTE_AO, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_AO, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_AO, digits = 2)["3rd Qu."], ")")
sum(SUCRAvsAvgTE_s<.9)
sum(SUCRAvsAvgTE_s<.9)/length(SUCRAvsAvgTE_s)
# save all pBV vs Avg TE in a vector separately for kendall, spearman and AP then store median and interquartile range
pBVvsAvgTE_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["pBV ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["pBV ranks","Avg TE ranks"], digits = 3)))
names(pBVvsAvgTE_s) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["pBV ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["pBV ranks","Avg TE ranks"], digits = 3)))
names(pBVvsAvgTE_k) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"pBV ranks"], con_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"pBV ranks"], bin_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)))
names(pBVvsAvgTE_AP) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"pBV ranks"]),order(con_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"pBV ranks"]),order(bin_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(pBVvsAvgTE_AO) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AO <- Filter(Negate(anyNA), pBVvsAvgTE_AO) ## exclude any NAs
head(pBVvsAvgTE_s)
head(pBVvsAvgTE_k)
head(pBVvsAvgTE_AP)
head(pBVvsAvgTE_AO)
results["Spearman rho","pBV vs ATE"] <- paste0(summary(pBVvsAvgTE_s, digits = 2)["Median"], " (", summary(pBVvsAvgTE_s, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_k, digits = 2)["Median"], " (", summary(pBVvsAvgTE_k, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_AP, digits = 2)["Median"], " (", summary(pBVvsAvgTE_AP, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_AO, digits = 2)["Median"], " (", summary(pBVvsAvgTE_AO, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_AO, digits = 2)["3rd Qu."], ")")
sum(pBVvsAvgTE_s<.9)
sum(pBVvsAvgTE_s<.9)/length(pBVvsAvgTE_s)
# save all SUCRA vs SUCRA jags in a vector separately for kendall and spearman, then store proportion of network with values >0.9
SUCRAvsSUCRAjags_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)))
names(SUCRAvsSUCRAjags_s) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)))
names(SUCRAvsSUCRAjags_k) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"SUCRA_ranks"], con_ranks[[i]][,"SUCRAjags ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"SUCRA_ranks"], bin_ranks[[i]][,"SUCRAjags ranks"], decreasing=F), digits = 3)))
names(SUCRAvsSUCRAjags_AP) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"SUCRA_ranks"]),order(con_ranks[[i]][,"SUCRAjags ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"SUCRA_ranks"]),order(bin_ranks[[i]][,"SUCRAjags ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(SUCRAvsSUCRAjags_AO) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AO <- Filter(Negate(anyNA), SUCRAvsSUCRAjags_AO) ## exclude any NAs
#only consider networks with original measures OR or SMD
SUCRAvsSUCRAjags_s <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_s[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_s))])
SUCRAvsSUCRAjags_k <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_k[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_k))])
SUCRAvsSUCRAjags_AP <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_AP[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_AP))])
SUCRAvsSUCRAjags_AO <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_AO[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_AO))])
head(SUCRAvsSUCRAjags_s)
head(SUCRAvsSUCRAjags_k)
head(SUCRAvsSUCRAjags_AP)
head(SUCRAvsSUCRAjags_AO)
# results["Spearman rho","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_s, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_s, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_s, digits = 2)["3rd Qu."], ")")
# results["Kendall tau","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_k, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_k, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_k, digits = 2)["3rd Qu."], ")")
# results["Yilmaz tauAP","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_AP, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_AP, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_AP, digits = 2)["3rd Qu."], ")")
# results["Average Overlap","SUCRA vs SUCRAjags"] <-paste0(summary(SUCRAvsSUCRAjags_AO, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_AO, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_AO, digits = 2)["3rd Qu."], ")")
# check agreement between the two SUCRA and save these in results for the paper
sumSUCRAs <- summary(SUCRAvsSUCRAjags_s, digits = 2)
sumSUCRAk <- summary(SUCRAvsSUCRAjags_k, digits = 2)
sumSUCRAap <- summary(SUCRAvsSUCRAjags_AP, digits = 2)
sumSUCRAao <- summary(SUCRAvsSUCRAjags_AO, digits = 2)
results["Spearman rho","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAs["Median"], " (", sumSUCRAs["1st Qu."], ", ", sumSUCRAs["3rd Qu."], ")")
results["Kendall tau","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAk["Median"], " (", sumSUCRAk["1st Qu."], ", ", sumSUCRAk["3rd Qu."], ")")
results["Yilmaz tauAP","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAap["Median"], " (", sumSUCRAap["1st Qu."], ", ", sumSUCRAap["3rd Qu."], ")")
results["Average Overlap","SUCRA vs SUCRAjags"] <-paste0(sumSUCRAao["Median"], " (", sumSUCRAao["1st Qu."], ", ", sumSUCRAao["3rd Qu."], ")")
SUCRAs_90 <- sum(SUCRAvsSUCRAjags_s<0.9)/length(SUCRAvsSUCRAjags_s) # % of networks with spearman correlation <0.9
SUCRAk_90 <- sum(SUCRAvsSUCRAjags_k<0.9)/length(SUCRAvsSUCRAjags_k) # % of networks with kendall correlation <0.9
SUCRAap_90 <- sum(SUCRAvsSUCRAjags_AP<0.9)/length(SUCRAvsSUCRAjags_AP) # % of networks with Yilmaz AP correlation <0.9
SUCRAao_90 <- sum(SUCRAvsSUCRAjags_AO<0.9)/length(SUCRAvsSUCRAjags_AO)# % of networks with AO <0.9
# export matrix results in table
write.xlsx(results, "agreement results.xlsx")
# store number of treatments in each network in a vector
ntreat <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["no. treatments"]),sapply(1:length(binary_rm), function(i) binary_rm[[i]]["no. treatments"]))
names(ntreat) <- as.character(c(continuousIDs,binaryIDs))
head(ntreat)
# store sample sizes in each network in a vector
samplesizes <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["sample size"]),sapply(1:length(binary_rm), function(i) binary_rm[[i]]["sample size"]))
names(samplesizes) <- as.character(c(continuousIDs,binaryIDs))
head(samplesizes)
# sample size over num of treatments
samp_nt <- as.numeric(samplesizes)/as.numeric(ntreat)
names(samp_nt) <- as.character(c(continuousIDs,binaryIDs))
head(samp_nt)
# store relative range of precision for Avg TE in each network in a vector
AvgTEprec_range <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["Avg TE prec range"]), sapply(1:length(binary_rm), function(i) binary_rm[[i]]["Avg TE prec range"]))
# store average precision for Avg TE in each network in a vector
AvgTEprec_avg <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["Avg TE prec avg"]), sapply(1:length(binary_rm), function(i) binary_rm[[i]]["Avg TE prec avg"]))
names(AvgTEprec_range) <- as.character(c(continuousIDs,binaryIDs))
head(AvgTEprec_range)
names(AvgTEprec_avg) <- as.character(c(continuousIDs,binaryIDs))
head(AvgTEprec_avg)
### graphs to show relationship between correlations and networks measures (avg sample size, avg precision)
source("plots.R")
# For Table1 for publication
## NUMBER OF TREATMENTS IS INCORRECT IN DATABASE - MUST BE TAKEN FROM NETMETA RESULTS
finalDB <- nmadb[nmadb$Record.ID %in% binaryIDs | nmadb$Record.ID %in% continuousIDs,]
finalDB <- finalDB[,!grepl("Country..choice",colnames(finalDB))]
finalRM <- c(continuous_rm, binary_rm)
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["no. treatments"])))
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["no. studies"])))
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["sample size"])))
table(finalDB$Harmful.Beneficial.Outcome)
table(droplevels(finalDB$Type.of.Outcome.))
table(droplevels(finalDB$Year))
table(droplevels(finalDB$Journal.Name))
table(finalDB$Ranking.metric..choice.Probability.of.being.the.best.)
table(finalDB$Ranking.metric..choice.Rankograms.)
table(finalDB$Ranking.metric..choice.Median.rank.)
table(finalDB$Ranking.metric..choice.Mean.rank.)
table(finalDB$Ranking.metric..choice.SUCRA.)
table(finalDB$Ranking.metric..choice.P.score.)
table(finalDB$Ranking.metric..choice.Other.)
table(finalDB$Ranking.metric..choice.None.)
esm-ispm-unibe-ch/rankingagreement documentation built on July 7, 2020, 2:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
rm(list=ls())
library(devtools)
library(nmadb)
library(readxl)
install_github("esm-ispm-unibe-ch/rankingagreement")
library(rankingagreement)
install_github("esm-ispm-unibe-ch/NMAJags")
library(NMAJags)
library(R2jags)
library(netmeta)
library(ircor)
install.packages("xlsx")
library(xlsx)
loadOrRun = function (filename, runfunction){
if(file.exists(filename)){
out = readRDS(filename)
}else{
out = runfunction()
}
return(out)
}
# access database or locally saved file
nmadb=loadOrRun("nmadb.RData",
function () {
ndb = getNMADB()
saveRDS(file = "./nmadb.RData", object = ndb)
return(ndb)}
)
#nmadb = getNMADB()
# nmadb_used = nmadb[nmadb$Verified=="True" & nmadb$Format!="iv" & (nmadb$Type.of.Outcome.=="Binary" | nmadb$Type.of.Outcome.=="Continuous"),]
# get IDs of networks separately for binary and continuous outcome; those with inverse-variance excluded
binaryIDs = nmadb[nmadb$Verified=="True" & nmadb$Type.of.Outcome.=="Binary" & nmadb$Format!="iv",]$Record.ID
continuousIDs = nmadb[nmadb$Verified=="True" & nmadb$Type.of.Outcome.=="Continuous" & nmadb$Format!="iv",]$Record.ID
# the following functions get the netmetas and the datas for all the networks with continuous outcomes
getContinuousNMAs = function() {
out = lapply(continuousIDs,
function(rid)
return(list(rid=rid,netobj=runnetmeta(rid)))
)
return(out)
}
continuousNetObjects = loadOrRun("continuousNetObjects.RData",
function () {
nmas = getContinuousNMAs()
saveRDS(file = "./continuousNetObjects.RData", object = nmas)
return(nmas)}
)
names(continuousNetObjects) <- continuousIDs
getContinuousDatasets = function() {
out = lapply(continuousIDs,
function(rid)
return(readByID(rid))
)
return(out)
}
continuousDatasets = loadOrRun("continuousDatasets.RData",
function () {
datas = getContinuousDatasets()
saveRDS(file = "./continuousDatasets.RData", object = datas)
return(datas)}
)
names(continuousDatasets) <- continuousIDs
# calculate ranking metrics for continuous outcome networks
getContinuous = function () {
out = lapply( 1:length(continuousNetObjects),
function(i) {
tryCatch({
nma = continuousNetObjects[[i]]$netobj
rid = continuousNetObjects[[i]]$rid
netd = continuousDatasets[[i]]
if (nmadb[nmadb$Record.ID==rid,]$Harmful.Beneficial.Outcome=="Beneficial"){
nmaranks = netmetaranks_B(nma,1000)
altnma = alternativenma(nma, small.values = "bad")
jagsranks = nmajagsranks_con_B(netd$data)
}
else {
nmaranks = netmetaranks_H(nma,1000)
altnma = alternativenma(nma)
jagsranks = nmajagsranks_con(netd$data)
}
return(list("no. treatments"=nma$n, "no. studies"=nma$k,"sample size"=sum(netd$data$n),
"ranking metrics"=cbind(nmaranks[order(as.numeric(rownames(nmaranks))),],
jagsranks, "Avg TE"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE,
"Avg TE ranks"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE_ranks,
"Avg Pscore"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$Pscoreaverage),
"Avg TE prec range"=(max(altnma$averages$seTE^2)-min(altnma$averages$seTE^2))/max(altnma$averages$seTE^2),
"Avg TE prec avg"=mean(altnma$averages$seTE^2)))
}, error=function(cond){
message(cond)
return(list(recid=rid,error=cond))
}
)
}
)
return(out)
}
pdf("traceplots_cont.pdf")
continuous_rm = loadOrRun("continuousRM.RData",
function () {
conts = getContinuous()
saveRDS(file = "./continuousRM.RData", object = conts)
return(conts)}
)
# continuous_rm = getContinuous()
dev.off()
names(continuous_rm) <- as.character(continuousIDs)
head(continuous_rm)
# the following functions get the netmetas and the datas for all the networks with continuous outcomes
getBinaryNMAs = function() {
out = lapply(binaryIDs,
function(rid)
return(list(rid=rid,netobj=runnetmeta(rid)))
)
return(out)
}
binaryNetObjects = loadOrRun("binaryNetObjects.RData",
function () {
nmas = getBinaryNMAs()
saveRDS(file = "./binaryNetObjects.RData", object = nmas)
return(nmas)}
)
names(binaryNetObjects) <- binaryIDs
getBinaryDatasets = function() {
out = lapply(binaryIDs,
function(rid)
return(readByID(rid))
)
return(out)
}
binaryDatasets = loadOrRun("binaryDatasets.RData",
function () {
datas = getBinaryDatasets()
saveRDS(file = "./binaryDatasets.RData", object = datas)
return(datas)}
)
names(binaryDatasets) <- binaryIDs
# calculate ranking metrics for binary outcome networks
getBinary = function () {
out = lapply( 1:length(binaryNetObjects),
function(i) {
tryCatch({
nma = binaryNetObjects[[i]]$netobj
rid = binaryNetObjects[[i]]$rid
netd = binaryDatasets[[i]]
if (nmadb[nmadb$Record.ID==rid,]$Harmful.Beneficial.Outcome=="Beneficial"){
nmaranks = netmetaranks_B(nma,1000)
altnma = alternativenma(nma, small.values = "bad")
jagsranks = nmajagsranks_bin_B(netd$data)
}
else {
nmaranks = netmetaranks_H(nma,1000)
altnma = alternativenma(nma)
jagsranks = nmajagsranks_bin(netd$data)
}
return(list("no. treatments"=nma$n, "no. studies"=nma$k,"sample size"=sum(netd$data$n),
"ranking metrics"=cbind(nmaranks[order(as.numeric(rownames(nmaranks))),],
jagsranks, "Avg TE"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE,
"Avg TE ranks"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$TE_ranks,
"Avg Pscore"=altnma$averages[order(as.numeric(rownames(altnma$averages))),]$Pscoreaverage),
"Avg TE prec range"=(max(altnma$averages$seTE^2)-min(altnma$averages$seTE^2))/max(altnma$averages$seTE^2),
"Avg TE prec avg"=mean(altnma$averages$seTE^2))) }, error=function(cond){
message(cond)
return(list(recid=rid,error=cond))
}
)
}
)
}
pdf("traceplots_bin.pdf")
binary_rm = loadOrRun("binaryRM.RData",
function () {
bin = getBinary()
saveRDS(file = "./binaryRM.RData", object = bin)
return(bin)}
)
# binary_rm = getBinary()
dev.off()
names(binary_rm) <- as.character(binaryIDs)
head(binary_rm)
# create lists with only ranks
con_ranks <- lapply(1:length(continuous_rm), function(i) continuous_rm[[i]][["ranking metrics"]][,grepl("rank",colnames(continuous_rm[[i]][["ranking metrics"]]))])
names(con_ranks) <- as.character(continuousIDs)
head(con_ranks)
bin_ranks <- lapply(1:length(binary_rm), function(i) binary_rm[[i]][["ranking metrics"]][,grepl("rank",colnames(binary_rm[[i]][["ranking metrics"]]))])
names(bin_ranks) <- as.character(binaryIDs)
head(bin_ranks)
# create lists with only ranking metrics values for spearman correlation
# con_values <- lapply(1:length(continuous_rm), function(i) continuous_rm[[i]][["ranking metrics"]][,c(1,3,5,7,8,9,11,13)])
# bin_values <- lapply(1:length(binary_rm), function(i) binary_rm[[i]][["ranking metrics"]][,c(1,3,5,7,8,9,11,13)])
# calculate kendall correlation
kendall_con <- lapply(con_ranks, cor, method="kendall")
names(kendall_con) <- as.character(continuousIDs)
head(kendall_con)
kendall_bin <- lapply(bin_ranks, cor, method="kendall")
names(kendall_bin) <- as.character(binaryIDs)
head(kendall_bin)
# calculate spearman correlation
spearman_con <- lapply(con_ranks, cor, method="spearman")
names(spearman_con) <- as.character(continuousIDs)
head(spearman_con)
spearman_bin <- lapply(bin_ranks, cor, method="spearman")
names(spearman_bin) <- as.character(binaryIDs)
head(spearman_bin)
# checks for P-score and SUCRA ranks
pSCOREvsSUCRA_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["Pscore_ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["Pscore_ranks","SUCRA_ranks"], digits = 3)))
names(pSCOREvsSUCRA_s) <- as.character(c(continuousIDs,binaryIDs))
res_pSCOREvsSUCRA_s <- paste0(summary(pSCOREvsSUCRA_s, digits = 3)["Median"], " (", summary(pSCOREvsSUCRA_s, digits = 3)["1st Qu."], ", ", summary(pSCOREvsSUCRA_s, digits = 3)["3rd Qu."], ")")
sum(pSCOREvsSUCRA_s>0.99)/length(pSCOREvsSUCRA_s) # % of networks with spearman correlation >0.99
print(pSCOREvsSUCRA_s[pSCOREvsSUCRA_s<0.99])
# prepare matrix to store results
results <- matrix(nrow = 4, ncol = 4,
dimnames = list(c("Spearman rho", "Kendall tau", "Yilmaz tauAP", "Average Overlap"),
c("pBV vs SUCRA", "SUCRA vs ATE", "pBV vs ATE", "SUCRA vs SUCRAjags")))
# save all pBV vs SUCRA in a vector separately for kendall, spearman and AP then store median and interquartile range
pBVvsSUCRA_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)))
names(pBVvsSUCRA_s) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["pBV ranks","SUCRA_ranks"], digits = 3)))
names(pBVvsSUCRA_k) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"pBV ranks"], con_ranks[[i]][,"SUCRA_ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"pBV ranks"], bin_ranks[[i]][,"SUCRA_ranks"], decreasing=F), digits = 3)))
names(pBVvsSUCRA_AP) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"pBV ranks"]),order(con_ranks[[i]][,"SUCRA_ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"pBV ranks"]),order(bin_ranks[[i]][,"SUCRA_ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(pBVvsSUCRA_AO) <- as.character(c(continuousIDs,binaryIDs))
pBVvsSUCRA_AO <- Filter(Negate(anyNA), pBVvsSUCRA_AO) ## exclude any NAs
head(pBVvsSUCRA_s)
head(pBVvsSUCRA_k)
head(pBVvsSUCRA_AP)
head(pBVvsSUCRA_AO)
results["Spearman rho","pBV vs SUCRA"] <- paste0(summary(pBVvsSUCRA_s, digits = 2)["Median"], " (", summary(pBVvsSUCRA_s, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_k, digits = 2)["Median"], " (", summary(pBVvsSUCRA_k, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_AP, digits = 2)["Median"], " (", summary(pBVvsSUCRA_AP, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","pBV vs SUCRA"] <-paste0(summary(pBVvsSUCRA_AO, digits = 2)["Median"], " (", summary(pBVvsSUCRA_AO, digits = 2)["1st Qu."], ", ", summary(pBVvsSUCRA_AO, digits = 2)["3rd Qu."], ")")
sum(pBVvsSUCRA_s<.9)
sum(pBVvsSUCRA_s<.9)/length(pBVvsSUCRA_s)
# save all SUCRA vs Avg TE in a vector separately for kendall, spearman and AP, then store median and interquartile range
SUCRAvsAvgTE_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)))
names(SUCRAvsAvgTE_s) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["SUCRA_ranks","Avg TE ranks"], digits = 3)))
names(SUCRAvsAvgTE_k) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"SUCRA_ranks"], con_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"SUCRA_ranks"], bin_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)))
names(SUCRAvsAvgTE_AP) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"SUCRA_ranks"]),order(con_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"SUCRA_ranks"]),order(bin_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(SUCRAvsAvgTE_AO) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsAvgTE_AO <- Filter(Negate(anyNA), SUCRAvsAvgTE_AO) ## exclude any NAs
head(SUCRAvsAvgTE_s)
head(SUCRAvsAvgTE_k)
head(SUCRAvsAvgTE_AP)
head(SUCRAvsAvgTE_AO)
results["Spearman rho","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_s, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_s, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_k, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_k, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","SUCRA vs ATE"] <- paste0(summary(SUCRAvsAvgTE_AP, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_AP, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","SUCRA vs ATE"] <-paste0(summary(SUCRAvsAvgTE_AO, digits = 2)["Median"], " (", summary(SUCRAvsAvgTE_AO, digits = 2)["1st Qu."], ", ", summary(SUCRAvsAvgTE_AO, digits = 2)["3rd Qu."], ")")
sum(SUCRAvsAvgTE_s<.9)
sum(SUCRAvsAvgTE_s<.9)/length(SUCRAvsAvgTE_s)
# save all pBV vs Avg TE in a vector separately for kendall, spearman and AP then store median and interquartile range
pBVvsAvgTE_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["pBV ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["pBV ranks","Avg TE ranks"], digits = 3)))
names(pBVvsAvgTE_s) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["pBV ranks","Avg TE ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["pBV ranks","Avg TE ranks"], digits = 3)))
names(pBVvsAvgTE_k) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"pBV ranks"], con_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"pBV ranks"], bin_ranks[[i]][,"Avg TE ranks"], decreasing=F), digits = 3)))
names(pBVvsAvgTE_AP) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"pBV ranks"]),order(con_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"pBV ranks"]),order(bin_ranks[[i]][,"Avg TE ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(pBVvsAvgTE_AO) <- as.character(c(continuousIDs,binaryIDs))
pBVvsAvgTE_AO <- Filter(Negate(anyNA), pBVvsAvgTE_AO) ## exclude any NAs
head(pBVvsAvgTE_s)
head(pBVvsAvgTE_k)
head(pBVvsAvgTE_AP)
head(pBVvsAvgTE_AO)
results["Spearman rho","pBV vs ATE"] <- paste0(summary(pBVvsAvgTE_s, digits = 2)["Median"], " (", summary(pBVvsAvgTE_s, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_s, digits = 2)["3rd Qu."], ")")
results["Kendall tau","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_k, digits = 2)["Median"], " (", summary(pBVvsAvgTE_k, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_k, digits = 2)["3rd Qu."], ")")
results["Yilmaz tauAP","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_AP, digits = 2)["Median"], " (", summary(pBVvsAvgTE_AP, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_AP, digits = 2)["3rd Qu."], ")")
results["Average Overlap","pBV vs ATE"] <-paste0(summary(pBVvsAvgTE_AO, digits = 2)["Median"], " (", summary(pBVvsAvgTE_AO, digits = 2)["1st Qu."], ", ", summary(pBVvsAvgTE_AO, digits = 2)["3rd Qu."], ")")
sum(pBVvsAvgTE_s<.9)
sum(pBVvsAvgTE_s<.9)/length(pBVvsAvgTE_s)
# save all SUCRA vs SUCRA jags in a vector separately for kendall and spearman, then store proportion of network with values >0.9
SUCRAvsSUCRAjags_s <- c(sapply(1:length(con_ranks), function(i) round(spearman_con[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(spearman_bin[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)))
names(SUCRAvsSUCRAjags_s) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_k <- c(sapply(1:length(con_ranks), function(i) round(kendall_con[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)),
sapply(1:length(bin_ranks), function(i) round(kendall_bin[[i]]["SUCRA_ranks","SUCRAjags ranks"], digits = 3)))
names(SUCRAvsSUCRAjags_k) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AP <- c(sapply(1:length(con_ranks), function(i) round(tauAP_b(con_ranks[[i]][,"SUCRA_ranks"], con_ranks[[i]][,"SUCRAjags ranks"], decreasing=F), digits = 3)),
sapply(1:length(bin_ranks), function(i) round(tauAP_b(bin_ranks[[i]][,"SUCRA_ranks"], bin_ranks[[i]][,"SUCRAjags ranks"], decreasing=F), digits = 3)))
names(SUCRAvsSUCRAjags_AP) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AO <- c(sapply(1:length(con_ranks), function(i) if(as.numeric(continuous_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(con_ranks[[i]][,"SUCRA_ranks"]),order(con_ranks[[i]][,"SUCRAjags ranks"]),floor(as.numeric(continuous_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}),
sapply(1:length(bin_ranks), function(i) if(as.numeric(binary_rm[[i]]["no. treatments"])>5)
{round(averageoverlap(order(bin_ranks[[i]][,"SUCRA_ranks"]),order(bin_ranks[[i]][,"SUCRAjags ranks"]),floor(as.numeric(binary_rm[[i]]["no. treatments"])/2)), digits = 3)}
else {NA}))
names(SUCRAvsSUCRAjags_AO) <- as.character(c(continuousIDs,binaryIDs))
SUCRAvsSUCRAjags_AO <- Filter(Negate(anyNA), SUCRAvsSUCRAjags_AO) ## exclude any NAs
#only consider networks with original measures OR or SMD
SUCRAvsSUCRAjags_s <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_s[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_s))])
SUCRAvsSUCRAjags_k <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_k[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_k))])
SUCRAvsSUCRAjags_AP <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_AP[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_AP))])
SUCRAvsSUCRAjags_AO <- Filter(Negate(anyNA),SUCRAvsSUCRAjags_AO[match(nmadb[nmadb$Effect.Measure=="odds ratio" | nmadb$Effect.Measure=="standardized mean difference","Record.ID"], names(SUCRAvsSUCRAjags_AO))])
head(SUCRAvsSUCRAjags_s)
head(SUCRAvsSUCRAjags_k)
head(SUCRAvsSUCRAjags_AP)
head(SUCRAvsSUCRAjags_AO)
# results["Spearman rho","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_s, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_s, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_s, digits = 2)["3rd Qu."], ")")
# results["Kendall tau","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_k, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_k, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_k, digits = 2)["3rd Qu."], ")")
# results["Yilmaz tauAP","SUCRA vs SUCRAjags"] <- paste0(summary(SUCRAvsSUCRAjags_AP, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_AP, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_AP, digits = 2)["3rd Qu."], ")")
# results["Average Overlap","SUCRA vs SUCRAjags"] <-paste0(summary(SUCRAvsSUCRAjags_AO, digits = 2)["Median"], " (", summary(SUCRAvsSUCRAjags_AO, digits = 2)["1st Qu."], ", ", summary(SUCRAvsSUCRAjags_AO, digits = 2)["3rd Qu."], ")")
# check agreement between the two SUCRA and save these in results for the paper
sumSUCRAs <- summary(SUCRAvsSUCRAjags_s, digits = 2)
sumSUCRAk <- summary(SUCRAvsSUCRAjags_k, digits = 2)
sumSUCRAap <- summary(SUCRAvsSUCRAjags_AP, digits = 2)
sumSUCRAao <- summary(SUCRAvsSUCRAjags_AO, digits = 2)
results["Spearman rho","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAs["Median"], " (", sumSUCRAs["1st Qu."], ", ", sumSUCRAs["3rd Qu."], ")")
results["Kendall tau","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAk["Median"], " (", sumSUCRAk["1st Qu."], ", ", sumSUCRAk["3rd Qu."], ")")
results["Yilmaz tauAP","SUCRA vs SUCRAjags"] <- paste0(sumSUCRAap["Median"], " (", sumSUCRAap["1st Qu."], ", ", sumSUCRAap["3rd Qu."], ")")
results["Average Overlap","SUCRA vs SUCRAjags"] <-paste0(sumSUCRAao["Median"], " (", sumSUCRAao["1st Qu."], ", ", sumSUCRAao["3rd Qu."], ")")
SUCRAs_90 <- sum(SUCRAvsSUCRAjags_s<0.9)/length(SUCRAvsSUCRAjags_s) # % of networks with spearman correlation <0.9
SUCRAk_90 <- sum(SUCRAvsSUCRAjags_k<0.9)/length(SUCRAvsSUCRAjags_k) # % of networks with kendall correlation <0.9
SUCRAap_90 <- sum(SUCRAvsSUCRAjags_AP<0.9)/length(SUCRAvsSUCRAjags_AP) # % of networks with Yilmaz AP correlation <0.9
SUCRAao_90 <- sum(SUCRAvsSUCRAjags_AO<0.9)/length(SUCRAvsSUCRAjags_AO)# % of networks with AO <0.9
# export matrix results in table
write.xlsx(results, "agreement results.xlsx")
# store number of treatments in each network in a vector
ntreat <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["no. treatments"]),sapply(1:length(binary_rm), function(i) binary_rm[[i]]["no. treatments"]))
names(ntreat) <- as.character(c(continuousIDs,binaryIDs))
head(ntreat)
# store sample sizes in each network in a vector
samplesizes <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["sample size"]),sapply(1:length(binary_rm), function(i) binary_rm[[i]]["sample size"]))
names(samplesizes) <- as.character(c(continuousIDs,binaryIDs))
head(samplesizes)
# sample size over num of treatments
samp_nt <- as.numeric(samplesizes)/as.numeric(ntreat)
names(samp_nt) <- as.character(c(continuousIDs,binaryIDs))
head(samp_nt)
# store relative range of precision for Avg TE in each network in a vector
AvgTEprec_range <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["Avg TE prec range"]), sapply(1:length(binary_rm), function(i) binary_rm[[i]]["Avg TE prec range"]))
# store average precision for Avg TE in each network in a vector
AvgTEprec_avg <- c(sapply(1:length(continuous_rm), function(i) continuous_rm[[i]]["Avg TE prec avg"]), sapply(1:length(binary_rm), function(i) binary_rm[[i]]["Avg TE prec avg"]))
names(AvgTEprec_range) <- as.character(c(continuousIDs,binaryIDs))
head(AvgTEprec_range)
names(AvgTEprec_avg) <- as.character(c(continuousIDs,binaryIDs))
head(AvgTEprec_avg)
### graphs to show relationship between correlations and networks measures (avg sample size, avg precision)
source("plots.R")
# For Table1 for publication
## NUMBER OF TREATMENTS IS INCORRECT IN DATABASE - MUST BE TAKEN FROM NETMETA RESULTS
finalDB <- nmadb[nmadb$Record.ID %in% binaryIDs | nmadb$Record.ID %in% continuousIDs,]
finalDB <- finalDB[,!grepl("Country..choice",colnames(finalDB))]
finalRM <- c(continuous_rm, binary_rm)
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["no. treatments"])))
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["no. studies"])))
summary(as.numeric(sapply(1:length(finalRM), function(i) finalRM[[i]]["sample size"])))
table(finalDB$Harmful.Beneficial.Outcome)
table(droplevels(finalDB$Type.of.Outcome.))
table(droplevels(finalDB$Year))
table(droplevels(finalDB$Journal.Name))
table(finalDB$Ranking.metric..choice.Probability.of.being.the.best.)
table(finalDB$Ranking.metric..choice.Rankograms.)
table(finalDB$Ranking.metric..choice.Median.rank.)
table(finalDB$Ranking.metric..choice.Mean.rank.)
table(finalDB$Ranking.metric..choice.SUCRA.)
table(finalDB$Ranking.metric..choice.P.score.)
table(finalDB$Ranking.metric..choice.Other.)
table(finalDB$Ranking.metric..choice.None.)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.