Nothing
# Functions for printing and summarising results from MBNMA models
# Author: Hugo Pedder
# Date created: 2019-07-24
get.beta.names <- function(mbnma) {
betanames <- list()
for (i in 1:4) {
if (!is.null(mbnma$model.arg[[paste0("beta.",i)]])) {
if (is.null(mbnma$model.arg$arg.params)) {
betanames[[paste0("beta.", i)]] <- paste0("beta.", i)
} else if (!is.null(mbnma$model.arg$arg.params)) {
temp <- mbnma$model.arg$arg.params$wrap.params[
mbnma$model.arg$arg.params$run.params==paste0("beta.", i)
]
betanames[[paste0("beta.", i)]] <- temp
}
}
}
return(betanames)
}
neatCrI <- function(vals, digits=3) {
vals <- signif(vals, digits = digits)
neat <- paste0(vals[2], " (", vals[1], ", ", vals[3], ")")
return(neat)
}
get.timeparam.str <- function(mbnma, beta=NULL, param="d") {
betanames <- get.beta.names(mbnma)
if (grepl("beta", betanames[[beta]])) {
temp <- strsplit(betanames[[beta]], split="\\.")[[1]][2]
} else {
temp <- betanames[[beta]]
}
match <- paste0("^", param, "\\.", temp, "(\\[[0-9]+\\])?")
sum.mat <- mbnma$BUGSoutput$summary[grepl(match, rownames(mbnma$BUGSoutput$summary)),
c(3,5,7)]
# Check for UME
# if (any(grepl("\\[[0-9]+,[0-9]+\\]", rownames(sum.mat)))) {
# sum.mat <- NULL
# }
# Check if size of sum.mat is too great for printing
if (is.matrix(sum.mat)) {
if (nrow(sum.mat) > length(mbnma$network$treatments)) {
sum.mat <- NULL
}
}
if (length(sum.mat)>0) {
tab.str <- c()
if (is.matrix(sum.mat)) {
if (any(grepl("^d\\..+\\[1\\]", rownames(sum.mat)[1]) |
grepl("^D\\..+\\[1\\]", rownames(sum.mat)[1]))) {
tab.str <- paste(tab.str,
paste(rownames(sum.mat)[1], "Network reference",
sep="\t"),
sep="\n"
)
count <- 2
} else {
tab.str <- paste(tab.str,
paste(rownames(sum.mat)[1],
neatCrI(sum.mat),
sep="\t"),
sep="\n"
)
count <- 1
}
for (i in count:nrow(sum.mat)) {
tab.str <- paste(tab.str,
paste(rownames(sum.mat)[i], neatCrI(sum.mat[i,]),
sep="\t"),
sep="\n"
)
}
} else if (is.vector(sum.mat)) {
tab.str <- paste(tab.str,
paste(rownames(mbnma$BUGSoutput$summary)[grepl(match, rownames(mbnma$BUGSoutput$summary))],
neatCrI(sum.mat),
sep="\t\t"),
sep="\n"
)
}
return(tab.str)
} else {
return(NULL)
}
}
treat.str <- function(mbnma, digits=4, ...) {
fun <- mbnma$model.arg$fun
treat.df <- mbnma$BUGSoutput$summary
treat.sect <- c()
for (i in seq_along(fun$apool)) {
cat(crayon::underline(crayon::bold(paste0(names(fun$apool)[i], " parameter\n"))))
paramdet <- vector()
if (fun$apool[i]=="rel") {
paramdet <- append(paramdet, paste0("Pooling: relative effects"))
} else if (fun$apool[i]=="abs") {
paramdet <- append(paramdet, paste0("Pooling: absolute effects"))
}
if (is.numeric(fun$amethod[i])) {
paramdet <- append(paramdet, paste0("Assigned a numeric value: ", fun$amethod[i]))
} else {
if (fun$amethod[i]=="common") {
paramdet <- append(paramdet, paste0("Method: common treatment effects"))
} else if (fun$amethod[i]=="random") {
paramdet <- append(paramdet, paste0("Method: random treatment effects"))
}
}
if (names(fun$apool)[i] %in% c("et50", "hill")) {
if (TRUE %in% fun$p.expon) {
paramdet <- append(paramdet, "Parameter modelled on exponential scale to ensure it takes positive values on the natural scale")
}
}
if (names(fun$amethod)[i] %in% names(mbnma$model.arg$class.effect)) {
paramdet <- append(paramdet, "Class effects modelled for this parameter")
}
if (any(c(names(fun$apool)[i], TRUE) %in% mbnma$model.arg$UME)) {
paramdet <- append(paramdet, "Unrelated Mean Effect results modelled for this parameter.\nToo many parameters to display in this summary")
cat(paste(paramdet, collapse="\n"))
} else {
cat(paste(paramdet, collapse="\n"))
# Select summary rows that correspond to parameter of interest
index <- grepl(paste0("^", names(fun$apool)[i]), rownames(treat.df)) |
grepl(paste0("^d\\.", i), rownames(treat.df)) |
grepl(paste0("^beta\\.", i), rownames(treat.df))
param.df <- data.frame(treat=mbnma$network$treatments,
param=rownames(treat.df)[index],
median=treat.df[index,'50%'],
l95=treat.df[index,'2.5%'],
u95=treat.df[index,'97.5%']
)
rownames(param.df) <- NULL
print(knitr::kable(param.df, col.names = c("Treatment", "Parameter", "Median", "2.5%", "97.5%"), digits=digits, ...))
}
# Add between-study SD for parameter
if (fun$amethod[i]=="random") {
cat("\nBetween-study SD modelled for this parameter:")
index <- grepl(paste0("^sd\\.", names(fun$apool)[i]), rownames(treat.df)) |
grepl(paste0("^sd\\.d\\.", i), rownames(treat.df)) |
grepl(paste0("^sd\\.beta\\.", i), rownames(treat.df))
param.df <- data.frame(param=rownames(treat.df)[index],
median=treat.df[index,'50%'],
l95=treat.df[index,'2.5%'],
u95=treat.df[index,'97.5%']
)
rownames(param.df) <- NULL
print(knitr::kable(param.df, col.names = c("Parameter", "Median", "2.5%", "97.5%"), digits=digits, ...))
}
cat("\n\n")
}
}
class.str <- function(mbnma, digits=4, ...) {
if (length(mbnma$model.arg$class.effect)>0) {
treat.df <- mbnma$BUGSoutput$summary
classes <- mbnma$model.arg$class.effect
cat(crayon::bold(crayon::underline("Class Effects\n")))
for (i in seq_along(classes)) {
cat(paste0("\nClass effects for ", names(classes)[i], "\n"))
if ("common" %in% classes[[i]]) {
cat("Common (fixed) class effects")
} else if ("random" %in% classes[[i]]) {
cat("Random (exchangeable) class effects")
}
# Select summary rows that correspond to parameter of interest
index <- grepl(paste0("^", toupper(names(classes)[i])), rownames(treat.df)) |
grepl(paste0("^D\\.", which(mbnma$model.arg$fun$params==names(classes)[i])), rownames(treat.df))
param.df <- data.frame(treat=mbnma$network$classes,
param=rownames(treat.df)[index],
median=treat.df[index,'50%'],
l95=treat.df[index,'2.5%'],
u95=treat.df[index,'97.5%']
)
rownames(param.df) <- NULL
print(knitr::kable(param.df, col.names = c("Class", "Parameter", "Median", "2.5%", "97.5%"), digits=digits, ...))
# Add between-study SD for parameter
if (classes[[i]]=="random") {
cat("\nWithin-class SD modelled for this parameter:")
index <- grepl(paste0("^sd\\.", toupper(names(classes)[i])), rownames(treat.df)) |
grepl(paste0("^sd\\.D\\.", which(mbnma$model.arg$fun$params==names(classes)[i])), rownames(treat.df))
param.df <- data.frame(param=rownames(treat.df)[index],
median=treat.df[index,'50%'],
l95=treat.df[index,'2.5%'],
u95=treat.df[index,'97.5%']
)
rownames(param.df) <- NULL
print(knitr::kable(param.df, col.names = c("Parameter", "Median", "2.5%", "97.5%"), digits=digits, ...))
}
cat("\n")
}
}
}
overall.str <- function(mbnma) {
# Print title
title <- crayon::bold("========================================\nTime-course MBNMA\n========================================\n")
# Prep time-course function
timefun <- mbnma$model.arg$fun$name
if (timefun %in% c("poly", "fpoly")) {
timefun <- paste0(timefun, " (degree = ", mbnma$model.arg$fun$nparam, ")")
}
if (timefun=="ns") {
timefun <- paste0("Natural cubic spline (knots = ", paste(mbnma$model.arg$fun$knots, collapse=", "), ")")
} else if (timefun=="ls") {
timefun <- paste0("Piecewise linear spline (knots = ", paste(mbnma$model.arg$fun$knots, collapse=", "), ")")
} else if (timefun=="bs") {
timefun <- paste0("B-spline (knots = ", paste(mbnma$model.arg$fun$knots, collapse=", "), "; degree = ",
mbnma$model.arg$fun$degree, ")")
}
# Print time-course function
overall.sect <- paste("Time-course function:", timefun, sep=" ")
# Add info on intercept
if (mbnma$model.arg$intercept==TRUE) {
overall.sect <- paste(overall.sect, "Data modelled with intercept", sep="\n")
} else if (mbnma$model.arg$intercept==FALSE) {
overall.sect <- paste(overall.sect, "Data modelled without intercept (change from baseline data assumed)", sep="\n")
}
# positive.scale
if (mbnma$model.arg$positive.scale==TRUE) {
overall.sect <- paste(overall.sect, "Responses restricted to taking positive values", sep="\n")
}
overall.sect <- paste(title, overall.sect, sep="\n")
cat(overall.sect)
}
cor.str <- function(mbnma, digits=4) {
if (!is.null(mbnma$model.arg$rho)) {
cat(paste0("\n", crayon::underline(crayon::bold("Correlation between time points")), "\n"))
cat(paste("Covariance structure:", mbnma$model.arg$covar, "\n", sep=" "))
if (is.numeric(mbnma$model.arg$rho)) {
cat(paste("Rho assigned a numeric value:", mbnma$model.arg$rho, sep=" "))
} else if (is.character(mbnma$model.arg$rho)) {
if ("rho" %in% mbnma$parameters.to.save) {
cat("Rho estimated from the data:")
rho <- mbnma$BUGSoutput$summary["rho",]
rho <- data.frame(param="rho",
median=rho['50%'],
l95=rho['2.5%'],
u95=rho['97.5%']
)
rownames(rho) <- NULL
print(knitr::kable(rho, col.names = c("Parameter", "Median", "2.5%", "97.5%"), digits=digits))
} else {
cat("<rho not monitored in parameters.to.save>")
}
}
cat("\n\n")
}
}
modfit.str <- function(x) {
totresdev.str <- c()
head <- crayon::bold("#### Model Fit Statistics ####\n")
# pD
pd.str <- "Effective number of parameters:"
if (x$model.arg$pd=="pv") {
pd <- "pD (pV) calculated using the rule, pD = var(deviance)/2 ="
} else if (x$model.arg$pd=="plugin") {
pd <- "pD calculated using the plug-in method ="
} else if (x$model.arg$pd=="pd.kl") {
pd <- "pD calculated using the Kullback-Leibler divergence ="
} else if (x$model.arg$pd=="popt") {
pd <- "pD calculated using an optimism adjustment ="
}
pd.str <- paste(pd.str, paste(pd, round(x$BUGSoutput$pD,0), sep=" "), sep="\n")
# Deviance
dev <- x$BUGSoutput$summary[
rownames(x$BUGSoutput$summary)=="deviance", 5]
dev.str <- paste("Deviance =", round(dev, 0), sep=" ")
# Totresdev
if ("totresdev" %in% x$parameters.to.save) {
totresdev <- round(
x$BUGSoutput$summary[
rownames(x$BUGSoutput$summary)=="totresdev", 5],
0)
} else {
totresdev <- "NOT MONITORED IN MODEL"
}
totresdev.str <- paste("Residual deviance =", totresdev, sep=" ")
dic <- x$BUGSoutput$DIC
dic.str <- paste("Deviance Information Criterion (DIC) =", round(dic, 0), "\n", sep=" ")
modfit.sect <- paste(head, pd.str, dev.str, totresdev.str, dic.str, sep="\n")
return(modfit.sect)
}
# Gets summary columns for parameters
getsum <- function(params, mbnma) {
if (length(params)==1) {
rows <- mbnma$BUGSoutput$summary[grepl(params, rownames(mbnma$BUGSoutput$summary)),]
} else if (length(params)==2) {
rows <- mbnma$BUGSoutput$summary[
grepl(params[1], rownames(mbnma$BUGSoutput$summary)) |
grepl(params[2], rownames(mbnma$BUGSoutput$summary))
,]
}
if (is.matrix(rows)) {
rows <- rows[,colnames(rows) %in% c("mean", "2.5%", "97.5%")]
} else if (is.vector(rows)) {
rows <- rows[names(rows) %in% c("mean", "2.5%", "97.5%")]
}
rows <- round(rows, digits = max(3, getOption("digits")-3))
return(rows)
}
rhat.warning <- function(mbnma, cutoff=1.02) {
rhats <- mbnma$BUGSoutput$summary[,colnames(mbnma$BUGSoutput$summary)=="Rhat"]
rhats <- names(rhats)[rhats>cutoff]
if (length(rhats)>0) {
msg <- paste0("The following parameters have Rhat values > ",
cutoff,
" which could be due to convergence issues:\n")
cat(crayon::red((paste0(msg, paste(rhats, collapse="\n")))))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.