R-scripts/Limit_value_exceedance_plot.R at main · PetitPascal/R-scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#------------------------------------------------------------------
#### Configurations ####

## Installing and loading packages
pack_needed<-c("tidyverse","scales","datasets")

is_installed<-pack_needed %in% rownames(installed.packages(all.available=TRUE))
if(any(is_installed == FALSE)){
  install.packages(pack_needed[!is_installed],repos = "http://cran.us.r-project.org")
}
invisible(lapply(pack_needed, library, character.only = TRUE))

#------------------------------------------------------------------
#### Data preparation ####

## Using the airquality dataset from the datasets package as example
test<-as_tibble(na.omit(datasets::airquality))

## Transforming the month variable into a factor
test$Month<-factor(test$Month,levels=sort(unique(test$Month)),labels=c("May","June","July","August","September"))

## Calculating summary statistics
test<-test %>% group_by(Month) %>% summarize(median=median(Ozone),      # median
                                             P5=quantile(Ozone,0.05),   # 5th percentile
                                             P95=quantile(Ozone, 0.95), # 95th percentile
                                             P99=quantile(Ozone, 0.99)) # 99th percentile

#------------------------------------------------------------------
#### Limit value exceedance assessment ####

## Setting a limit value for ozone concentration
limit_value<-75

## Determining if the limit value is exceeded
test<-test %>% mutate(comp=case_when(P95>limit_value~"Yes",
                                     P95<=limit_value&P99>limit_value~"Inconclusive",
                                     P99<=limit_value~"No",
                                     T~"No"))

## Transforming the limit value exceedance variable into a factor
test$comp<-factor(test$comp,levels=unique(test$comp))

#------------------------------------------------------------------
#### Plot ####

test %>% ggplot() +
  geom_point(aes(x=Month,y=median),size=2) +
  geom_errorbar(aes(x=Month,ymin=P95,ymax=P99))+
  geom_crossbar(aes(x=Month, y=median,ymin = P5, ymax = P95,fill=comp),width=0.2)+
  scale_x_discrete("Month") +
  scale_y_log10("Ozone concentration",breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x)))+
  scale_fill_manual("Exceedance of the limit value?",
                    values=c("Yes"="red","No"="#009E73","Inconclusive"="gold"))+
  geom_hline(yintercept=limit_value,linetype="dashed",col="navyblue")+
  annotate(geom="text",x=0.5,y=limit_value*1.1,label='limit value',col="navyblue")+
  theme_bw()+
  theme(strip.text.x = element_text(size = 11, face="bold",colour = "black", angle = 0,margin = margin(0.3,0,0.3,0, "cm")),
        strip.background = element_rect(fill=alpha('#009E73', 0.2),
                                        colour="black", size=1),
        axis.text = element_text(size=12,color="black"),
        axis.title=element_text(size=12,face="bold",color="black"),
        legend.text = element_text(size = 12, face = "bold"),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        axis.line = element_line(color = "black",size = 0.1, linetype =
                                   "solid"))+
  coord_flip()