Торнадо / двусторонний горизонтальный гистограмма в R с осями диаграммы пересекается с заданным значением (вместо пересечения в нуле) - PullRequest
0 голосов
/ 18 апреля 2019

Я хотел бы построить диаграмму торнадо (двусторонний горизонтальный столбец) в R для детерминированного анализа чувствительности, для которого я пробовал несколько кодов, но не получил желаемого результата.

Я бы хотел достичь: -

  1. График должен быть расположен в порядке убывания чувствительных параметров (т. Е. Самый широкий интервал должен быть представлен в верхней части диаграммы - чтобы получить чувствительность, мы сначала вычисляем абсолютную разницу в значениях нижней границы и значения верхней границы, которые Я назвал "UL_Difference" в моем коде фрейма данных).

  2. Центр не должен быть равен нулю, но должен иметь определенное значение, называемое «базовый случай» или основной / конечный результат моей таблицы результатов (на которой мы хотим проверить влияние различных фиксированных параметров используя нижнюю границу и верхнюю границу значений параметров и генерируя основной результат для нижней грани и верхней граничной величины). Пример кода в Excel VBA:

  3. Сюжет должен иметь заголовок «Сюжет Торнадо для препарата А против Р».

Я пробовал много кодов. Ниже приведен пример, который дает мне сюжет торнадо, но не совсем то, что я хотел получить из R.

Base_Result <- results.table[5,4] # Base/Core result (which I have not used in my codes below yet)

Drug_AP <- seq(1, 48, 4)
D_AP <- data.frame(OWSA[Drug_AP,]) # OWSA[] is a 10x3 matrix with 'Lower_Bound', 'Upper_Bound' and Absolute Difference of the LB and UB termed as 'UL_Difference' (row names are parameters)
DSA_Drug_AP <- D_AP[order(D_AP$UL_Difference, decreasing = T),] # Ordering the data.frame above in Descending order of 'UL_Difference'
cat("DSA Table: Drug A vs P \n")
library(formattable)
print(accounting(as.matrix(DSA_Drug_AP), digits = 0, format = "f", big.mark = ","), right = T) # Just printing the above data.frame

Я попробовал приведенные ниже коды для нанесения торнадо: -

(я не уверен, должен ли я сделать приведенный ниже кадр данных, возможно, это одна из причин, по которой я не получаю желаемый вывод)

dat <- data.frame(Group = c(rep("Lower_Bound", 12), rep("Upper_Bound", 12)), 
                  Parameters = rep(rownames(DSA_Drug_AP), 2), 
                  UL = c(-DSA_Drug_AP[,1], DSA_Drug_AP[,2]))

(Наконец, я построил вышеупомянутый фрейм данных, используя "ggplot", как показано ниже)

library(ggplot2)
ggplot(dat, aes(x = Parameters, y = UL, fill = Group)) + 
    coord_flip() + 
    geom_bar(stat = "identity", position = "identity", width = 0.525) +
    theme(legend.position="top", axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5, size = 10))

И получить вывод, как показано ниже: -

enter image description here

Ниже приведен вывод, который я хотел бы получить (точки 1 и 2 достигнуты; график генерируется из Excel).

enter image description here

# Also, the data I'm using is shown below: -

Base_Result <- 9,504  # Value of results.table[5,4] on which I get 'lower' and 'upper' limit values below (and want tornado with the origin at this base_result).

# My data.frame "D_AP" will look like (I just renamed my parameters to 1(to)12)

           Lower_Bound  Upper_Bound UL_Difference
Parameter_01     8,074      11,181   3,108 
Parameter_02     8,177      11,007   2,831 
Parameter_03     8,879      10,188   1,308 
Parameter_04     4,358      18,697   14,339 
Parameter_05     9,073      10,087   1,013 
Parameter_06     12,034      7,572   4,462 
Parameter_07     11,357      7,933   3,423 
Parameter_08     9,769       9,202   567 
Parameter_09     8,833      10,403   1,570 
Parameter_10     13,450      4,219   9,231 
Parameter_11     10,691      7,915   2,776 
Parameter_12     10,036      8,792   1,244 

# Once, I did sort in descending order then it will be data.frame "DSA_Drug_AP" as below: -

            Lower_Bound Upper_Bound UL_Difference
Parameter_04     4,358      18,697   14,339 
Parameter_10     13,450      4,219   9,231 
Parameter_06     12,034      7,572   4,462 
Parameter_07     11,357      7,933   3,423 
Parameter_01     8,074      11,181   3,108 
Parameter_02     8,177      11,007   2,831 
Parameter_11     10,691      7,915   2,776 
Parameter_09     8,833      10,403   1,570 
Parameter_03     8,879      10,188   1,308 
Parameter_12     10,036      8,792   1,244 
Parameter_05     9,073      10,087   1,013 
Parameter_08     9,769       9,202   567 

