# I pick "Outfall.Numbers" , "Number.of.Overflow.Events" and "Number.of.Permitted.Outfalls" as key variables.
library(expss)
data(nysdata)
cro_cpct(nysdata$Number.of.Overflow.Events, list(total(), nysdata$Number.of.Permitted.Outfalls, nysdata$Outfall.Number))
 #Total   1   2   3   4   5   6   8   9   10   11   12   13   16   22   32   33   36   44   47   49   52   62   79   001A   002A   002B   003A   003B   005A   005B   005E   006A   007A   007E   008A   008B   013A   013B   017E   023A   023B   029A   029A-E   030A   046A   046B   068 (PSO)   069 (PSO)   070 (PSO)   083 (PSO)   1   10   100   108   11   110   116   117   118   119   12   124   125   127   13   139   14   140   142   143   15   16   17   18   19   2   20   21   22   23   24   25   26   27   28   29   3   30   31   32   33   34   35   36   37   38   39   4   40   41   42   43   44   45   46   47   48   49   5   50   51   52   53   54   55   56   57   58   59   6   60   61   62   63   64   65   66   67   68   69   7   70   71   72   73   74   75   76   77   78   79   8   80   81   82   83   85   87   88   9   T001   T002   T003 
 nysdata$Number.of.Overflow.Events 
   0  5.6 2.9 5.4 22.2 23.5 8 20.3 6.4 19.2 5.6 10.5 13.3 25 6.7 13.3 8.3 9.1 12.5 8.3 7.7 18.2 3.1 10 14.3 20 6.9 25 33.3 12.9 20 33.3 33.3 33.3 7.1 3.6 5.9 100 10
   1  1.3 10 5.9 5.1 6.4 1.9 8.7 5.3 100 8.3 5.4 33.3 3.6 33.3
   1.2  0.1 20 3.2
   1.4  0.1 1.7 6.7
   10  1.2 12.5 2.7 2.9 4.2 1.7 2.1 5.8 25 8.3 2.7 11.1 3.4 25 25 33.3 3.6
   11  0.9 2.7 5.1 4.3 1.9 33.3 8.3 3.1 11.1 3.4 3.2 5.9
   12  1.7 12.5 16.7 2.9 2.7 11.8 4 1.7 2.1 2.1 1.9 100 8.3 6.7 5.4 3.4 25 33.3 3.2 3.6 7.1 5
   13  0.9 2.9 5.9 4.2 2.1 3.8 5.3 100 8.3 11.1 3.4 33.3 5.9
   136  0.1 2.9 3.1
   14  0.4 4.3 1.9 100 33.3 20
   15  0.7 12.5 5.6 6.4 6.7 10 3.1 3.2 3.6
   16  0.5 12.5 20 2.9 2.1 8.7 9.1 3.4
   17  0.3 2.7 2.1 5.4
   18  0.4 4.3 2.1 7.7 8.3 3.4
   19  0.3 10 4.2 6.7 2.7
   2  0.8 12.5 2.7 4.2 6.4 4.3 6.7 100 8.3 2.7 100
   20  0.1 2.1 25
   21  0.7 1.7 2.1 6.2 8.3 12.5 16.7 25 3.6
   22  0.3 1.7 2.1 10 3.6
   23  0.7 16.7 5.9 2.9 2.1 2.1 100 6.2 2.7 25 5
   24  0.7 20 5.6 2.9 2.1 2.1 100 2.7 3.4 25 3.6
   25  0.1 2.1 6.7
   26  0.7 8.1 2.9 2.1 100 2.7 3.1 3.6 3.6
   264  0.1 2.9 3.6
   27  0.3 5.9 4.2 8.3 2.7
   28  0.4 16.7 5.6 2.1 2.7 9.1 3.4
   29  0.3 12.5 2.1 2.7 20
   3  0.9 11.8 1.7 4.3 100 6.7 6.7 8.3 8.3 3.1 3.2
   30  0.3 5.6 2.1 3.1 20
   31  0.4 16.7 4 2.1 8.3 3.1 3.4
   32  0.3 4 1.7 3.6 5
   33  0.5 5.6 8 2.1 4.3 6.2 2.7 7.7
   34  0.7 16.7 8 2.1 1.9 5.6 6.7 16.7 3.1 20
   35  0.1 1.7 3.1
   36  0.1 2.9 5.3
   37  0.5 4.2 4 4.2 20 20 5
   38  0.1 2.9 6.7
   39  0.4 4.2 4 1.9 7.7 3.1 5.9
   4  1.5 5.4 2.9 10.6 5.8 6.7 3.4 3.2 33.3 25 25 20 25 33.3 5 100
   40  0.8 5.9 1.7 2.1 6.2 100 100 6.7 6.7 3.4 3.6
   41  0.5 16.7 2.1 3.8 5.3 2.7 25 25
   42  0.4 5.9 4 1.9 5.3 6.7 3.6
   43  0.4 5.6 2.9 4 3.1 3.2 3.6
   44  0.7 20 4.2 3.8 5.6 5.3 8.3 3.1 25
   45  0.5 4.2 1.7 4.2 6.7 2.7 11.1 25
   46  0.1 2.1 3.4
   47  0.1 2.1 6.7
   48  0.1 2.1 6.7
   49  0.3 3.8 6.7 8.3
   5  0.9 20 4 3.4 5.8 4.3 5.6 5.3 3.4 25 25 3.6
   50  0.4 6.2 16.7 11.1 5
   51  0.5 4.2 6.2 100 9.1 12.5 3.6
   52  0.4 5.9 4.2 5.6 12.5 5.9
   53  0.4 4.2 4.2 10 3.1 14.3
   55  0.1 2.1 3.2
   56  0.4 4.2 1.9 14.3 33.3 3.6
   57  0.3 4.2 5.6 8.3
   58  0.1 4.2 6.2
   59  0.1 4 3.2
   6  1.1 10 5.9 8.5 2.1 1.9 100 100 8.3 2.7 3.1 12.5 25 5.9
   61  0.1 1.7 4.3
   62  0.1 2.1 10
   63  0.1 1.9 10
   65  0.4 2.1 3.8 100 20 3.6
   68  0.1 2.7 2.7
   69  0.1 1.9 11.1
   7  0.7 5.9 2.1 3.8 6.2 6.7 3.1 25 3.6
   78  0.1 5.6 3.6
   8  0.4 2.7 2.9 2.1 50 5.6 100
   88  0.1 2.9 5.9
   9  0.5 4 4.3 1.9 11.1 20 25 3.6
   < 6 annually  1.2 52.9 5.6 5.3 6.7 2.7 3.1 3.2 3.6 3.6 5
   Average of <4 overflows annually  1.2 36 100 100 100 100 100 100 4.3 3.6 5
   Eliminated  0.3 4.2 2.1 5.3 100
   No Data  11.8 50 25 55.9 45.9 11.8 38.9 2.9 50 18.6 8.5 17.3 50 33.3 50 100 100 100 100 100 8.7 16.7 100 5.3 100 6.2 20 8.3 6.7 13.3 16.7 8.3 21.6 9.1 8.3 30 7.7 16.7 10 11.1 9.1 12.5 12.5 12.5 10 12.5 20 17.2 16.7 20 19.4 14.3 50 33.3 17.9 33.3 50 17.6 15 100
   Visit facility website  39.7 20 29.4 16.2 100 28.8 100 100 100 100 100 100 25 33.3 100 100 100 52.2 33.3 31.6 33.3 37.5 46.7 50 46.7 46.7 58.3 50 21.6 54.5 75 50 60 53.8 58.3 50 50 55.6 54.5 31.2 50 62.5 50 66.7 62.5 57.1 71.4 55.6 40 40 27.6 50 50 40 40 25 50 33.3 25 25 33.3 32.3 25 25 40 40 50 50 33.3 33.3 33.3 32.1 25 25 50 33.3 33.3 33.3 33.3 33.3 100 100 28.6 100 100 100 50 50 50 33.3 50 50 35.3 50 100 100 100 100 100 30
   http://cso.savetherain.us/ ; 8.8 100 100 100 100 100 100 5.6 5.3 6.2 6.7 8.3 6.7 6.7 8.3 9.1 12.5 8.3 7.7 16.7 11.1 9.1 3.1 12.5 12.5 10 11.1 12.5 14.3 14.3 11.1 20 20 3.4 16.7 25 20 20 25 25 33.3 25 25 3.2 33.3 25 25 20 20 33.3 33.3 33.3 25 25 33.3 33.3 33.3 33.3 3.6 50 50 50 33.3 50 100 5.9 50 5
   #Total cases  754 10 8 6 34 5 37 17 18 11 34 24 25 59 22 34 35 37 43 47 48 52 66 82 4 1 1 3 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 23 18 1 1 19 1 1 1 1 1 15 1 1 1 16 1 15 1 1 1 12 15 15 12 12 37 11 8 12 10 13 12 10 6 9 11 32 8 8 10 9 8 7 7 9 5 5 29 6 4 5 5 4 4 3 4 4 3 31 3 4 4 5 5 4 4 3 3 3 28 4 4 2 3 3 3 3 3 1 1 28 1 1 1 2 2 2 3 2 2 1 17 2 1 1 1 1 1 1 20 1 1 1 
