data-viz-proj.Rmd

---
title: "Data Vizualization Project"
author: Andrii Yurkiv
---

```{r}
library(ggplot2)
library(readxl)
library(tidyr)
library(dplyr)
library(maps)
library(maptools)
library(data.table)
library(ggpubr)
library(ggrepel)
library(ellipse)
```

```{r}
garbadge_patch <- read_excel('papers/data/garbadge_patch/Lebreton2018_HistoricalDataset.xls')
```

```{r}
garbadge_patch$`Sampling Month` <- as.integer(garbadge_patch$`Sampling Month`)

garbadge_patch$date <- as.Date(paste(1, garbadge_patch$`Sampling Month`, garbadge_patch$`Sampling Year`, sep = '-'), format = '%d-%m-%Y')
```

```{r}
names(garbadge_patch) <- c('id', 'data_origin', 'data_reference', 'year', 'month', 'lon', 'lat', 'sampling_type', 'mesh_size_microns', 'net_tow_depth', 'microplastic_numerical_concentration_n_km2', 'mictoplastic_mass_concentration_g_km2', 'inside_GPGP', 'date')
```

```{r}
garbadge_patch$data_origin <- NULL
garbadge_patch$data_reference <- NULL
garbadge_patch$year <- NULL
garbadge_patch$month <- NULL
garbadge_patch$sampling_type <- NULL
garbadge_patch$mesh_size_microns <- NULL
garbadge_patch$net_tow_depth <- NULL
```

```{r}
min(garbadge_patch$lat)
max(garbadge_patch$lat)

min(garbadge_patch$lon)
max(garbadge_patch$lon)
```

```{r}
garbadge_patch$microplastic_numerical_concentration_n_km2 <- as.integer(garbadge_patch$microplastic_numerical_concentration_n_km2)
```



```{r}
garbadge_patch <- garbadge_patch %>% mutate_if(is.numeric, ~replace(., is.na(.), 0)) %>% mutate_if(is.integer, ~replace(., is.na(.), 0))
```

```{r}
garbadge_patch
```



```{r}
grouped_by_date_patch <- garbadge_patch %>%
  group_by(date, inside_GPGP) %>% 
  summarize(mmnc = mean(microplastic_numerical_concentration_n_km2, na.rm = TRUE), mmmc = mean(mictoplastic_mass_concentration_g_km2, na.rm = TRUE))
```




```{r}
grouped_by_date_patch$inside_GPGP <- as.factor(grouped_by_date_patch$inside_GPGP)
```