# Please note that I need to plot the 1st and 2nd column of values 
# (shown in above table in order of 3rd column as a tornado plot).
# The parameter-## names will come to the left vertical line of plot.

Заранее спасибо!

Ответы [ 2 ]

2 голосов
/ 19 апреля 2019

Я пытался сделать это некоторое время назад, используя geom_bar(), и это было не весело.geom_bar() по умолчанию складывает столбцы с нулем в качестве ссылки.Я должен был создать пустые части в столбцах, чтобы получить (вроде) желаемый эффект.

Лучший способ сделать это - использовать geom_rect().Вам просто нужно немного помассировать фрейм данных, чтобы получить требуемые эстетики xmin, xmax, ymin, ymax (гораздо меньше работы, чем пытаться обойти проблемы с geom_bar())

Поскольку вы не опубликовали свой набор данных, я создал очень простой.Но, надеюсь, структура достаточно близка к вашей


РЕДАКТИРОВАТЬ : я изменил код, чтобы включить фрейм данных в ваш пример.

library(ggplot2)
library(plyr)
library(dplyr)
library(tidyverse)

# this is throwing some warnings in my computer, but it is reading the data frame correctly
df <- '
Parameter Lower_Bound Upper_Bound UL_Difference
Parameter01 8074 11181 3108 
Parameter02 8177 11007 2831 
Parameter03 8879 10188 1308 
Parameter04 4358 18697 14339 
Parameter05 9073 10087 1013 
Parameter06 12034 7572 4462 
Parameter07 11357 7933 3423 
Parameter08 9769 9202 567 
Parameter09 8833 10403 1570 
Parameter10 13450 4219 9231 
Parameter11 10691 7915 2776 
Parameter12 10036 8792 1244
' %>% read_table2()

# original value of output
base.value <- 9504

# get order of parameters according to size of intervals
# (I use this to define the ordering of the factors which I then use to define the positions in the plot)
order.parameters <- df %>% arrange(UL_Difference) %>%
  mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
  select(Parameter) %>% unlist() %>% levels()

# width of columns in plot (value between 0 and 1)
width <- 0.95

# get data frame in shape for ggplot and geom_rect
df.2 <- df %>% 
  # gather columns Lower_Bound and Upper_Bound into a single column using gather
  gather(key='type', value='output.value', Lower_Bound:Upper_Bound) %>%
  # just reordering columns
  select(Parameter, type, output.value, UL_Difference) %>%
  # create the columns for geom_rect
  mutate(Parameter=factor(Parameter, levels=order.parameters),
         ymin=pmin(output.value, base.value),
         ymax=pmax(output.value, base.value),
         xmin=as.numeric(Parameter)-width/2,
         xmax=as.numeric(Parameter)+width/2)

# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
png(width = 960, height = 540)
ggplot() + 
  geom_rect(data = df.2, 
            aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin, fill=type)) +
  theme_bw() + 
  theme(axis.title.y=element_blank(), legend.position = 'bottom',
        legend.title = element_blank()) + 
  geom_hline(yintercept = base.value) +
  scale_x_continuous(breaks = c(1:length(order.parameters)), 
                     labels = order.parameters) +
  coord_flip()
dev.off()

Here is my resulting plot

0 голосов
/ 19 мая 2019

Вот функция для генерации графиков торнадо с использованием ggplot2 :: geom_col () вместе с некоторыми примерами использования. Надеюсь, это поможет ....

# Tornado Plot using ggplot2(), 2019/05/19.
# See Wikipedia: ["Tornado diagram"](https://en.wikipedia.org/wiki/Tornado_diagram).

library( magrittr )
library( tidyverse )