## We can view the relation between those 3 variables from the table.
#Simple Example
nysdata %>% 
    tab_cells(Number.of.Overflow.Events) %>% 
    tab_cols(total(), Number.of.Permitted.Outfalls) %>% 
    tab_stat_cpct() %>% 
    tab_pivot()
 #Total     Number.of.Permitted.Outfalls 
   1   2   3   4   5   6   8   9   10   11   12   13   16   22   32   33   36   44   47   49   52   62   79 
 Number.of.Overflow.Events 
   0  5.6   2.9 5.4 22.2 23.5 8 20.3 6.4 19.2
   1  1.3   10 5.9 5.1 6.4 1.9
   1.2  0.1   20
   1.4  0.1   1.7
   10  1.2   12.5 2.7 2.9 4.2 1.7 2.1 5.8
   11  0.9   2.7 5.1 4.3 1.9
   12  1.7   12.5 16.7 2.9 2.7 11.8 4 1.7 2.1 2.1 1.9
   13  0.9   2.9 5.9 4.2 2.1 3.8
   136  0.1   2.9
   14  0.4   4.3 1.9
   15  0.7   12.5 5.6 6.4
   16  0.5   12.5 20 2.9 2.1
   17  0.3   2.7 2.1
   18  0.4   4.3 2.1
   19  0.3   10 4.2
   2  0.8   12.5 2.7 4.2 6.4
   20  0.1   2.1
   21  0.7   1.7 2.1 6.2
   22  0.3   1.7 2.1
   23  0.7   16.7 5.9 2.9 2.1 2.1
   24  0.7   20 5.6 2.9 2.1 2.1
   25  0.1   2.1
   26  0.7   8.1 2.9 2.1
   264  0.1   2.9
   27  0.3   5.9 4.2
   28  0.4   16.7 5.6 2.1
   29  0.3   12.5 2.1
   3  0.9   11.8 1.7 4.3
   30  0.3   5.6 2.1
   31  0.4   16.7 4 2.1
   32  0.3   4 1.7
   33  0.5   5.6 8 2.1
   34  0.7   16.7 8 2.1 1.9
   35  0.1   1.7
   36  0.1   2.9
   37  0.5   4.2 4 4.2
   38  0.1   2.9
   39  0.4   4.2 4 1.9
   4  1.5   5.4 2.9 10.6 5.8
   40  0.8   5.9 1.7 2.1 6.2
   41  0.5   16.7 2.1 3.8
   42  0.4   5.9 4 1.9
   43  0.4   5.6 2.9 4
   44  0.7   20 4.2 3.8
   45  0.5   4.2 1.7 4.2
   46  0.1   2.1
   47  0.1   2.1
   48  0.1   2.1
   49  0.3   3.8
   5  0.9   20 4 3.4 5.8
   50  0.4   6.2
   51  0.5   4.2 6.2
   52  0.4   5.9 4.2
   53  0.4   4.2 4.2
   55  0.1   2.1
   56  0.4   4.2 1.9
   57  0.3   4.2
   58  0.1   4.2
   59  0.1   4
   6  1.1   10 5.9 8.5 2.1 1.9
   61  0.1   1.7
   62  0.1   2.1
   63  0.1   1.9
   65  0.4   2.1 3.8
   68  0.1   2.7
   69  0.1   1.9
   7  0.7   5.9 2.1 3.8
   78  0.1   5.6
   8  0.4   2.7 2.9 2.1
   88  0.1   2.9
   9  0.5   4 4.3 1.9
   < 6 annually  1.2   52.9
   Average of <4 overflows annually  1.2   36
   Eliminated  0.3   4.2 2.1
   No Data  11.8   50 25.0 55.9 45.9 11.8 38.9 2.9 50.0 18.6 8.5 17.3
   Visit facility website  39.7   20 29.4 16.2 100 28.8 100 100 100 100 100 100
   http://cso.savetherain.us/ ; 8.8   100
   #Total cases  754   10 8 6 34 5 37 17 18 11 34 24 25 59 22 34 35 37 43 47 48 52 66 82 
