FunctionsToPlot.R
In htx-r/doseresNMA:
# I modified absolutePredplot() && sucraPlot() works when we have a class effect but ask Georgia why predicting on larger scale make a big difference
#!!!!!! histDRparam(), I might remove placebo from the figure
#!!!!!! add OR-dose plot (ask Georgia)
# Note: most plots are build assuming we got estimates per-drug (change when makeJAGSmodel is changed)
# These functions are needed to create the following
# 1. 3 Tables of characterstics
# 2. net plot
# 3. forest plot
# 4. Assess model fit plot
# 5. Absoluteresponse-dose curve
# 6. consistency-incosistency agreement plot
# 7. covergance table: rhat and Effective sample size
# 8. convergance histograms of coefficeints
# 9. Sucra plot
# 1. Tables of characterstics + # 2. net plot I code them directly in app.R (I use functions from BUGSnet and MBNMAdose)
# 3. forest plot: dose-response coefficients
forestplot <- function(x=x, ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,metareg=F,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- c('b1.','b2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
plotdata$doseparam <- rep(param.lab,each=nclass)
plotdata$drug <- rep(class.lab,length(param.lab))
plotdata <- plotdata[plotdata[,'drug']!=refclass.lab,]
}else{
ndrugs <- length(drug.lab)
# param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
plotdata$doseparam <- rep(param.lab,each=ndrugs)
plotdata$drug <- rep(drug.lab,length(param.lab))
plotdata <- plotdata[plotdata[,'drug']!=ref.lab,]
}
theme_set(
theme_minimal() +
theme(legend.position = "right",axis.text.x = element_text(size=12),axis.text.y = element_text(size=12))
)
g <- ggplot2::ggplot(plotdata, ggplot2::aes(y=plotdata$X50., x=plotdata$drug)) +
ggplot2::geom_point() +
ggplot2::geom_errorbar(ggplot2::aes(ymin=plotdata$X2.5., ymax=plotdata$X97.5.),width=0.4,size=0.5) +
ggplot2::coord_flip()
g <- g + ggplot2::facet_wrap(~doseparam, scales="free_x")
# Axis labels
g <- g + ggplot2::xlab("Drug") +
ggplot2::ylab("log (OR)")
g <-g+theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
# axis.text.x = element_text(face='bold',size=14),axis.text.y = element_text(face='bold',size=14),
axis.title.x=element_text(size=16,face = "bold"),axis.title.y=element_text(size=16,face = "bold"),
# plot.title = element_text(size = 16, face = "bold"),
strip.background =element_rect(fill="snow3"),
strip.text.x = element_text(size = 14))
return(g)
}
# 4. fit plot
fitplot<- function(x){
jagssamples <- as.mcmc(x)
samples = do.call(rbind, jagssamples) %>% data.frame()
#dev <- samples %>% select(., starts_with("dev"))
dev <- samples %>% select(., starts_with("resdev"))
totresdev <- samples$totresdev %>% mean()
pmdev <- colMeans(dev)
rhat <- samples %>% select(., starts_with("rhat"))
r <- samples %>% select(., starts_with("r."))
n <- samples %>% select(., starts_with("n"))
rtilde <- rhat %>%
colMeans() %>%
matrix(nrow=1, ncol=ncol(rhat)) %>%
data.frame() %>%
slice(rep(1:n(), each = nrow(rhat)))
pmdev_fitted <- 2*(r*log(r/rtilde)+(n-r)*log((n-r)/(n-rtilde)))[1,]
leverage = pmdev-pmdev_fitted
DIC = sum(leverage,na.rm = T) + totresdev
sign = sign(colMeans(r)-colMeans(rhat))
w = sign*sqrt(as.numeric(pmdev))
pD = sum(leverage,na.rm = T)
eq = function(x,c){c-x^2}
x=seq(-3, 3, 0.001)
c1=eq(x, c=rep(1, 6001))
plot(w, leverage, xlab=expression('w'[ik]), ylab=expression('leverage'[ik]),
ylim=c(0, max(c1+3, na.rm=TRUE)*1.15), xlim=c(-3,3))
points(x, ifelse(c1 < 0, NA, c1), lty=1, col="firebrick3", type="l")
points(x, ifelse(c1 < -1, NA, c1)+1, lty=2, col="chartreuse4", type="l")
points(x, ifelse(c1 < -2, NA, c1)+2, lty=3, col="mediumpurple3", type="l")
points(x, ifelse(c1 < -3, NA, c1)+3, lty=4, col="deepskyblue3", type="l")
text(2, 4.3, paste("pD=", round(pD, 2)), cex = 0.8)
text(2, 3.9, paste("Dres=", round(totresdev,2)), cex = 0.8)
text(2, 3.5, paste("DIC=", round(DIC,2)), cex = 0.8)
}
# 5. absolute response VS dose - curve
absolutePredplot <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- 'p.drug'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
ndose <- sapply((1:nclass)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(class.lab[class.lab!=refclass.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}else{
ndrugs <- length(drug.lab)
param.lab <- 'p.drug.'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
ndose <- sapply((1:ndrugs)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}
plotdata$p.placebo<- exp(x$BUGSoutput$mean$Z)/(1+exp(x$BUGSoutput$mean$Z))
plotdata$lp.ci <- exp(x$BUGSoutput$summary['Z','2.5%'])/(1+exp(x$BUGSoutput$summary['Z','2.5%']))
plotdata$up.ci <- exp(x$BUGSoutput$summary['Z','97.5%'])/(1+exp(x$BUGSoutput$summary['Z','97.5%']))
# plot arguments
ylim <- round(max(plotdata$X97.5.),1)
theme_set(
theme_minimal() +
theme(legend.position = "right",axis.text.x = element_text(size=12),axis.text.y = element_text(size=12))
)
g <- ggplot2::ggplot(plotdata, ggplot2::aes(x=dose)) +
geom_line(aes(y=X50.,color='treatment'))+
geom_line(aes(y=p.placebo,color='placebo'))+
scale_color_manual(values=c('steelblue','coral4'),name="")+
guides(color = guide_legend(override.aes = list(size = 1.5)))+
ggplot2::geom_smooth(ggplot2::aes(x=dose,y=X50.,ymin=X2.5.,ymax=X97.5.),color='coral4',fill='coral1',
data=plotdata, stat="identity")+
geom_smooth(aes(x=dose, y=p.placebo, ymin=lp.ci,
ymax=up.ci),color='steelblue',fill='steelblue2',
data=plotdata, stat="identity")+
coord_cartesian(ylim = c(0, ylim))
g <- g + ggplot2::facet_wrap(~drug,scales = 'free_x')
g <- g + ggplot2::labs(y="Predicted absolute response", x="Dose")
g <- g + ggplot2::scale_linetype_manual(name="")
g+theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
# axis.text.x = element_text(face='bold',size=14),axis.text.y = element_text(face='bold',size=14),
axis.title.x=element_text(size=16,face = "bold"),axis.title.y=element_text(size=16,face = "bold"),
# plot.title = element_text(size = 16, face = "bold"),
strip.background =element_rect(fill="snow3"),
strip.text.x = element_text(size = 14))
}
# 6. consistency-incosistency agreement plot
cons.icons.Plot <- function(cons.model, incons.model,direct.lab){
myd <- data.frame(y=cons.model$BUGSoutput$summary[direct.lab,'mean'],x=incons.model$BUGSoutput$summary[direct.lab,'mean'])
ggplot(myd,aes(y=y,x=x))+
ylim(-0.5,0.5)+
xlim(-0.5,0.5)+
geom_point()+
geom_abline(slope=1,intercept=0)+
ggplot2::labs(y="Consistency model", x="Inconsistency model")+
theme_bw()
}
# 7. Histogram of DR parameters
histDRparam <- function(x=x,metareg=F,class.effect=F){ # , drug.lab=drug.lab
require(tidyr)
#param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
sims <- x$BUGSoutput$sims.array
sims.mat <- apply(sims,3L,c) # merge the 3 chains (columns into one)
sims.param <- sims.mat%>%data.frame()%>%select(starts_with(param.lab))
plotdata <- data.frame(sim=unlist(c(sims.param)),param=rep(colnames(sims.param),each=nrow(sims.param)))
g <- ggplot(plotdata, aes(x=sim)) +
geom_histogram(alpha=0.6, binwidth = 0.001,color="darkblue", fill="lightblue")
g <- g + ggplot2::facet_wrap(~plotdata$param,scales = 'free')
g <- g + ggplot2::labs(y="Frequency", x="iterations")
g+ theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+theme_bw()
}
# 8. table of rhat and effective sample size
converge.table <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
# paramlabs <- c(paste0('b1[',(1:nclass)[-ref.index],']'),paste0('b2[',(1:nclass)[-ref.index],']'))
paramlabs <- c(paste0('s.beta.1[',(1:nclass)[-ref.index],']'),paste0('s.beta.2[',(1:nclass)[-ref.index],']'))
labs <- class.lab[class.lab!=refclass.lab]
}else{
ndrugs <- length(drug.lab)
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
#paramlabs <- c(paste0('b1[',(1:ndrugs)[-ref.index],']'),paste0('b2[',(1:ndrugs)[-ref.index],']'))
paramlabs <- c(paste0('s.beta.1[',(1:ndrugs)[-ref.index],']'),paste0('s.beta.2[',(1:ndrugs)[-ref.index],']'))
labs <- drug.lab[drug.lab!=ref.lab]
}
xx <- x[["BUGSoutput"]][["summary"]][paramlabs,][,c('Rhat', 'n.eff')]
rownames(xx) <- c(paste0('b1.',labs),paste0('b2.',labs))
tibble(Parameters=rownames(xx), Rhat=xx[,'Rhat'],'Effective Sample Size'=as.integer(xx[,'n.eff']))
}
# 9. Sucra plot
sucraPlot <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- 'p.drug'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
ndose <- sapply((1:nclass)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(class.lab[class.lab!=refclass.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}else{
ndrugs <- length(drug.lab)
param.lab <- 'p.drug.'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
ndose <- sapply((1:ndrugs)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}
s.plot <- ggplot(data=plotdata, aes(x=dose, y=mean, group=drug)) +
geom_line(aes(color=drug)) +
theme_bw()+
ggplot2::labs(y="Predicted absolute response", x="Dose")
s.plot
}
# table of network characterstics per drug
library(dplyr)
tab.net <- function(data,drug,r){
drug.name <- levels(factor(eval(substitute(drug), data)))
data.drug <- sapply(drug.name, function(d){
data_per_drug <- data%>%filter(data$drug==d)
n.events <- data_per_drug%>%select(r)%>%sum()
n.studies <-data_per_drug%>%select(studyid)%>%unique()%>%nrow()
n.doses <-data_per_drug%>%nrow()
n.patients <- data_per_drug%>%select(n)%>%sum()
return(c(n.events,
n.patients,
n.studies,
n.doses
))
},simplify = F)
#
tbl <- do.call(rbind,data.drug)
colnames(tbl) <- c(
"Number of events",
"Number of patients",
"Number of studies",
"Number of non-zero doses")
add_tab <- as_tibble(tbl,rownames = NA)
#rownames(add_tab) <- rownames(tbl)
return(add_tab)
# to make nice tables in App
# add.tble <- tibble(Characteristic=c("Agent",
# "Number of events",
# "Number of patients",
# "Number of studies",
# "Number of non-zero doses"
# ),
# Value= c(drug.name,
# tbl[,1],
# tbl[,2],
# tbl[,3],
# tbl[,4]
# ))
}
# print the summary table of estimated coefficients
est.tab <- function(x=x, ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,metareg=F,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- c('b1.','b2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
est.data <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
doseparam.drug <- paste0(rep(param.lab,each=ndrugs),drug.lab)
est.tab <- with(est.data,cbind(
mean=mean,
sd=sd,
'%2.5'=X2.5.,
'%25'=X25.,
'%50'=X50.,
'%75'=X75.,
'%97.5'=X97.5.))
rownames(est.tab) <-doseparam.drug
}else{
ndrugs <- length(drug.lab)
# param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
est.data <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
doseparam.drug <- paste0(rep(param.lab,each=ndrugs),drug.lab)
est.tab <- with(est.data,cbind(
mean=mean,
sd=sd,
'%2.5'=X2.5.,
'%25'=X25.,
'%50'=X50.,
'%75'=X75.,
'%97.5'=X97.5.))
rownames(est.tab) <-doseparam.drug
}
return(est.tab)
}
htx-r/doseresNMA documentation built on Jan. 28, 2021, 5:32 a.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.
# I modified absolutePredplot() && sucraPlot() works when we have a class effect but ask Georgia why predicting on larger scale make a big difference
#!!!!!! histDRparam(), I might remove placebo from the figure
#!!!!!! add OR-dose plot (ask Georgia)
# Note: most plots are build assuming we got estimates per-drug (change when makeJAGSmodel is changed)
# These functions are needed to create the following
# 1. 3 Tables of characterstics
# 2. net plot
# 3. forest plot
# 4. Assess model fit plot
# 5. Absoluteresponse-dose curve
# 6. consistency-incosistency agreement plot
# 7. covergance table: rhat and Effective sample size
# 8. convergance histograms of coefficeints
# 9. Sucra plot
# 1. Tables of characterstics + # 2. net plot I code them directly in app.R (I use functions from BUGSnet and MBNMAdose)
# 3. forest plot: dose-response coefficients
forestplot <- function(x=x, ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,metareg=F,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- c('b1.','b2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
plotdata$doseparam <- rep(param.lab,each=nclass)
plotdata$drug <- rep(class.lab,length(param.lab))
plotdata <- plotdata[plotdata[,'drug']!=refclass.lab,]
}else{
ndrugs <- length(drug.lab)
# param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
plotdata$doseparam <- rep(param.lab,each=ndrugs)
plotdata$drug <- rep(drug.lab,length(param.lab))
plotdata <- plotdata[plotdata[,'drug']!=ref.lab,]
}
theme_set(
theme_minimal() +
theme(legend.position = "right",axis.text.x = element_text(size=12),axis.text.y = element_text(size=12))
)
g <- ggplot2::ggplot(plotdata, ggplot2::aes(y=plotdata$X50., x=plotdata$drug)) +
ggplot2::geom_point() +
ggplot2::geom_errorbar(ggplot2::aes(ymin=plotdata$X2.5., ymax=plotdata$X97.5.),width=0.4,size=0.5) +
ggplot2::coord_flip()
g <- g + ggplot2::facet_wrap(~doseparam, scales="free_x")
# Axis labels
g <- g + ggplot2::xlab("Drug") +
ggplot2::ylab("log (OR)")
g <-g+theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
# axis.text.x = element_text(face='bold',size=14),axis.text.y = element_text(face='bold',size=14),
axis.title.x=element_text(size=16,face = "bold"),axis.title.y=element_text(size=16,face = "bold"),
# plot.title = element_text(size = 16, face = "bold"),
strip.background =element_rect(fill="snow3"),
strip.text.x = element_text(size = 14))
return(g)
}
# 4. fit plot
fitplot<- function(x){
jagssamples <- as.mcmc(x)
samples = do.call(rbind, jagssamples) %>% data.frame()
#dev <- samples %>% select(., starts_with("dev"))
dev <- samples %>% select(., starts_with("resdev"))
totresdev <- samples$totresdev %>% mean()
pmdev <- colMeans(dev)
rhat <- samples %>% select(., starts_with("rhat"))
r <- samples %>% select(., starts_with("r."))
n <- samples %>% select(., starts_with("n"))
rtilde <- rhat %>%
colMeans() %>%
matrix(nrow=1, ncol=ncol(rhat)) %>%
data.frame() %>%
slice(rep(1:n(), each = nrow(rhat)))
pmdev_fitted <- 2*(r*log(r/rtilde)+(n-r)*log((n-r)/(n-rtilde)))[1,]
leverage = pmdev-pmdev_fitted
DIC = sum(leverage,na.rm = T) + totresdev
sign = sign(colMeans(r)-colMeans(rhat))
w = sign*sqrt(as.numeric(pmdev))
pD = sum(leverage,na.rm = T)
eq = function(x,c){c-x^2}
x=seq(-3, 3, 0.001)
c1=eq(x, c=rep(1, 6001))
plot(w, leverage, xlab=expression('w'[ik]), ylab=expression('leverage'[ik]),
ylim=c(0, max(c1+3, na.rm=TRUE)*1.15), xlim=c(-3,3))
points(x, ifelse(c1 < 0, NA, c1), lty=1, col="firebrick3", type="l")
points(x, ifelse(c1 < -1, NA, c1)+1, lty=2, col="chartreuse4", type="l")
points(x, ifelse(c1 < -2, NA, c1)+2, lty=3, col="mediumpurple3", type="l")
points(x, ifelse(c1 < -3, NA, c1)+3, lty=4, col="deepskyblue3", type="l")
text(2, 4.3, paste("pD=", round(pD, 2)), cex = 0.8)
text(2, 3.9, paste("Dres=", round(totresdev,2)), cex = 0.8)
text(2, 3.5, paste("DIC=", round(DIC,2)), cex = 0.8)
}
# 5. absolute response VS dose - curve
absolutePredplot <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- 'p.drug'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
ndose <- sapply((1:nclass)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(class.lab[class.lab!=refclass.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}else{
ndrugs <- length(drug.lab)
param.lab <- 'p.drug.'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
ndose <- sapply((1:ndrugs)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}
plotdata$p.placebo<- exp(x$BUGSoutput$mean$Z)/(1+exp(x$BUGSoutput$mean$Z))
plotdata$lp.ci <- exp(x$BUGSoutput$summary['Z','2.5%'])/(1+exp(x$BUGSoutput$summary['Z','2.5%']))
plotdata$up.ci <- exp(x$BUGSoutput$summary['Z','97.5%'])/(1+exp(x$BUGSoutput$summary['Z','97.5%']))
# plot arguments
ylim <- round(max(plotdata$X97.5.),1)
theme_set(
theme_minimal() +
theme(legend.position = "right",axis.text.x = element_text(size=12),axis.text.y = element_text(size=12))
)
g <- ggplot2::ggplot(plotdata, ggplot2::aes(x=dose)) +
geom_line(aes(y=X50.,color='treatment'))+
geom_line(aes(y=p.placebo,color='placebo'))+
scale_color_manual(values=c('steelblue','coral4'),name="")+
guides(color = guide_legend(override.aes = list(size = 1.5)))+
ggplot2::geom_smooth(ggplot2::aes(x=dose,y=X50.,ymin=X2.5.,ymax=X97.5.),color='coral4',fill='coral1',
data=plotdata, stat="identity")+
geom_smooth(aes(x=dose, y=p.placebo, ymin=lp.ci,
ymax=up.ci),color='steelblue',fill='steelblue2',
data=plotdata, stat="identity")+
coord_cartesian(ylim = c(0, ylim))
g <- g + ggplot2::facet_wrap(~drug,scales = 'free_x')
g <- g + ggplot2::labs(y="Predicted absolute response", x="Dose")
g <- g + ggplot2::scale_linetype_manual(name="")
g+theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
# axis.text.x = element_text(face='bold',size=14),axis.text.y = element_text(face='bold',size=14),
axis.title.x=element_text(size=16,face = "bold"),axis.title.y=element_text(size=16,face = "bold"),
# plot.title = element_text(size = 16, face = "bold"),
strip.background =element_rect(fill="snow3"),
strip.text.x = element_text(size = 14))
}
# 6. consistency-incosistency agreement plot
cons.icons.Plot <- function(cons.model, incons.model,direct.lab){
myd <- data.frame(y=cons.model$BUGSoutput$summary[direct.lab,'mean'],x=incons.model$BUGSoutput$summary[direct.lab,'mean'])
ggplot(myd,aes(y=y,x=x))+
ylim(-0.5,0.5)+
xlim(-0.5,0.5)+
geom_point()+
geom_abline(slope=1,intercept=0)+
ggplot2::labs(y="Consistency model", x="Inconsistency model")+
theme_bw()
}
# 7. Histogram of DR parameters
histDRparam <- function(x=x,metareg=F,class.effect=F){ # , drug.lab=drug.lab
require(tidyr)
#param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
sims <- x$BUGSoutput$sims.array
sims.mat <- apply(sims,3L,c) # merge the 3 chains (columns into one)
sims.param <- sims.mat%>%data.frame()%>%select(starts_with(param.lab))
plotdata <- data.frame(sim=unlist(c(sims.param)),param=rep(colnames(sims.param),each=nrow(sims.param)))
g <- ggplot(plotdata, aes(x=sim)) +
geom_histogram(alpha=0.6, binwidth = 0.001,color="darkblue", fill="lightblue")
g <- g + ggplot2::facet_wrap(~plotdata$param,scales = 'free')
g <- g + ggplot2::labs(y="Frequency", x="iterations")
g+ theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+theme_bw()
}
# 8. table of rhat and effective sample size
converge.table <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
# paramlabs <- c(paste0('b1[',(1:nclass)[-ref.index],']'),paste0('b2[',(1:nclass)[-ref.index],']'))
paramlabs <- c(paste0('s.beta.1[',(1:nclass)[-ref.index],']'),paste0('s.beta.2[',(1:nclass)[-ref.index],']'))
labs <- class.lab[class.lab!=refclass.lab]
}else{
ndrugs <- length(drug.lab)
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
#paramlabs <- c(paste0('b1[',(1:ndrugs)[-ref.index],']'),paste0('b2[',(1:ndrugs)[-ref.index],']'))
paramlabs <- c(paste0('s.beta.1[',(1:ndrugs)[-ref.index],']'),paste0('s.beta.2[',(1:ndrugs)[-ref.index],']'))
labs <- drug.lab[drug.lab!=ref.lab]
}
xx <- x[["BUGSoutput"]][["summary"]][paramlabs,][,c('Rhat', 'n.eff')]
rownames(xx) <- c(paste0('b1.',labs),paste0('b2.',labs))
tibble(Parameters=rownames(xx), Rhat=xx[,'Rhat'],'Effective Sample Size'=as.integer(xx[,'n.eff']))
}
# 9. Sucra plot
sucraPlot <- function(x=x,ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- 'p.drug'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(class.lab))==refclass.lab)
ndose <- sapply((1:nclass)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(class.lab[class.lab!=refclass.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}else{
ndrugs <- length(drug.lab)
param.lab <- 'p.drug.'
plotdata <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
ndose <- sapply((1:ndrugs)[-ref.index], function(i) ncol(plotdata %>% t()%>% data.frame()%>%select(ends_with(paste0('.',i,'.')))))
plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],times=ndose)
plotdata$dose <- unlist(sapply(1:length(ndose), function(i) 1:ndose[i],simplify = FALSE))
}
s.plot <- ggplot(data=plotdata, aes(x=dose, y=mean, group=drug)) +
geom_line(aes(color=drug)) +
theme_bw()+
ggplot2::labs(y="Predicted absolute response", x="Dose")
s.plot
}
# table of network characterstics per drug
library(dplyr)
tab.net <- function(data,drug,r){
drug.name <- levels(factor(eval(substitute(drug), data)))
data.drug <- sapply(drug.name, function(d){
data_per_drug <- data%>%filter(data$drug==d)
n.events <- data_per_drug%>%select(r)%>%sum()
n.studies <-data_per_drug%>%select(studyid)%>%unique()%>%nrow()
n.doses <-data_per_drug%>%nrow()
n.patients <- data_per_drug%>%select(n)%>%sum()
return(c(n.events,
n.patients,
n.studies,
n.doses
))
},simplify = F)
#
tbl <- do.call(rbind,data.drug)
colnames(tbl) <- c(
"Number of events",
"Number of patients",
"Number of studies",
"Number of non-zero doses")
add_tab <- as_tibble(tbl,rownames = NA)
#rownames(add_tab) <- rownames(tbl)
return(add_tab)
# to make nice tables in App
# add.tble <- tibble(Characteristic=c("Agent",
# "Number of events",
# "Number of patients",
# "Number of studies",
# "Number of non-zero doses"
# ),
# Value= c(drug.name,
# tbl[,1],
# tbl[,2],
# tbl[,3],
# tbl[,4]
# ))
}
# print the summary table of estimated coefficients
est.tab <- function(x=x, ref.lab='placebo',refclass.lab='placebo',drug.lab=drug.lab, class.lab=NULL,metareg=F,class.effect=F, ...){
if(class.effect){
nclass <- length(class.lab)
param.lab <- c('b1.','b2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
est.data <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
doseparam.drug <- paste0(rep(param.lab,each=ndrugs),drug.lab)
est.tab <- with(est.data,cbind(
mean=mean,
sd=sd,
'%2.5'=X2.5.,
'%25'=X25.,
'%50'=X50.,
'%75'=X75.,
'%97.5'=X97.5.))
rownames(est.tab) <-doseparam.drug
}else{
ndrugs <- length(drug.lab)
# param.lab <- c('b1.','b2.')
param.lab <-c('s.beta.1.','s.beta.2.')
if(metareg){
param.lab0 <- 'g.'
param.lab <- append(param.lab,param.lab0)
}
est.data <- x$BUGSoutput$summary %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
doseparam.drug <- paste0(rep(param.lab,each=ndrugs),drug.lab)
est.tab <- with(est.data,cbind(
mean=mean,
sd=sd,
'%2.5'=X2.5.,
'%25'=X25.,
'%50'=X50.,
'%75'=X75.,
'%97.5'=X97.5.))
rownames(est.tab) <-doseparam.drug
}
return(est.tab)
}
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.