```{r fig.width=10}
ggplot(data = grouped_by_date_patch, aes(x = date, y = mmmc, color = inside_GPGP)) +
  geom_point(size = 5, alpha = 0.8) +
  scale_x_date(date_breaks = "5 years",  date_labels = "%Y") +
  scale_y_continuous(limits = c(0, 2100)) +
  scale_color_manual(name = "", values = c("0"="lightblue", "1"="darkred"), labels = c("Outside GPGP", "Inside GPGP")) +
  theme_minimal() +
  #ggtitle("Average month microplastic mass concentration", "among all available observation points in Pacific ocean") +
  ylab("g/km2") +
  xlab("") +
  geom_text(data = data.frame(x = as.Date('1975-01-01'), y = 1750, label = "Average month microplastic\n    mass concentration"), 
            aes(x = x, y = y, label = label), size = 8, family = "Ubuntu Mono", color = "black", hjust = 0) +
  geom_text(data = data.frame(x = as.Date('1972-12-01'), y = 1375, label = "(among all available observation points in Pacific ocean)"), 
            aes(x = x, y = y, label = label), size = 4, family = "Ubuntu Mono", color = "grey30", hjust = 0) +
  theme(
    axis.title = element_text(family = "Ubuntu Mono", size = 15),
    axis.title.x = element_text(vjust = -2),
    axis.title.y = element_text(vjust = 2),
    plot.title = element_text(family = "Ubuntu Mono", size = 20),
    plot.subtitle = element_text(family = "Ubuntu Mono", color = "grey30"),
    legend.title = element_text(family = "Ubuntu Mono"),
    legend.text = element_text(family = 'Ubuntu Mono', size = 15),
    panel.grid.minor.x = element_blank(),
    legend.position = c(0.3,0.53),
    legend.background = element_rect(colour = 'gray90'),
    legend.margin = margin(t = -0.2, r = 0.5, b = 0.3, l = 0.5, unit='cm')
    #plot.margin = unit(c(2,2,2,2), "cm")
  ) + 
  coord_cartesian(clip = 'off')


ggsave('final/mmmc.png', last_plot(), width = 10)
```

  
```{r fig.width=10}
#ggplot(data = grouped_by_date_patch, aes(x = date, y = mmnc, color = inside_GPGP)) + 
# geom_point(size = 2.5) + 
#  ggtitle("Microplastic numerical concentration")

ggplot(data = grouped_by_date_patch, aes(x = date, y = mmnc, color = inside_GPGP)) +
  geom_point(size = 5, alpha = 0.8) +
  scale_x_date(date_breaks = "5 years",  date_labels = "%Y") +
  scale_y_continuous(limits = c(0, 1.1e6), breaks = c(0, 2e5, 4e5, 6e5, 8e5, 1e6), labels = c("0", "2e5", "4e5", "6e5", "8e5", "10e5")) +
  scale_color_manual(name = "", values = c("0"="lightblue", "1"="darkred"), labels = c("Outside GPGP", "Inside GPGP")) +
  #ggtitle("Microplastic mass concentration") +
  theme_minimal() +
  #ggtitle("Average month microplastic numerical concentration", "among all available observation points in Pacific ocean") +
  ylab("n/km2") +
  xlab("") +
  geom_text(data = data.frame(x = as.Date('1975-01-01'), y = 900000, label = "Average month microplastic\n numerical concentration"), 
            aes(x = x, y = y, label = label), size = 8, family = "Ubuntu Mono", color = "black", hjust = 0) +
  geom_text(data = data.frame(x = as.Date('1973-01-01'), y = 750000, label = "(among all available observation points in Pacific ocean)"), 
            aes(x = x, y = y, label = label), size = 4, family = "Ubuntu Mono", color = "grey30", hjust = 0) +
  theme(
    axis.title = element_text(family = "Ubuntu Mono", size = 15),
    axis.text = element_text(family = 'Ubuntu Mono'),
    axis.title.x = element_text(vjust = -2),
    axis.title.y = element_text(vjust = 2),
    plot.title = element_text(family = "Ubuntu Mono", size = 20),
    plot.subtitle = element_text(family = "Ubuntu Mono", color = "grey20"),
    legend.title = element_text(family = "Ubuntu Mono"),
    legend.position = c(0.3,0.55),
    legend.text = element_text(family = 'Ubuntu Mono', size = 15),
    panel.grid.minor.x = element_blank(),
    legend.background = element_rect(colour = 'gray90'),
    legend.margin = margin(t = -0.2, r = 0.5, b = 0.3, l = 0.5, unit='cm')
    #plot.margin = unit(c(2,2,2,2), "cm")
  ) + 
  coord_cartesian(clip = 'off')
 

ggsave('final/mmnc.png', last_plot(), width = 10)
```

```{r}
worldmap <- map("world", fill=TRUE, plot=FALSE)
WGS84 <- CRS("+proj=longlat +datum=WGS84")
worldmapPolys <- map2SpatialPolygons(worldmap, 
                                     IDs=sapply(strsplit(worldmap$names, ":"), "[", 1L), 
                                     proj4string=WGS84)
shift <- 180+20

XY <- data.table(map_data(as(worldmapPolys, "SpatialPolygonsDataFrame")))
XY[, long.new := long + shift]
XY[, long.new := ifelse(long.new > 180, long.new-360, long.new)]
XY[, to.split := sum(diff(long.new) > 300, na.rm=TRUE) > 0, by=group]
XY[, gr.split := ifelse(to.split & long.new < 0, paste0(group, ".", 1), group)]
```