#Table with Caption
nysdata %>% 
    tab_cells(DEC.Region) %>%
    tab_cols(total(), Number.of.Permitted.Outfalls) %>% 
    tab_stat_mean_sd_n() %>%
    tab_last_sig_means(subtable_marks = "both") %>% 
    tab_pivot() %>% 
    set_caption("Table with summary statistics and significance marks")
Table with summary statistics and significance marks
 #Total     Number.of.Permitted.Outfalls 
   1     2     3     4     5     6     8     9     10     11     12     13     16     22     32     33     36     44     47     49     52     62     79 
   A     B     C     D     E     F     G     H     I     J     K     L     M     N     O     P     Q     R     S     T     U     V     W 
 DEC.Region 
   Mean  4.4    6.2 > D E G I N O P Q R T W < U   5.9 > D E G I M N O P Q R T W < H U V   5.0 > D I N O P Q R W < H K U V   2.9 > I N O P Q R W < A B C E F G H J K L M S T U V   4.0 > D < A B F H J K   5.5 > D E G I M N O P Q R T W < H U V   4.0 > D < A B F H J K   8.0 > B C D E F G I J K L M N O P Q R S T V W < U   2.0 < A B C D F H J K L M   6.1 > D E G I L M N O P Q R T W < H U V   6.7 > C D E G I L M N O P Q R T W < H U   4.8 > D I N O P Q R W < H J K S U V   4.4 > D I N O P Q R W < B F H J K S U V   2.0 < A B C D F H J K L M   2.0 < A B C D F H J K L M   2.0 < A B C D F H J K L M   2.0 < A B C D F H J K L M   2.0 < A B C D F H J K L M   6.0 > D L M < H   4.0 > D < A B F H J K   9.0 > A B C D F H J K L M   7.0 > B C D F J L M < H   2.0 < A B C D F H J K L M
   Std. dev.  2.5    3.0     0.8     1.1     0.8     0.0     2.6     0.0     1.0     0.0     1.5     2.5     2.4     1.9     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0     0.0  
   Unw. valid N  754.0    10.0     8.0     6.0     34.0     5.0     37.0     17.0     18.0     11.0     34.0     24.0     25.0     59.0     22.0     34.0     35.0     37.0     43.0     47.0     48.0     52.0     66.0     82.0   
# ggplot graph
library(tidyverse)
ggplot(nysdata, aes(County, Number.of.Permitted.Outfalls)) +
stat_summary(fun.y="sum",geom="bar",alpha=.4,fill="cornflowerblue")+
theme_minimal()+
coord_flip()+
ggtitle("Permitted Outfalls by County")+
xlab("County")+ylab("Number of Permitted Outfalls")

#Barplot
ggplot(nysdata, aes(x=Number.of.Permitted.Outfalls, y=1))+
geom_bar(stat='identity')

#Boxplot
ggplot(nysdata, aes(Number.of.Permitted.Outfalls, Number.of.Overflow.Events, label=rownames(nysdata))) +
geom_boxplot(fill = "grey80", colour = "#3366FF") +
ggtitle("Figure 1: Miles per Gallon by Cylinders")

# Include at least one meaningful annotation on the graphs to highlight a key feature of the results
p = ggplot(nysdata, aes(x = Number.of.Overflow.Events, y = Number.of.Permitted.Outfalls)) + geom_point()
p + annotate("rect", xmin = 70, xmax = 80, ymin = 60, ymax = 80,alpha = .2)

## Shaded area represents outliers.