doseresNMA antidep/drnma.plot.forest.R

In htx-r/doseresNMA:

# forest plot: dose-response coefficients
require(ggplot2)
require(dplyr) 
forestplot <- function(m=m, 
                       ref.lab='placebo',
                       drug.lab=levels(antidep$drug),
                       metareg=FALSE,...){
  #-- data to plot --
  param.lab <-c('beta1','beta2')
  #param.lab.exp <-c(expression(beta[1]^2),expression(beta[2]^2)) # greek beta's
    if(metareg){
      param.lab0 <- 'g.'
      param.lab <- append(param.lab,param.lab0)
    }
    ndrugs <- length(drug.lab)
    
    plotdata <- m$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,]

    #-- plot --
    
    # theme
  theme_set(
    theme_minimal() +
      theme(legend.position = "right",
            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"),
            strip.background =element_rect(fill="snow3"),
            strip.text.x = element_text(size = 14)
            )
  )
  
  g_crI <- ggplot(plotdata, 
              aes(y=plotdata$X50., x=plotdata$drug)) +
    geom_point()+ # point estimate
    geom_errorbar(aes(ymin=plotdata$X2.5., ymax=plotdata$X97.5.), # 95%CrI
                  width=0.4,size=0.5) +
    coord_flip() 
  
  g_split <- g_crI + facet_wrap(~doseparam, scales="free_x")
  
  g <- g_split + 
    xlab("Drug") +
    ylab("logOR")
  g 
  
}