```{r}
garbadge_patch$lon <- 180 + garbadge_patch$lon
```

```{r}
garbadge_patch$inside_GPGP <- as.factor(garbadge_patch$inside_GPGP)
```


```{r}
garbadge_patch[order(garbadge_patch$lon, decreasing = TRUE), ]
```



```{r}
temp <- garbadge_patch[garbadge_patch$lon < 300, ]
```


```{r}
breaks = c(-Inf, 100, 1000, 10000, 100000, 1000000, Inf)

temp$mnc_factors <- cut(temp$microplastic_numerical_concentration_n_km2, breaks = breaks)
```

```{r}
levels(temp$mnc_factors)[1] = "<100"
levels(temp$mnc_factors)[2] = "100-1000"
levels(temp$mnc_factors)[3] = "1000-10000"
levels(temp$mnc_factors)[4] = "10000-100000"
levels(temp$mnc_factors)[5] = "100000-1000000"
levels(temp$mnc_factors)[6] = ">1000000"
```


```{r}
temp
```

```{r fig.width=14}
ggplot() + 
    geom_polygon(data=XY, 
                 aes(x=long.new, y=lat, group=gr.split), 
                 colour="grey70", fill="grey70", size = 0.25) + 
  scale_x_continuous(limits = c(-30, 110), breaks = c(-20, 0, 20, 40, 60, 80, 100), labels = c("160°", "0°", "-160°", "-140°", "-120°", "-100°", "-80°")) +
  scale_y_continuous(limits = c(-30, 85), breaks = c(-20, 0, 20, 40, 60, 80), labels = c("-20°", "0°", "20°", "40°", "60°", "80°")) +
  geom_point(data = temp, 
             aes(x = lon, 
                 y = lat, 
                 colour = mnc_factors), 
             alpha = 0.8,
             size = 2) +
  scale_color_manual("n/km2", values = c("<100"="lightblue", "100-1000"="green", "1000-10000"="yellow", "10000-100000"="orange", "100000-1000000"="red", ">1000000"="darkred")) +
  theme_light() +
  coord_fixed() +
  theme() +
  geom_text(data = data.frame(x = -12, y = 85, label = "Microplastic numerical concentration in Pacific ocean"), 
            aes(x = x, y = y, label = label), size = 6, family = "Ubuntu Mono", color = "black", hjust = 0) +
  geom_text(data = data.frame(x = -10, y = 15, label = "Great Paciific\nGarbage Patch"), aes(x = x, y = y, label = label), size = 4, family = "Ubuntu Mono", color = "black") +
   geom_curve(aes(x = -1, y = 19, xend = 15, yend = 30),
             arrow = arrow(length = unit(0.2, 'cm')), curvature = -0.1, size = 0.4, colour = 'black') + 
  theme(
    axis.text = element_text(family = 'Ubuntu Mono'),
    #axis.title = element_text(family = 'Ubuntu Mono'),
    plot.title = element_text(family = "Ubuntu Mono", size = 15),
    legend.title = element_text(family = "Ubuntu Mono", size = 12),
    legend.text = element_text(family = "Ubuntu Mono"),
    axis.title = element_blank(),
    legend.position = c(0.85, 0.68),
    legend.background = element_blank(),
    legend.key = element_blank(),
    plot.margin = unit(c(1,1,1,1), "cm")
  )

#ggsave("final/pacific_garbage_patch.png", last_plot(), width = 10, height = 10)
```

