The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
library(RaCE.NMA)
library(ggplot2)
library(dplyr)
library(cowplot)
library(gridExtra)
library(reshape2)In this vignette, we provide code to reproduce the analyses presented in the paper associated with this R package. Specifically, we reproduce two simulation studies (Sections 3.1 and 3.2) and a case study (Section 4). We assume knowledge of the associated R package functions; please refer to the “Tutorial” and “Reference” pages for more information. No figures are produced; see paper.
The following code reproduces Figure 1.
s <- 0.1
g1a<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
theme_bw()+labs(x=" ",y="Density",subtitle=expression(hat(sigma)~"=0.1"))+
theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
s <- 0.3
g1b<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
theme_bw()+labs(x="Posterior Intervention Effects",y=NULL,subtitle=expression(hat(sigma)~"=0.3"))+
theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
s <- 0.5
g1c<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
theme_bw()+labs(x=" ",y=NULL,subtitle=expression(hat(sigma)~"=0.5"))+
theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
grid.arrange(g1a,g1b,g1c,nrow=1)The following code reproduces Figure 2: First, we perform a simulation study. Second, we create the figure.
## Perform simulation study
set.seed(1)
results <- matrix(NA,nrow=0,ncol=6)
for(iter in 1:20){
for(J in c(6,12,18)){
for(K in c(J/3,2*J/3,J)){
for(s in c(0.1,0.3,0.5)){
if(K==J){
mu_hat <- 1:K
}else{
mu_hat <- sample(1:K,J,replace=T)
while(length(unique(mu_hat))<K){mu_hat <- sample(1:K,J,replace=T)}
}
mcmc <- mcmc_raceNMA(mu_hat = mu_hat,s=rep(s,J),mu=mean(mu_hat),
sigma0=max(1,4*sd(mu_hat)),nu0=mu_hat,
tau=1,iter=3000,nu_iter=2,chains=4,
burn_prop=0.5,thin=3,verbose = FALSE)
equal <- posterior_equal <- matrix(NA,nrow=J,ncol=J)
for(i in 1:(J-1)){for(j in (i+1):J){
equal[i,j] <- ifelse(mu_hat[i]==mu_hat[j],1,0) # test if treatments are truly equal in mean
posterior_equal[i,j] <- mean( # assess if treatments are rank-clustered
mcmc[,paste0("G",i)] == mcmc[,paste0("G",j)]
)
}}
results <- rbind(
results,
c(iter,J,K,s,
mean(posterior_equal[equal==1],na.rm=T),
mean(posterior_equal[equal==0],na.rm=T))
)
}
}
}
}## Plotting results from simulation study
results <- as.data.frame(results)
names(results) <- c("Iteration","J","K","s","Prob_Clustered","Prob_Distinct")
results$s<-factor(results$s, levels=c(.1,.3,.5),
labels=c(expression(hat(sigma)~"=0.1"),
expression(hat(sigma)~"=0.3"),
expression(hat(sigma)~"=0.5")))
results$J <- factor(results$J,levels=c(6,12,18),
labels=c(expression(J~"= 6"),
expression(J~"= 12"),
expression(J~"= 18")))
results_melt <- melt(results,id.vars=1:4)
results_melt <- results_melt[!is.nan(results_melt$value),] # drop "cluster" cases when K=J
ggplot(results_melt,aes(x=factor(K),y=value,
color=factor(variable,labels=c("Rank-Clustered","Distinct"))))+
facet_grid(s~J,scales="free_x",labeller=label_parsed)+
geom_boxplot(outlier.alpha=0.5,position="identity")+theme_bw()+
scale_color_manual(values=c("skyblue","darkblue"))+
labs(x="Number of Rank-Clusters, K",
y="Posterior Rank-Clustering Probability",
color=NULL)+
theme(legend.position = "bottom",
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank())The following code reproduces Figure 3.
J <- 4
mu_hat <- c(-1,0,1,0)
s <- c(.1,.1,.1,1)
set.seed(2)
data <- matrix(data = rnorm(10000*4,mean=mu_hat,sd=s), ncol = 4, byrow = TRUE)
forestplot_muhat(data = data, order_by_average = FALSE)The following code reproduces Table 1: First, we fit the RaCE model to the simulated data. Second, we calculate and display the values in Table 1.
mcmc <- mcmc_raceNMA(mu_hat = mu_hat, s = s,
iter = 50000, nu_iter = 2, chains = 4,
burn_prop = 0.5,thin = 1,verbose = FALSE)
posterior_equal <- matrix(NA,nrow=J,ncol=J)
for(i in 1:(J-1)){for(j in (i+1):J){
posterior_equal[i,j] <- mean(mcmc[,paste0("G",i)] == mcmc[,paste0("G",j)])
}}
round(posterior_equal,2)The following code reproduces Figure 4: First, we calculate the posterior rank-clustering probability under a traditional model (i.e., assumption of distinct ranks). Second, we calculate the posterior rank-clustering probability based on the RaCE model. Third, we produce the figure.
wang_ranks <- t(apply(data,1,function(mu){rank(mu)}))
wang_ranks_probs <- apply(wang_ranks,2,function(rank){
sapply(1:J,function(j){mean(rank==j)})
})
race_ranks <- t(apply( mcmc[,paste0("mu",1:J)], 1, function(mu){
rank(mu,ties.method="min")
}))
race_ranks_probs <- apply(race_ranks,2,function(rank){
sapply(1:J,function(j){mean(rank==j)})
})
g4a<-ggplot(melt(wang_ranks_probs),
aes(x=Var1,y=factor(Var2),fill=value)) +
geom_tile() + theme_minimal() +
scale_y_discrete(limits=rev) +
scale_x_continuous(breaks=1:4,limits=c(.5,4.5)) +
scale_fill_gradient(low="white",high="black",limits=c(0,1))+
labs(x="Rank",y="Treatment",fill="Probability")+
theme(panel.grid = element_blank(),legend.position = "bottom")+
geom_text(aes(x=Var1,y=factor(Var2),label=round(value,2)),
color=ifelse(melt(wang_ranks_probs)$value>0.4,"white","black"))
g4b<-ggplot(melt(race_ranks_probs),
aes(x=Var1,y=factor(Var2,levels=paste0("mu",1:J),
labels=paste0(1:J)),fill=value))+
geom_tile() + theme_minimal() +
scale_y_discrete(limits=rev) +
scale_x_continuous(breaks=1:4,limits=c(.5,4.5)) +
scale_fill_gradient(low="white",high="black",limits=c(0,1)) +
labs(x="Rank",y="Treatment",fill="Probability")+
theme(panel.grid = element_blank(),legend.position = "bottom")+
geom_text(aes(x=Var1,y=factor(Var2,levels=paste0("mu",1:J),labels=paste0(1:J)),
label=round(value,2)),
color=ifelse(melt(race_ranks_probs)$value>0.4,"white","black"))
plot_grid(plot_grid(g4a+theme(legend.position = "none"),
g4b+theme(legend.position = "none"),
labels = c('A', 'B'), label_size = 12),
get_plot_component(g4a, 'guide-box-bottom', return_all = TRUE),
nrow=2,rel_heights = c(.9,.1) )We first load the posteriors from Wang et al. (2022) and calculate the mean and covariances of the relative treatment effects.
data("wang_posterior") # loads posterior of non-baseline treatments
# define assumed mean and variance for baseline treatment (R-CHOP)
mu_hat_baseline <- 0
var_baseline <- min(apply(wang_posterior,2,var))/10
# calculate summary statistics, mu_hat and cov, for all treatments
mu_hat <- c( mu_hat_baseline, apply(wang_posterior,2,mean) )
cov <- cbind( c(var_baseline, rep(0, 10) ),
rbind(0, cov(wang_posterior)) )
# store treatment names
treatments <- c("R-CHOP", names(wang_posterior))The following code chunk reproduces Figure 5.
forestplot_muhat(data=wang_posterior,names=treatments[-1])Next, we fit the RaCE model to estimated relative treatment effects from Wang et al. (2022).
mcmc <- mcmc_raceNMA(mu_hat = mu_hat, cov = cov,
mu0 = mean(mu_hat), sigma0 = sqrt(10*var(mu_hat)),
iter = 50000, nu_iter = 5, chains = 4,
seed = 1, verbose = FALSE)The following code chunk reproduces Figure 6, Figure 7, and Table 2: First, we calculate the values underpinning each graphic. Then, we create each in turn.
# Figure 6
g9a <- clusterplot_ranks(data=cbind(0,wang_posterior),
names=treatments,label_ranks=1)+
theme(legend.position = "bottom")
g9b <- clusterplot_ranks(mcmc=mcmc,names=treatments,label_ranks = 1)
plot_grid(plot_grid(g9a+theme(legend.position = "none"),
g9b+theme(legend.position = "none"),
labels = c('A', 'B'), label_size = 12),
get_plot_component(g9a, 'guide-box-bottom', return_all = TRUE),
nrow=2,rel_heights = c(.9,.1))
# Figure 7
g10a <- cumulativeprobplot_ranks(data=cbind(0,wang_posterior),names=treatments)+
theme(legend.position = "bottom") +
guides(color = guide_legend(nrow = 2))
g10b <- cumulativeprobplot_ranks(mcmc=mcmc,names=treatments)
plot_grid(plot_grid(g10a+theme(legend.position = "none"), g10b+theme(legend.position = "none"),
labels = c('A', 'B'), label_size = 12),
get_plot_component(g10a, 'guide-box-bottom', return_all = TRUE),
nrow=2,rel_heights = c(.85,.15) )
# Table 2
res_WANG <- calculate_SUCRA_MNBT(data = cbind(0,wang_posterior),
level = 0.50, names = treatments)
res_RaCENMA <- calculate_SUCRA_MNBT(mcmc = mcmc, level = 0.50, names = treatments)
table2 <- left_join(res_WANG,res_RaCENMA,by="Treatment")[,c(1,2,4,3,5)]
names(table2) <- c("Treatment","SUCRA (Wang)","SUCRA (RaCE)","MNBT (Wang)","MNBT (RaCE)")
print(table2)The following code chunk reproduces the Figures 8–10 (Appendix B).
traceplot_K(mcmc)+
scale_x_continuous(breaks=seq(125000,250000,by=25000),
labels=paste0(seq(125,250,by=25),"k"))
traceplot_mu(mcmc,names=treatments)+
scale_x_continuous(breaks=seq(125000,250000,by=25000),
labels=paste0(seq(125,250,by=25),"k"))
forestplot_muhat(mcmc=mcmc,names=treatments)These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.