doseresNMA antidep/drnma.plot.2absolute.R

In htx-r/doseresNMA:

# Plot: 
# y: absolute response: 1. model 1, 2. model2 and 3. placebo,
# x=dose 
# + rug of obseravtions
require(dplyr)
require(ggplot2)
absolutePredplot2 <- function(m=m,
                             m2=m2,
                             ref.lab='placebo',
                             data=antidep, ...){
  
  #-- Prepare data to plot
  # reference numeric index
  drug.lab=levels(data$drug)
  nd <- 200
  ref.index <- which(levels(as.factor(drug.lab))==ref.lab)
  
  # 1. y-axis: p.drug M1
  plotdata.drug <- m$BUGSoutput$summary %>% 
    t() %>% 
    data.frame() %>% 
    select(starts_with('p.drug'))%>%
    t()%>% 
    data.frame() 
  # 2. y-axis: p.drug M2
  plotdata2 <- m2$BUGSoutput$summary %>% 
    t() %>% 
    data.frame() %>% 
    select(starts_with('p.drug'))%>%
    t()%>% 
    data.frame()
  
  # 3. add placebo data
  data.placebo <- m$BUGSoutput$summary%>% 
    t() %>% 
    data.frame() %>% 
    select(starts_with('p.placebo'))%>%
    t()%>% 
    data.frame() 
  plotdata <- plotdata.drug%>%
    mutate(X50.placebo=data.placebo$X50.,
           X2.5.placebo=data.placebo$X2.5.,
           X97.5.placebo=data.placebo$X97.5.)
  # 4. x-axis: doses per drug
  ndrugs <- length(drug.lab)-1 # number of drugs
  max.dose <- unlist(data%>%
                       group_by(drug)%>%
                       group_map(~round(max(.x$dose,na.rm = T))))[-ref.index]
  plotdata$dose <- unlist(
    sapply(1:ndrugs, 
           function(i) seq(0,max.dose[i],l=nd),
           simplify = FALSE
    )
  )
  
  
  # 5. labels: drug name
  plotdata$drug <- rep(drug.lab[drug.lab!=ref.lab],each=nd)
  
  # 6. rug: observed doses
  data0 <- data %>% # add a placebo indicator column
    group_by(studyid)%>%
    group_map(~mutate(.x,
                      is.placebo=ifelse(sum(.x$drug=='placebo')>0,0,NA)
    )
    ) # add 0 if we have a trial compare drug with placebo OR NA if not
  data0 <- data.frame(do.call(rbind,data0))
  
  plotdata$obs.dose <- unlist(sapply(1:ndrugs, 
                                     function(i){
                                       dose_plc_drug <- unique(unlist( # a vector of placebo(0/NA) and drug doses
                                         data0[data0$drug%in%drug.lab[-ref.index][i],][c('is.placebo','dose')]
                                       ))
                                       
                                       dose_plc_drug_nd <-rep(dose_plc_drug[!is.na(dose_plc_drug)],
                                                              l=nd) # repeat the vector until nd length to fit the data 
                                       dose_plc_drug_nd
                                     }
                                     ,simplify = F
  )
  )
  
  #-- To plot     
  # y axis limits
  ymax <- round(max(plotdata$X97.5.),1)
  ymin <- round(min(plotdata$X2.5.),1)
  
  # theme
  theme_set(
    theme_minimal() +
      theme(legend.position = "right",
            panel.background = element_rect(fill = 'snow1',colour = 'white'),
            axis.text.x = element_text(size=10),axis.text.y = element_text(size=10),
            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))
  )
  
  # plot
  g_X50 <- # the median line of drug and placebo
    ggplot(plotdata,aes(x=dose)) +       
    geom_line(aes(y=X50.,color='treatment M1'))+       # p.drug median
    geom_line(aes(y=plotdata2[,'X50.'],color='treatment M3'))+
    geom_line(aes(y=X50.placebo,color='placebo'))+  # p.placebo: median
    scale_color_manual(values=c('steelblue','coral4','darkgreen'),name="")+ 
    guides(color = guide_legend(override.aes = list(size = 1.5)))
  
  g_crI <-  # add 95% CrI of p.drug and p.placebo
    g_X50+geom_smooth(                                    
      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=X50.placebo, ymin=X2.5.placebo,ymax=X97.5.placebo),
      color='steelblue',fill='steelblue2',
      data=plotdata, stat="identity")+
    coord_cartesian(ylim = c(ymin, ymax))#+
  g_rug <- # add rug of the observed doses
    g_crI+geom_rug(#data=plotdata,
      mapping=aes(x=obs.dose),
      inherit.aes = F,
      color='red',
      length = unit(0.05, "npc"))
  
  g_split <- g_rug + facet_wrap(~drug,scales = 'free_x')
  
  g_labs <- g_split + ggplot2::labs(y="Predicted absolute response", x="Dose")
  
  g <- g_labs + ggplot2::scale_linetype_manual(name="")
  
  g
}