```{r fig.width=14}
ggplot() + 
    geom_polygon(data=XY, 
                 aes(x=long.new, y=lat, group=gr.split), 
                 colour="#9e9c8e", fill="#9e9c8e", size = 0.25) + 
  scale_x_continuous(limits = c(-30, 110), breaks = c(-20, 0, 20, 40, 60, 80), labels = c("160°", "0°", "-160°", "-140°", "-120°", "-100°")) +
  scale_y_continuous(limits = c(-30, 85), breaks = c(-20, 0, 20, 40, 60, 80), labels = c("-20°", "0°", "20°", "40°", "60°", "80°")) +
  geom_point(data = temp[temp$mictoplastic_mass_concentration_g_km2 > 0, ], 
             aes(x = lon, 
                 y = lat, 
                 color = mictoplastic_mass_concentration_g_km2, 
                 size = mictoplastic_mass_concentration_g_km2), 
             alpha = 0.3) +
  scale_color_gradient2(low = "white", mid = "#02f4cc", high = "#8c2d04", midpoint = 5000) +
  theme_minimal() +
  coord_fixed() +
  theme() +
  xlab("Longtitude") +
  ylab("Latitude") +
  theme(
    axis.text = element_text(family = 'Ubuntu Mono'),
    axis.title = element_text(family = 'Ubuntu Mono')
  )
```


```{r}
marine_pollution <- read_excel('papers/data/marine_pollution.xlsx')
```

```{r}
names(marine_pollution) <- c("country", "economic_status", "coastal_population", "waste_generation_rate_kg_person_day", "%_of_plastic_in_waste_stream", 
                             "%_of_indequately_managed_waste", "%_of_littered_waste", "waste_generation_kg_day", "plastic_waste_generation_kg_day",
                             "inadequately_managed_plastic_waste_kg_day", "plastic_waste_littered_kg_day", "mismanaged_plastic_waste_kg_person_day",
                             "mismanaged_plastic_waste_in_2010_tonnes", "mismanaged_plastic_waste_in_2025_tonnes")
```



```{r}
world <- map_data('world')
```

```{r}
hf <- left_join(XY, marine_pollution, by = c('region' = 'country'))
```

```{r}
hf$mismanaged_plastic_waste_kg_person_day <- as.double(hf$mismanaged_plastic_waste_kg_person_day)
```

```{r}
hf %>% group_by(region) %>% summarise(cos_pop = min(coastal_population))
```

```{r}
summary(hf$mismanaged_plastic_waste_in_2010_tonnes)
```



```{r}
breaks = c(-Inf, 5000, 10000, 50000, 100000, 500000, 1000000, Inf)

hf$factors_mismanaged_2010_tonnes <- cut(hf$mismanaged_plastic_waste_in_2010_tonnes, breaks = breaks)
```

```{r}
levels(hf$factors_mismanaged_2010_tonnes)[1] = "<5000"
levels(hf$factors_mismanaged_2010_tonnes)[2] = "5000-10000"
levels(hf$factors_mismanaged_2010_tonnes)[3] = "10000-50000"
levels(hf$factors_mismanaged_2010_tonnes)[4] = "50000-100000"
levels(hf$factors_mismanaged_2010_tonnes)[5] = "100000-500000"
levels(hf$factors_mismanaged_2010_tonnes)[6] = "500000-1000000"
levels(hf$factors_mismanaged_2010_tonnes)[7] = ">1000000"
```


```{r}
hf
```
```{r}
levels(hf$factors_mismanaged_2010_tonnes)
```

```{r}
countries <- hf[!(hf$region %in% c("Antarctica", "Greenland", "Iceland")), ] 

# aggregate names
cnames <- aggregate(cbind(long.new, lat) ~ mismanaged_plastic_waste_kg_person_day, data = countries, FUN = function(x)mean(range(x)))
```

```{r}
cnames$mismanaged_plastic_waste_kg_person_day <- round(cnames$mismanaged_plastic_waste_kg_person_day, 3)
```

