R/ForestPlot.R

In lsyang99/changhua.sensor: Functions specifically designed for analyzing the "ens" data.

#' A Self-made Function for forest plot.
#'
#' This is an extended function of the function "SummaryEns". It uses the summary sheet from "SummaryEns" to do forest plot. The function setwd() is recommended to use beforhand to specify the path to the preferred folder.
#' @param range The range of the Y axis. The default is from -5 to 5.
#' @param by The increment of the sequence. The default is 0.5.
#' @param summary_path The path of the xlsx file which is the output of the function called SummaryEns.
#' @param subtitle Additional subtitle of the plot.
#' @param type The function can build either a forest plot or a box plot. The default "f" will lead to a forest plot, while "b" will give a box plot.

ForestPlot <- function(summary_path , range = c(-5 , 5) , by = 0.5 , subtitle = NULL , type = 'f'){ # 建議先 setwd() 以建立輸出之圖檔路徑
  require(data.table) ; require(readxl) ; require(magrittr) ; require(ggplot2) ; require(tidyr) ; require(ggpubr)
  sheets <- excel_sheets(path = summary_path)
  summary <- # 辨識excel檔中全部的sheet名稱，逐一呼叫後儲存成一list
    lapply(sheets , function(x){
      read_excel(summary_path , sheet = x)
    })

  # between variability
  # data fromat for multiple x scales
  X <- lapply(1:length(summary) , function(i){
    x1 <- cbind.data.frame(
      value = summary[[i]][7 , - c(1:2)] %>% t ,
      type = 'RSQ' ,
      period = sheets[i]
    )

    x2 <- cbind.data.frame(
      value = summary[[i]][8 , - c(1:2)] %>% t ,
      type = 'slope' ,
      period = sheets[i]
    )

    x3 <- cbind.data.frame(
      value = summary[[i]][9 , - c(1:2)] %>% t ,
      type = 'intercept' ,
      period = sheets[i]
    )

    Reduce(
      function(a , b) merge(a , b , all = T) ,
      list(x1 , x2 , x3)
    )
  })

  B.table <- Reduce(
    function(a , b) merge(a , b , all = T) ,
    X) %>%
    group_by(type , period) %>%  # 新增上界與下界範圍
    mutate(
      L.limit = min(value) ,
      U.limit = max(value)
    )

  # within variability
  XX <- suppressWarnings(
    lapply(1:length(summary) , function(i){
      raw <- summary[[i]][-c(1:9) , ]
      raw1 <- raw[seq(1 , 42 , 3) , ] %>% unlist %>% as.numeric  # RSQ
      RSQ <- raw1[is.na(raw1) == FALSE & raw1 != 1]  # 去除NA與1
      N <- 91 - length(RSQ)
      RSQ <- c(RSQ , rep(1 , N))

      raw2 <- raw[seq(2 , 42 , 3) , ] %>% unlist %>% as.numeric  # slope
      slope <- raw2[is.na(raw2) == FALSE & raw2 != 1]  # 去除NA與1
      N <- 91 - length(slope)
      slope <- c(slope , rep(1 , N))  # 如有例外出現1的，會不小心去除，所以最後要補回來

      raw3 <- raw[seq(3 , 42 , 3) , ] %>% unlist %>% as.numeric  # intercept
      intercept <- raw3[is.na(raw3) == FALSE & raw3 != 0]  # 去除NA與0
      N <- 91 - length(intercept)
      intercept <- c(intercept , rep(1 , N))  # 如有例外出現0的，會不小心去除，所以最後要補回來

      cbind.data.frame(RSQ , slope , intercept) %>%
        gather(. , key = 'type' , value = 'value') %>%
        mutate(. , period = sheets[i])
    })
  )
  W.table <- Reduce(
    function(a , b) merge(a , b , all = T) ,
    XX) %>%
    group_by(type , period) %>%  # 新增上界與下界範圍
    mutate(
      L.limit = min(value) ,
      U.limit = max(value)
    )



  # ggplot
  if(type == 'f'){
    plot1 <- ggplot(B.table , aes(x = period , y = value , ymin = L.limit , ymax = U.limit)) +
      geom_errorbar(aes(ymin = L.limit , ymax = U.limit)  , width = 0.25) +  # forest plot
      labs(title = expression('Between Variability (EPA ' * PM[2.5] * ' vs. Sensor ' * PM[2.5] * ')') ,
           subtitle = subtitle , y = 'Linaer Regression Parameters and Their Ranges') +
      scale_y_continuous(breaks = seq(range[1] , range[2] , by)) +
      coord_cartesian(ylim = c(range[1] , range[2]) , expand = FALSE , clip = 'off') +
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'point' , size = 4.5) +  # 標上平均值
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'text' , size = 6 , vjust = -1 , aes( label = round(..y.. , digits = 3 )) , family = 'Times') +
      theme_bw() +
      theme(text = element_text(family = 'Times' , size = 21) ,
            axis.title = element_text(size = 21) ,
            plot.title = element_text(size = 25 , hjust = 0.5 , face = 'bold') ,
            plot.subtitle = element_text(size = 21 , hjust = 0.5)) +
      facet_wrap( ~ type , scales = 'fixed' , ncol = 3)

    plot2 <- ggplot(W.table , aes(x = period , y = value , ymin = L.limit , ymax = U.limit)) +
      geom_errorbar(aes(ymin = L.limit , ymax = U.limit)  , width = 0.25) +  # forest plot
      labs(title = expression('Within Variability (EPA ' * PM[2.5] * ' vs. Sensor ' * PM[2.5] * ')') ,
           subtitle = subtitle , y = 'Linaer Regression Parameters and Their Ranges') +
      scale_y_continuous(breaks = seq(range[1] , range[2] , by)) +
      coord_cartesian(ylim = c(range[1] , range[2]) , expand = FALSE , clip = 'off') +
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'point' , size = 4.5) +  # 標上平均值
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'text' , size = 6 , vjust = -1 , aes( label = round(..y.. , digits = 3 )) , family = 'Times') +
      theme_bw() +
      theme(text = element_text(family = 'Times' , size = 21) ,
            axis.title = element_text(size = 21) ,
            plot.title = element_text(size = 25 , hjust = 0.5 , face = 'bold') ,
            plot.subtitle = element_text(size = 21 , hjust = 0.5)) +
      facet_wrap( ~ type , scales = 'fixed' , ncol = 3)
  }else{
    plot1 <- ggplot(B.table , aes(x = period , y = value)) +  # lex.order: x 軸排列按照輸入順序
      geom_boxplot(width = 0.5 , outlier.shape = 3 , outlier.size = 1.8) +
      labs(title = expression('Between Variability (EPA ' * PM[2.5] * ' vs. Sensor ' * PM[2.5] * ')') ,
           subtitle = subtitle , y = 'Linaer Regression Parameters and Their Ranges') +
      scale_y_continuous(breaks = seq(range[1] , range[2] , by)) +
      coord_cartesian(ylim = c(range[1] , range[2]) , expand = FALSE , clip = 'off') +

      stat_summary(fun = mean, colour = 'darkblue' , geom = 'point' , size = 4.5) +
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'text' , size = 6 , vjust = -1 , aes( label = round(..y.. , digits = 2 )) , family = 'Times') +
      theme_bw() +
      theme(text = element_text(family = 'Times' , size = 21) ,
            axis.title = element_text(size = 21) ,
            plot.title = element_text(size = 25 , hjust = 0.5 , face = 'bold') ,
            plot.subtitle = element_text(size = 21 , hjust = 0.5)) +
      facet_wrap( ~ type , scales = 'fixed' , ncol = 3)


    plot2 <- ggplot(W.table , aes(x = period , y = value)) +  # lex.order: x 軸排列按照輸入順序
      geom_boxplot(width = 0.5 , outlier.shape = 3 , outlier.size = 1.8) +
      labs(title = expression('Within Variability (EPA ' * PM[2.5] * ' vs. Sensor ' * PM[2.5] * ')') ,
           subtitle = subtitle , y = 'Linaer Regression Parameters and Their Ranges') +
      scale_y_continuous(breaks = seq(range[1] , range[2] , by)) +
      coord_cartesian(ylim = c(range[1] , range[2]) , expand = FALSE , clip = 'off') +

      stat_summary(fun = mean, colour = 'darkblue' , geom = 'point' , size = 4.5) +
      stat_summary(fun = mean, colour = 'darkblue' , geom = 'text' , size = 6 , vjust = -1 , aes( label = round(..y.. , digits = 2 )) , family = 'Times') +
      theme_bw() +
      theme(text = element_text(family = 'Times' , size = 21) ,
            axis.title = element_text(size = 21) ,
            plot.title = element_text(size = 25 , hjust = 0.5 , face = 'bold') ,
            plot.subtitle = element_text(size = 21 , hjust = 0.5)) +
      facet_wrap( ~ type , scales = 'fixed' , ncol = 3)

  }
  ggsave(plot1 , file = 'between.png' , width = 40 , height = 20 , units = 'cm')
  ggsave(plot2 , file = 'within.png' , width = 40 , height = 20 , units = 'cm')

  return( ggarrange(plot1 , plot2 , ncol = 1) )

}