# Function tornado_plot() produces a "tornado plot" given the sensitivity
# analysis results in data_frame df. It plots green bars indicating the levels
# of the response variable when each **x input variable** is moved to its maximum
# level while holding all other variables constant.  Similarly, the red bars are
# the outputs when each **x input variable** is moved to its minimum value while
# holding all other variables constant. The input variable to which the output
# is most sensitive is shown at the top of the plot. And the bars are stacked
# from most sensitive to least sensitive, fancifully yielding the shape of a
# tornado.
tornado_plot <-
  function(
    df,
    var_names_col,
    min_level_col,
    min_output_col,
    max_level_col,
    max_output_col,
    base_level_col,
    baseline_output,
    title_str      = "Tornado Plot",
    subtitle_str   = "",
    caption_str    = "",
    ylab_str       = "output",
    baseline_label = "",
    scale_breaks   = NULL,
    limits         = NULL
  ) {
    # + The argument df must be a tidyverse::tibble with columns referred to by all of the
    #   other arguments having "col" in their names.
    # + The var_names_col argument must be an unquoted column name that contains characters
    #    naming the variables that were varied in the sensitivity analysis.
    # + The level column arguments -- min_level_col, max_level_col and
    #   base_level_col -- must be unquoted column names that contain characters to be
    #   used in forming labels for each variable bar of the plot.
    # + The output column arguments -- min_output_col and max_output_col -- must
    #   be unquoted column names that contain numerical values to be plotted as the
    #   extents of the bars in the plot.
    # + The baseline_output argument is the numeric value of the output (response) variable
    #   produced by setting all of the variables to their base levels.

    var_names_col  <- enquo( var_names_col )
    min_level_col  <- enquo( min_level_col )
    max_level_col  <- enquo( max_level_col )
    base_level_col <- enquo( base_level_col )
    min_output_col <- enquo( min_output_col )
    max_output_col <- enquo( max_output_col )

    have_custom_y_breaks <- !any( is.null(scale_breaks) )

    # Create a generic tibble as the data source for plotting.
    # Sorts variables from the one to which the output was least sensitive
    # to the one to which the output was most sensitive.
    # Then creates labels for each variable capturing the min, base, and max
    # levels of that variable.
    # Finally, it centers all outputs around the baseline output so thta the
    # ggplot2::geom_col() function can still work with zero-based bars.
    plt_df <- df %>% 
      mutate(del = abs(!!max_output_col - !!min_output_col) ) %>% 
      arrange(del) %>% 
      mutate(
        names = sprintf(
          "%s\n(min=%s; base=%s; max=%s)",
          !!var_names_col,
          !!min_level_col,
          !!base_level_col,
          !!max_level_col
        ),
        names = factor(names,names),
        min   = !!min_output_col,
        max   = !!max_output_col
      ) %>% 
      dplyr::select(names,min,max) %>% 
      gather( key = Level, value = output, -names) %>% 
      mutate( output = output - baseline_output, Level = factor(Level,c("min","max")) ) #%T>% print()

    # Generate the tornado plot.
    plt <- plt_df %>% 
      {
        ggplot(., aes( fill = Level, x = names, y = output )) + 
          geom_hline(yintercept = 0, linetype = 1, size = 2, color = "darkgray") +
          geom_col( alpha = 0.4, width = 0.98) + 
          coord_flip() + #*** NOTE THE COORDINATE FLIP ***
          geom_text(aes(y = 0, label = names), size = 4, fontface = "bold" ) +
          scale_x_discrete( expand = expand_scale(add = 1 ) ) +
          scale_fill_manual(values = c(min = "red", max = "green") ) +
          ylab( ylab_str ) +
          theme( # **Hmmm, references the ACTUAL plotted (post-flipped) x-y axes. **
            axis.ticks.y = element_blank(),
            axis.text.y  = element_blank(),
            axis.title.y = element_blank(),
            panel.grid.major.y = element_blank(), # Remove horizontal grid lines
            panel.grid.minor.y = element_blank(),
            axis.text.x  = element_text( size = 14 ),
            axis.title.x = element_text( size = 16 ),
            title        = element_text( size = 18 ),
            legend.position = "bottom"
          ) +
          labs( title = title_str, subtitle = subtitle_str, caption = caption_str )
      }
    # Set the pre-flipped y-axis (which gets flipped to be the x-axis in the final plot).
    if( !is.null(limits) ){
      y_limits = limits
    } else {
      y_limits = c(-max(abs(plt_df$output)),max(abs(plt_df$output)))
    }
    if( have_custom_y_breaks ){
      plt <- plt + scale_y_continuous(
        limits = y_limits,
        breaks = scale_breaks,
        labels = names(scale_breaks)
      )
    } else {
      plt <- plt + scale_y_continuous(
        limits = y_limits,
        labels = function(x) baseline_output + x
      )
    }
    # Add the baseline output label, if any
    if(baseline_label != ""){
      return(
        plt + 
          geom_label(
            data = tibble( x = 0.25, y = 0, label = baseline_label),
            mapping = aes( x = x, y = y, label = label),
            fontface = "bold",
            show.legend = FALSE,
            inherit.aes = FALSE
          )
      )
    } else {
      return( plt )
    }
  }

#--------------------------------------------------------------------------------     