```{r}
cnames$mismanaged_plastic_waste_kg_person_day
```
```{r}
SAG <- read.csv("dots/SAG.csv")
SPG <- read.csv("dots/SPG.csv")
NPG0 <- read.csv("dots/NPG0.csv")
NPG00 <- read.csv("dots/NPG0.csv")
NPG2 <- read.csv("dots/NPG2.csv")
NPG1 <- read.csv("dots/NPG1.csv")
NAG <- read.csv("dots/NAG.csv")
IOG <- read.csv("dots/IOG.csv")

llt <- c(as.integer(runif(n = 20, min = -35, max = -27)),as.integer(runif(n = 20, min = -32, max = -20)),as.integer(runif(n = 10, min = -28, max = -25)))
llg <- c(as.integer(runif(n = 20, min = -90, max = -78)),as.integer(runif(n = 20, min = -82, max = -70)),as.integer(runif(n = 10, min = -80, max = -75)))
buyos2 <- data.frame("platform_code"=seq(1,50))
buyos2$longitude <- llg
buyos2$latitude <- llt
buyos2$Set <- as.integer(runif(n = 50, min = 1, max = 4))
```


```{r fig.width=14}
ggplot(data = countries, aes(x = long.new, y = lat)) +
  geom_polygon(aes(group = group, fill = factors_mismanaged_2010_tonnes)) +
  scale_fill_manual("", 
                    values = c("<5000"="#9bf902", 
                               "5000-10000"="#ccf902", 
                               "10000-50000"="#edf902", 
                               "50000-100000"="#f9d802", 
                               "100000-500000"="#f9bc02", 
                               "500000-1000000"="#f99202",
                               ">1000000"="#f96d02"),
                    labels = c("<5000", "5000-10000", "10000-50000", "50000-100000", "100000-500000", "500000-1000000", ">1000000", "No data"),
                    na.value = "grey90") +
  theme_minimal() +
  coord_fixed() + 
  scale_x_continuous(breaks = c(-160, -140, -120, -100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100, 120, 140, 160),
                     labels = c("20°", "","60°", "", "100°", "", "140°", "", "0°", "", "-140°", "", "-100°", "", "-60°", "", "-20°")) +
  scale_y_continuous(breaks = c(-60, -40, -20, 0, 20, 40, 60, 80),
                     limits = c(-65,95),
                     labels = c("", "-40°", "", "0°", "", "40°", "", "80°")) +
  theme(
    #panel.grid.major = element_line(size = 0.2),
    #panel.grid.minor = element_line(size = 0.1),
    legend.text = element_text(family = "Ubuntu Mono"),
    #axis.text = element_text(family = "Ubuntu Mono"),
    legend.position = c(0.06, 0.45),
    legend.title = element_text(family = "Ubuntu Mono"),
    plot.title = element_text(family = "Ubuntu"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank()
  ) +
  ylab("") +
  xlab("") +
  #ggtitle("Mismatched plastic waste per coutry per year, tonnes") +
  geom_text(data = data.frame(x = -20, y = 95, label = "Mismatched plastic waste per coutry per year, tonnes"), 
            aes(x = x, y = y, label = label), size = 5, family = "Ubuntu Mono", color = "black") +
  geom_point(data = SAG, 
              aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19) +
  stat_ellipse(data = SAG,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  geom_point(data = SPG, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = SPG,
              aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f96d02", level = 0.99) + 
  geom_point(data = NPG0, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = NPG0,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f96d02", level = 0.99) + 
  geom_point(data = NPG00, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = NPG00,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  geom_point(data = NPG1, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = NPG1,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  geom_point(data = NPG2, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = NPG2,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  geom_point(data = NAG, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = NAG,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  geom_point(data = IOG, 
             aes(x = longitude, y = latitude, color=as.factor(Set)), na.rm=TRUE, size=0.5, shape=19)+
  stat_ellipse(data = IOG,
               aes(x = longitude, y = latitude), type = 'norm', alpha = 0.3, color = "#f90212", level = 0.99) + 
  scale_colour_manual(values = c("grey50", "grey60", "grey70"))+ guides(color = FALSE) +
  geom_rect(aes(xmin=-10, xmax=90, ymin=-10, ymax=65), color="black", alpha=0, size = 0.05, linetype = 3) +
  geom_text(x=35, y=12, size=3, label="GREAT PACIFIC \nGARBAGE PATCH", family = "Ubuntu Mono", color="grey40") +
  geom_text(x=-174, y=-50, size=3, label="SOUTH ATLANTIC \nGARBAGE PATCH", family = "Ubuntu Mono", color="grey40") +
  geom_text(x=-80, y=-49, size=3, label="INDIAN OCEAN \nGARBAGE PATCH", family = "Ubuntu Mono", color="grey40") +
  geom_text(x=70, y=-49, size=3, label="SOUTH PACIFIC\nGARBAGE PATCH", family = "Ubuntu Mono", color="grey40") +
  geom_text(x=152, y=13, size=3, label="NORTH \nATLANTIC\nGARBAGE\nPATCH", family = "Ubuntu Mono", color="grey40")


# + 
 # geom_text_repel(data=cnames[cnames$mismanaged_plastic_waste_kg_person_day >= 0.1, ], 
  #                aes(long.new, lat, label = mismanaged_plastic_waste_kg_person_day), 
   #               size=3, 
    #              family = "Ubuntu Mono")


ggsave("final/mismatched.png", last_plot( ), width = 14)
```

```{r}
marine_pollution = marine_pollution[1:192,]
```

```{r}
marine_pollution$economic_status[marine_pollution$economic_status == "HIC"] <- "High-income"
marine_pollution$economic_status[marine_pollution$economic_status == "UMI"] <- "Upper mid-income"
marine_pollution$economic_status[marine_pollution$economic_status == "LMI"] <- "Lower mid-income"
marine_pollution$economic_status[marine_pollution$economic_status == "LI"] <- "Low income"
```



```{r}
grouped_by_status <- marine_pollution %>% group_by(economic_status) %>% summarize(avg__pl_wast_gen_kg_day = mean(plastic_waste_generation_kg_day), avg_inad_man_pl_wast_kg_day = mean(inadequately_managed_plastic_waste_kg_day))
```

```{r}
grouped_by_status$economic_status <- factor(grouped_by_status$economic_status, levels = grouped_by_status$economic_status[order(grouped_by_status$avg__pl_wast_gen_kg_day)])
```

```{r}
names(grouped_by_status) = c('economic_status', " Average plastic waste generation per day, kg", " Average amount of inadequately manged plastic waste per day, kg")
```

```{r fig.width=10}
grouped_by_status %>% 
  gather("metric", "value", -economic_status) %>%
  ggplot(aes(economic_status, value, fill = metric)) +
  scale_y_continuous(breaks = c(0, 200000, 400000, 600000, 800000, 1000000, 1200000, 1400000, 1600000, 1800000),
                     labels = c("0" ,"2e5", "4e5", "6e5", "8e5", "1e6", "1.2e6", "1.4e6", "1.6e6", "1.8e6")) + 
  geom_col(position = "dodge", width = 10) +
  scale_fill_manual(values = c(rep(c("#f9bc02", "#9bf902")))) +
  facet_wrap(~economic_status, scales = "free_x", nrow = 1) +
  theme_minimal() +
  theme(
    strip.placement = "bottom",
    legend.position = c(0.5, 1.2),
    legend.title = element_blank(),
    legend.text = element_text("Ubuntu mono"), 
    axis.title = element_text("Ubuntu mono"),
    strip.text.x = element_blank(),
    panel.spacing.x=unit(4, "lines"),
    panel.spacing.y=unit(1, "lines"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_text(family = "Ubuntu Mono", color = "black", size = 11),
    axis.title.x.bottom = element_text(family = 'Ubuntu Mono', vjust = -5, size = 13),
    axis.title.y.left = element_text(family = 'Ubuntu Mono', vjust = 5, size = 13),
    plot.margin = unit(c(3,3,3,3), "cm")
  ) + 
  ylab("kg/day") +
  xlab("Economic status")

ggsave('final/histogram.png', last_plot(), width = 10)
```