# USAGE EXAMPLE:
# Hypothetical Investment Strategy Analysis:
# These are data from a sensitivity analysis on an investment strategy that invests in an
# an S&P 500 index fund and a "safety" value-store (a 0%-real-return investment); 
# protecting winnings from market with transfer to safety when strategy criteria are met. 
# Disregards taxes and fees. Real values (i.e., inflation-adjusted).
sensitivity_df <- tribble(
  ~variable,                            ~min,  ~base,   ~max, ~Total_at_min, ~Total_base, ~Total_at_max,             ~Time_period,
  "Start Value",                           0,   2000, 100000,        239600,      245900,        554800, "start: 1980, end: 2005",
  "Monthly Investment",                    0,    500,   1000,          6300,      245900,        485600, "start: 1980, end: 2005",
  "Allocation to Safety",                  0,    0.3,    0.5,        277800,      245900,        224700, "start: 1980, end: 2005",
  "Annual Increase in Mo. Investment",     0,   0.01,   0.03,        222700,      245900,        303800, "start: 1980, end: 2005",
  "Protection Rate",                       0, 0.0025,   0.03,        310300,      245900,        199500, "start: 1980, end: 2005",

  "Start Value",                           0,   2000, 100000,        174300,      175900,        253300, "start: 1910, end: 1935",
  "Monthly Investment",                    0,    500,   1000,          1600,      175900,        350100, "start: 1910, end: 1935",
  "Allocation to Safety",                  0,    0.3,    0.5,        177700,      175900,        174600, "start: 1910, end: 1935",
  "Annual Increase in Mo. Investment",     0,   0.01,   0.03,        155600,      175900,        227100, "start: 1910, end: 1935",
  "Protection Rate",                       0, 0.0025,   0.03,        171800,      175900,        176000, "start: 1910, end: 1935"
) %>%  # Add x-input level labels (overwriting reals min, base, max with character values through mutate_at()).
  mutate_at(vars(contains("Total")), ~{100*round(./100)}) %>%
  mutate_at(
    vars( min, base, max), 
    ~ { 
      ifelse(
        abs(.) >= 1000,
        paste0("$",formatC(.,big.mark = ",",format = "f",digits = 0)),
        sprintf(
          c( "$%.0f", "$%.0f", "%.0f%%", "%.1f%%", "%.2f%%" ), 
          . * c(1,1,100,100,100)
        )
      )
    } 
  )

# Generate the tornado plot with generic labeling and axis.
sensitivity_df %>%
  filter( grepl("1980.+2005", Time_period ) ) %>%
  tornado_plot(
    var_names_col   = variable,
    min_level_col   = min,
    min_output_col  = Total_at_min,
    max_level_col   = max,
    max_output_col  = Total_at_max,
    base_level_col  = base,
    baseline_output = .$Total_base[[1]]
  ) %>% print()


# Generate the tornado plot with customized labeling and axis.
scl_limits = c(0, 6.0e5 )
sensitivity_df %>%
  filter( grepl("1980.+2005", Time_period ) ) %>%
  tornado_plot(
    var_names_col   = variable,
    min_level_col   = min,
    min_output_col  = Total_at_min,
    max_level_col   = max,
    max_output_col  = Total_at_max,
    base_level_col  = base,
    baseline_output = .$Total_base[[1]],
    title_str       = "Sensitivity of Total Value to Strategy Variables",
    subtitle_str    = sprintf( "Time period %s", .$Time_period[[1]] ),
    caption_str     = "Assuming S&P 500 index & 0%-real-return 'safe harbor'",
    ylab_str        = "Total Value",
    baseline_label  = paste0("Base Case:\n$",format(100*round(.$Total_base[[1]]/100,0),big.mark = ",")),
    scale_breaks    = setNames(
      seq(min(scl_limits), max(scl_limits), 1e5) - .$Total_base[[1]], 
      paste0("$",formatC(
        seq(min(scl_limits), max(scl_limits), 1e5),big.mark = ",",format = "f",digits = 0)
      )
    ),
    limits          = scl_limits - .$Total_base[[1]]
  ) %>% print()

# Generate the tornado plot for another time period, with scaling
# to be comparable with the first time period.
sensitivity_df %>%
  filter( grepl("1910.+1935", Time_period ) ) %>%
  tornado_plot(
    var_names_col   = variable,
    min_level_col   = min,
    min_output_col  = Total_at_min,
    max_level_col   = max,
    max_output_col  = Total_at_max,
    base_level_col  = base,
    baseline_output = .$Total_base[[1]],
    title_str       = "Sensitivity of Total Value to Strategy Variables",
    subtitle_str    = sprintf( "Time period %s", .$Time_period[[1]] ),
    caption_str     = "Assuming S&P 500 index & 0%-real-return 'safe harbor'",
    ylab_str        = "Total Value",
    baseline_label  = paste0("Base Case:\n$",format(100*round(.$Total_base[[1]]/100,0),big.mark = ",")),
    scale_breaks    = setNames(
      seq(min(scl_limits), max(scl_limits), 1e5) - .$Total_base[[1]], 
      paste0("$",formatC(
        seq(min(scl_limits), max(scl_limits), 1e5),big.mark = ",",format = "f",digits = 0)
      )
    ),
    limits          = scl_limits - .$Total_base[[1]]
  ) %>% print()

Общий сюжет торнадо

Участок торнадо с пользовательскими метками, ось

Сюжет торнадо с пользовательскими метками и в том же масштабе, что и предыдущий

...