exercises_rlab_01.Rmd

---
title: "Exercises Laboratory Session 01"
author: "Nicola Zomer"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document: 
    default
---

I need to use this setup to obtain nice graphics when I convert R markdown to HTML. 

```{r setup-chunk}
knitr::opts_chunk$set(dev = "ragg_png")

```


```{r, message=FALSE}
# tidyverse
library(tidyverse)
# library(ggplot2)
# library(readr)
# library(tibble)
library(glue)
# library(dplyr)

# others
library(gridExtra)
library(latex2exp)
library(scales)
library(kableExtra)

```

In the whole notebook I will use *kable* to show the tables, as it allows to keep the table formatting also in the html file, which I use to share my work on GitHub. For a similar reason, I will never print an entire table. In the future exercises I will define a function to do it. 

# Exercise 1 - Vectors and data frames

```{r}

filepath = "../../data/lochs_of_Scotland.csv"

df_lakes <- read.csv(filepath, header=TRUE, sep=",")
df_lakes %>%
  kbl() %>%
  kable_styling()

```

```{r}

# Remove columns not containing volume or area and in [mi]
df_lakes <- df_lakes[, c(1, 2, 4)]

# Rename the columns
colnames(df_lakes) <- c('Loch', 'Volume', 'Area')
df_lakes %>%
  kbl() %>%
  kable_styling()

```

## 1. Evaluate the highest and lowest volume and area lake

```{r}

# indices 
idx_max_vol <- which.max(df_lakes$Volume)
idx_min_vol <- which.min(df_lakes$Volume)
idx_max_area <- which.max(df_lakes$Area)
idx_min_area <- which.min(df_lakes$Area)

# results 

cat(
  'Highest volume lake:', df_lakes$Loch[idx_max_vol], ', with volume ', df_lakes$Volume[idx_max_vol], 'km^3\n',
  'Lowest volume lake:', df_lakes$Loch[idx_min_vol], ', with volume ', df_lakes$Volume[idx_min_vol], 'km^3\n', 
  'Highest area lake:', df_lakes$Loch[idx_max_area], ', with area ', df_lakes$Area[idx_max_area], 'km^2\n',
  'Lowest area lake:', df_lakes$Loch[idx_min_area], ', with area ', df_lakes$Area[idx_min_area], 'km^2'
  
  )

```

## 2. Order the frame with respect to the area and determine the two largest area lakes

Ordered data frame:

```{r}

lakes_byarea <- df_lakes[order(df_lakes$Area, decreasing=TRUE), ]
lakes_byarea %>%
  kbl() %>%
  kable_styling()

```

```{r}

largest_area_2 = lakes_byarea$Loch[1:2]
glue('2 largest area lakes: {largest_area_2[1]} and {largest_area_2[2]}')

```

## 3. By summing up the areas occupied by the lakes, determine the area of Scotland covered by water

```{r}

area_water <- sum(df_lakes$Area)
glue('Area of Scotland covered by water: {area_water} km^2')

```

Reference and data: <https://en.wikipedia.org/wiki/List_of_lochs_of_Scotland>

# Exercise 2 - Crude Oil Production

## 1. Write R code that is able to read the file and import it in a data frame structure

The last column of the data frame contains data on crude oil prices from 1861 to 2020, measured in US dollars per barrel.

```{r}

filepath="../../data/crude-oil-prices.csv"

df_oilprices <- read.csv(filepath, header=TRUE, sep=",")
colnames(df_oilprices) <- c(names(df_oilprices[1:3]), 'Price') 

str(df_oilprices)

```

## 2. Produce a plot with the Oil price as a function of the year

```{r}

gg <- ggplot(df_oilprices, aes(x=Year, y=Price))+
  geom_line(col='navyblue', size=0.8) +
  labs(title="Crude Oil Prices from 1861 to 2020 ($/barrel)", 
       y="Price",
       x='Year', 
       caption = "Source: https://ourworldindata.org/grapher/crude-oil-prices") +
  scale_y_continuous(
    breaks = seq(0, 120, 15),
    minor_breaks = NULL
  ) + 
  scale_x_continuous(
    breaks = seq(1850, 2050, 25),
    minor_breaks = NULL,
    limits=c(1850, 2025)
  )+
  theme_bw()

plot(gg)

```

## 3. Which is the highest price in history ? When did it occur ?

```{r}

highest_price <- max(df_oilprices$Price)
highest_price_year <- df_oilprices$Year[which.max(df_oilprices$Price)]

glue('Highest price in hystory: {format(highest_price, digits=5)} $/barrel.\n
  It occured in {highest_price_year}.')

```

## 4. Plot the derivative of the curve, simply evaluated with the finite difference formula

$$
\frac{\partial price}{\partial year} = price_{j+1}-price{j}
$$

```{r}

prices <- df_oilprices$Price
years <- df_oilprices$Year
derivatives <- prices[2:length(prices)]-prices[1:length(prices)-1]

df_derivatives <- data.frame(years[1:length(years)-1], derivatives)
colnames(df_derivatives) <- c('Year', 'Derivative')

str(df_derivatives)

```

```{r}

gg <- ggplot(df_derivatives, aes(x=Year, y=Derivative))+
  geom_line(col='navyblue', size=0.8) + 
  labs(title="Annual Variation of Crude Oil Prices from 1861 to 2020 ($/barrel)", 
       y=TeX("$\\Delta$Price"),
       x='Year', 
       caption = "Source: https://ourworldindata.org/grapher/crude-oil-prices") +
  scale_x_continuous(
    breaks = seq(1850, 2050, 25),
    minor_breaks = NULL,
    limits=c(1850, 2025)
  ) +
  scale_y_continuous(
    breaks = seq(-45, 45, 15),
    minor_breaks = NULL,
  )+
  theme_bw()

plot(gg)

```

Reference and data: <https://ourworldindata.org/grapher/crude-oil-prices>

# Exercise 3 - World Coal Production

## 1. Write R code that is able to read the file and import it in a tibble structure

```{r}

filepath = "../../data/coal-production-by-country.csv"

coal_prod <- read_csv(filepath, col_names=TRUE)
colnames(coal_prod) <- c(names(coal_prod[1:3]), 'Production')

glimpse(coal_prod)

```

```{r}
is_tibble(coal_prod)

```

## 2. Count the number of countries available in the file and produce a barplot with the number of entries for each country

```{r}
countries <- unique(coal_prod$Entity)
cat('Number of countries:', length(countries))

```
Notice that if we use `Code` instead of `Entity` we get a different result.

```{r}
cat("Number of countries using 'Code' instead of 'Entity':", length(unique(coal_prod$Code)))

```
This is due to the fact that some entities, such as continents, do not have a corresponding code and this results in a `NA` entry in the tibble. As I want to consider also the continents in the following analysis, I prefer to define the countries vector using the `Entity` attribute.

```{r, fig.height = 25, fig.width = 7}

gg <- ggplot(data=coal_prod, aes(x=Entity)) +
  geom_bar(stat = "count", width=0.7, fill="dodgerblue3", alpha=0.8) + 
  coord_flip() +
  scale_x_discrete(limits=rev) +
  labs(title="Number of entries for each country", 
       x="Country",
       y="Count", 
       caption = "Source: https://ourworldindata.org/grapher/coal-production-by-country") +
  theme_bw()

gg


```

For the following items select only the years $$\geq$$ 1970.

## 3. Selecting only the year after 1970, determine the total integrated production for each country and print the top 5 countries with highest coal productions

```{r}

# select only the years after 1970
coal_prod_1970 <- coal_prod[coal_prod['Year']>1970, ]

# total integrated production for each country
int_prod_1970 <- aggregate(coal_prod_1970$Production, by=list(coal_prod_1970$Entity), FUN=sum)
colnames(int_prod_1970) <- c('Country', 'IntProd')

# top 5 countries
int_prod_1970[order(int_prod_1970$IntProd, decreasing = TRUE)[1:5], ] %>%
  kbl() %>%
  kable_styling()

```

As some elements are not countries, I want to remove these rows, and repeat the process until I get 5 countries.

```{r}

removed_elements = c('World', 'Asia Pacific', 'Asia and Oceania', 'OECD',
                     'North America', 'Eurasia', 'Europe', 'EU-28', 'CIS', 'Africa',
                     'South Africa'
                     ) # not all non-countries entries

int_prod_1970 <- int_prod_1970[!int_prod_1970$Country %in% removed_elements, ]

# top 5 countries
int_prod_1970[order(int_prod_1970$IntProd, decreasing = TRUE)[1:5], ] %>%
  kbl() %>%
  kable_styling()

top5_countries <- int_prod_1970$Country[order(int_prod_1970$IntProd, decreasing = TRUE)[1:5]]


```

## 4. For the 5 top Countries, create a plot of production as a function of time

```{r}

# create the reduced dataframe
df_top5 <- coal_prod_1970[coal_prod_1970$Entity %in% top5_countries, ]

# inspect the structure of the dataframe
str(df_top5)

```

```{r}

colors <- c('darkred', 'deepskyblue', 'darkgreen', 'darkgoldenrod', 'darkblue')

ggplot(df_top5, aes(x = Year, y = Production)) + 
  geom_line(aes(color = Entity), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title="Production vs time for the top 5 countries", 
       y="Coal Production (TWh)",
       caption = "Source: https://ourworldindata.org/grapher/coal-production-by-country") +
  theme_bw()

```

## 5. Generate a plot with the cumulative sum of the World's coal production over the years

```{r}

world_cum <- cumsum(coal_prod_1970$Production[coal_prod_1970$Entity=='World'])
years <- coal_prod_1970$Year[coal_prod_1970$Entity=='World']

ggplot() + 
  geom_line(aes(x=years, y=world_cum), col='navyblue', size=1) + 
  labs(title="Cumulative World's coal production", 
       x="Year",
       y="Production", 
       caption = "Source: https://ourworldindata.org/grapher/coal-production-by-country") +
  theme_bw()

```

Reference and data: <https://ourworldindata.org/grapher/coal-production-by-country>

# Exercise 4 - Covid19 Vaccine data

## File 'vaccinations-by-manufacturer.csv'

### 1. Filter() the original tibble by selecting the following countries: Italy

```{r}

# load the data
urlfile = "https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/vaccinations/vaccinations-by-manufacturer.csv"
vacc_bymanu <- read_csv(url(urlfile), col_names=TRUE)

glimpse(vacc_bymanu)

```

```{r}

# select only data related to Italy
vacc_Italy <- vacc_bymanu[vacc_bymanu$location=='Italy',]
head(vacc_Italy) %>%
  kbl() %>%
  kable_styling()

```

### 2. Plot the number of vaccines given as a function of time for the different vaccine manufacturer

```{r}

# inspect the manufacturers
unique(vacc_Italy$vaccine)

```

```{r}

colors <- c('darkred', 'deepskyblue', 'darkgreen', 'darkgoldenrod', 'darkblue')

gg <- ggplot(vacc_Italy, aes(x = date, y = total_vaccinations)) + 
  geom_line(aes(color = vaccine), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title="Number of vaccines for different manufacturers in Italy", 
       x='Date',
       y='Number of vaccinations',
       color='Vaccine', 
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_y_continuous(trans='log10') +
  scale_x_date(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
    date_labels="%b-%y")+
  theme_bw()
  
plot(gg)

```

### 3. From the same tibble plot the total number of vaccines shot per day in Italy

```{r}

tot_vacc_Italy <- aggregate(vacc_Italy$total_vaccinations, by=list(vacc_Italy$date), FUN=sum)
colnames(tot_vacc_Italy) <- c('date', 'tot_vaccinations')


gg <- ggplot(tot_vacc_Italy, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Total number of vaccines in Italy", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()
  
plot(gg)

```

Probably due to the the missing information about Orford/AstraZeneca, after a certain date, the cumulative number of vaccines decreases at a certain point and does strange oscillations. This is absurd and cannot be justified.

This implies that if we try to get the number of vaccines per day by taking the difference between consecutive days, we get negative results, that is also absurd. Let's fix this.

First of all, I want to find the exact dates on which the number of total vaccinations decreases.

```{r}
# tail(vacc_Italy$date[vacc_Italy$vaccine=='Oxford/AstraZeneca'], 1)

j=0
problematic_dates = c()
for (i in seq_along(tot_vacc_Italy$tot_vaccinations)){
  if (tot_vacc_Italy$tot_vaccinations[i]<j){
    problematic_dates <- append(problematic_dates, tot_vacc_Italy$date[i])
  }
  j=tot_vacc_Italy$tot_vaccinations[i]
}
kable(problematic_dates, col.names="Problematic Dates") %>%
  kable_styling()

```

Let's see if we know the total Oxford/Astrazeneca vaccinations on these days.

```{r}

for (i in seq_along(problematic_dates)){
  tmp <- vacc_Italy[(vacc_Italy$date==problematic_dates[i]) & (vacc_Italy$vaccine=='Oxford/AstraZeneca'), ]
  is_not_empty <- as.logical(nrow(tmp))
  print(paste(problematic_dates[i], ': is it NOT empty?', is_not_empty))
}

```

```{r}
print(paste('Last information about Oxford/AstraZeneca vaccine in ', 
            tail(vacc_Italy$date[vacc_Italy$vaccine=='Oxford/AstraZeneca'], 1)))
```

As expected I have no information on Oxford/AstraZeneca vaccines on these dates and after Jan 22. This causes the graph of the total number of vaccinations to drop. Moreover, notice how the same problem arises even after this date, a sign that some data relating to other pharmaceutical companies is also missing. However, I don't want to take these into consideration too. A possible solution would be to define a new dataframe of zeros, containing all the dates, and fill it with the information provided, expanding them in the missing days. 

**Instead, this is what I will try to fix the problem, considering only the missing data of AstraZeneca:**

1.  **Break the plot in January 2022.** After this date in my opinion it has no meaning to show the total number of vaccinations. One possibility would be to extend the AstraZeneca plot, but this would mean assuming zero vaccinations of that type after January 2022. Since this is a strong assumption, I prefer to stop the study first, using only data that I believe to be reliable.
2.  **Add the total number of vaccinations in the dates before January 2022 in which I have no information**. I set them to be equal to the value on the day before. Hopefully, being single days, this does not alter the results too much.

```{r}

# cut the dataframe
vacc_IT_fixed = vacc_Italy[vacc_Italy$date < '2022-01-01', ]

# add missing data
for (date in seq(as.Date("2021-01-14"), as.Date("2021-12-31"), by='day')){
  if (nrow(vacc_IT_fixed[(vacc_IT_fixed$date==date) & (vacc_IT_fixed$vaccine=='Oxford/AstraZeneca'),])==0){
    
    tmp <- vacc_IT_fixed[(vacc_IT_fixed$date==date-1) & (vacc_IT_fixed$vaccine=='Oxford/AstraZeneca'),]
    tmp$date = tmp$date+1
    vacc_IT_fixed[nrow(vacc_IT_fixed)+1, ] <- tmp

  }
}

# check
for (i in seq_along(problematic_dates)[1:11]){
 tmp <- vacc_IT_fixed[
   (vacc_IT_fixed$date==problematic_dates[i]) & 
     (vacc_IT_fixed$vaccine=='Oxford/AstraZeneca'), ]
 is_not_empty <- as.logical(nrow(tmp))
 print(paste(problematic_dates[i], ': is it NOT empty?', is_not_empty))
}

```

Let's see how the cumulative distribution changed.

```{r}

tot_vacc_Italy <- aggregate(vacc_IT_fixed$total_vaccinations, by=list(vacc_IT_fixed$date), FUN=sum)
colnames(tot_vacc_Italy) <- c('date', 'tot_vaccinations')


gg <- ggplot(tot_vacc_Italy, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Total number of vaccines in Italy", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-01-01"), by = "quarter"),
  ) +
  theme_bw()
  
plot(gg)
    
```

It seems OK. Now I can evaluate the number of vaccinations per day.

```{r}

tot_vacc_Italy_perday <- diff(tot_vacc_Italy$tot_vaccinations)
dates <- tot_vacc_Italy$date[2:length(tot_vacc_Italy$date)]

tot_vacc_Italy_perday <- data.frame(dates, tot_vacc_Italy_perday)
colnames(tot_vacc_Italy_perday) <- c('date', 'tot_vaccinations')

gg <- ggplot(tot_vacc_Italy_perday, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Number of vaccines per day in Italy", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-01-01"), by = "quarter"),
  ) +
  theme_bw()
  
plot(gg)

```

### 4. Do the same exercise for the following countries: Germany and United States of America

#### Germany

```{r, message=FALSE}

# select only data related to Germany
vacc_Germany <- vacc_bymanu[vacc_bymanu$location=='Germany',]


# number of vaccines for the different manufacturers
colors <- c('darkred', 'deepskyblue', 'darkgreen', 'darkgoldenrod', 'darkblue')

G_num <- ggplot(vacc_Germany, aes(x = date, y = total_vaccinations)) + 
  geom_line(aes(color = vaccine), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title="Number of vaccines for different manufacturers in Germany", 
       x='Date',
       y='Number of vaccinations',
       color='Vaccine', 
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_y_continuous(trans='log10') +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()



# total number of vaccines per day in Germany
tot_vacc_Germany <- aggregate(vacc_Germany$total_vaccinations, by=list(vacc_Germany$date), FUN=sum)
colnames(tot_vacc_Germany) <- c('date', 'tot_vaccinations')

tot_vacc_Germany_perday <- diff(tot_vacc_Germany$tot_vaccinations)
dates <- tot_vacc_Germany$date[2:length(tot_vacc_Germany$date)]

tot_vacc_Germany_perday <- data.frame(dates, tot_vacc_Germany_perday)
colnames(tot_vacc_Germany_perday) <- c('date', 'tot_vaccinations')

G_perday <- ggplot(tot_vacc_Germany_perday, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Number of vaccines per day in Germany", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()



# plot the results
plot(G_num)
plot(G_perday)

```

#### United States of America

```{r}

# select only data related to Germany
vacc_USA <- vacc_bymanu[vacc_bymanu$location=='United States',]


# number of vaccines for the different manufacturers
colors <- c('darkred', 'deepskyblue', 'darkgreen', 'darkgoldenrod', 'darkblue')

USA_num <- ggplot(vacc_USA, aes(x = date, y = total_vaccinations)) + 
  geom_line(aes(color = vaccine), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title="Number of vaccines for different manufacturers in USA", 
       x='Date',
       y='Number of vaccinations',
       color='Vaccine', 
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_y_continuous(trans='log10') +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()



# total number of vaccines per day in USA
tot_vacc_USA <- aggregate(vacc_USA$total_vaccinations, by=list(vacc_USA$date), FUN=sum)
colnames(tot_vacc_USA) <- c('date', 'tot_vaccinations')

tot_vacc_USA_perday <- diff(tot_vacc_USA$tot_vaccinations)
dates <- tot_vacc_USA$date[2:length(tot_vacc_USA$date)]

tot_vacc_USA_perday <- data.frame(dates, tot_vacc_USA_perday)
colnames(tot_vacc_USA_perday) <- c('date', 'tot_vaccinations')

USA_perday <- ggplot(tot_vacc_USA_perday, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Number of vaccines per day in USA", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()



# plot the results
plot(USA_num)
plot(USA_perday)

```


Notice that around March 2022 there are some outliers for Johnson & Johnson. These cause an incorrect behavior in the vaccine per day curve, which shows negative values over that period, due to a negative difference in the total number of vaccine given. To solve this problem, I will consider the trend of the total number of Johnson & Johnson vaccines linear over that time interval.

```{r}

# find the initial and final wrong samples
df_JJ <- vacc_USA[vacc_USA$vaccine=='Johnson&Johnson', ]

date = as.Date('2022-03-01')
first_date = as.Date('2022-03-01')
last_date = as.Date('2022-03-01')

first_found = FALSE
last_found = FALSE

while (first_found == FALSE && last_found == FALSE && date<tail(df_JJ$date, 1)){
  if (df_JJ$total_vaccinations[df_JJ$date==date+1]<df_JJ$total_vaccinations[df_JJ$date==date]){
    if (first_found==FALSE){
      first_date = date
      vacc_broken = df_JJ$total_vaccinations[df_JJ$date==date]
      first_found=TRUE
    }
  }
  else if (first_found==TRUE){
    if (df_JJ$total_vaccinations[df_JJ$date==date+1] > vacc_broken){
      last_date = date+1
      print(last_date)
      last_found=TRUE
    }
  }
  date = date + 1
}

# show results
if (first_date == '2022-03-01') {
  print('Strange, there is no problem')
} else if (last_date == '2022-03-01') {
  corr_date = first_date+1
  print(paste('There is only one corrupted sample:', corr_date))
} else {
  print(paste('Problems begin after: ', first_date))
  print(paste('Problems end at: ', last_date))
}


```

```{r}

# fix it
tmp <- vacc_USA[(vacc_USA$vaccine=='Johnson&Johnson') & (vacc_USA$date==corr_date), ]
tmp$total_vaccinations <-
  (vacc_USA$total_vaccinations[(vacc_USA$vaccine=='Johnson&Johnson') & (vacc_USA$date==corr_date-1)]+
  vacc_USA$total_vaccinations[(vacc_USA$vaccine=='Johnson&Johnson') & (vacc_USA$date==corr_date+1)])/2

vacc_USA[(vacc_USA$vaccine=='Johnson&Johnson') & (vacc_USA$date==corr_date), ] <- tmp

# PLOT THE CORRECTED DISTRIBUTIONS
# number of vaccines for the different manufacturers
USA_num <- ggplot(vacc_USA, aes(x = date, y = total_vaccinations)) + 
  geom_line(aes(color = vaccine), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title="Number of vaccines for different manufacturers in USA", 
       x='Date',
       y='Number of vaccinations',
       color='Vaccine', 
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_y_continuous(trans='log10') +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()

# total number of vaccines per day in USA
tot_vacc_USA <- aggregate(vacc_USA$total_vaccinations, by=list(vacc_USA$date), FUN=sum)
colnames(tot_vacc_USA) <- c('date', 'tot_vaccinations')

tot_vacc_USA_perday <- diff(tot_vacc_USA$tot_vaccinations)
dates <- tot_vacc_USA$date[2:length(tot_vacc_USA$date)]

tot_vacc_USA_perday <- data.frame(dates, tot_vacc_USA_perday)
colnames(tot_vacc_USA_perday) <- c('date', 'tot_vaccinations')

USA_perday <- ggplot(tot_vacc_USA_perday, aes(x = date, y = tot_vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Number of vaccines per day in USA", 
       x='Date',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()


# plot the results
plot(USA_num)
plot(USA_perday)

```

Data: <https://github.com/owid/covid-19-data/blob/master/public/data/vaccinations/vaccinations-by-manufacturer.csv>

## File 'vaccinations.csv'

### 1. Selecting all the European countries in the tibble, plot the number of daily vaccinations per million as a function of date

```{r}

# load the data
urlfile = "https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/vaccinations/vaccinations.csv"
df_vacc <- read_csv(url(urlfile), col_names=TRUE)

# european countries
df_vacc_europe = df_vacc[df_vacc$iso_code=='OWID_EUR', ]

# number of vaccinations per million per date
vacc_permillion = aggregate(df_vacc_europe$daily_vaccinations_per_million, list(df_vacc_europe$date), FUN=sum)
colnames(vacc_permillion) <- c('date', 'vaccinations')

gg1 <- ggplot(vacc_permillion, aes(x = date, y = vaccinations)) + 
  geom_line(col='dodgerblue3', size=0.8) + 
  labs(title="Number of daily vaccinations per millions in Europe", 
       x='Date',
       y='Number of vaccinations/day/million',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()

gg2 <- ggplot(vacc_permillion, aes(x = date, y = vaccinations)) + 
  geom_col(fill='dodgerblue3', width=0.6) + 
  labs(title="Number of daily vaccinations per millions in Europe", 
       x='Date',
       y='Number of vaccinations/day/million',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()

# plot(gg1)
plot(gg2)

```

It can be noted that this plot presents the same general trend as the graphs relating to Italy and Germany. Let's plot over this graph the number of daily vaccinations per million of each European country. 

```{r}

# import iso-codes of Europian countries from csvfile

EU_isocodes <- read_csv("../../data/eur_country_codes.csv", col_names=TRUE )

df_vacc_eucountry <- df_vacc[df_vacc$iso_code %in% EU_isocodes$`Alpha-3`, ]

unique(df_vacc_eucountry$location) %>%
  kable(col.names='European Countries in the vaccinations tibble') %>%
  kable_styling()

```

```{r, message=FALSE}

# number of vaccinations per million per date

vacc_permillion_alleu <- summarise(group_by(df_vacc_eucountry, date, iso_code), vacc_permill= sum(daily_vaccinations_per_million))
# head(vacc_permillion_alleu) %>%
#   kable() %>%
#   kable_styling()


gg <- ggplot() +
  geom_line(data=vacc_permillion_alleu, aes(x = date, y = vacc_permill, color=iso_code), size=0.4) +
  scale_colour_grey(start=0, end=0.9)+
  geom_line(data=vacc_permillion, aes(x = date, y = vaccinations), col='firebrick', size=1.4) +
  labs(title="Number of daily vaccinations per millions in Europe",
       x='Date',
       y='Number of vaccinations/day/million',
       color='ISO-CODE',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()

plot(gg)

```



### 2. Study the data structure and produce few relevant plots of your taste

```{r}
cat('Data structure:\n\n')
glimpse(df_vacc)

```

First of all I want to study for the European countries what is **today** the number of: 
- people vaccinated per hundred
- people fully vaccinated per hundred
- total boosters per hundred


```{r}

last_date = tail(df_vacc_eucountry$date,1)
df_today <- na.omit(df_vacc_eucountry[df_vacc_eucountry$date==last_date, ])

plot_1 <- df_today %>%
  mutate(location = fct_reorder(location, people_vaccinated_per_hundred)) %>%
  ggplot(aes(location, people_vaccinated_per_hundred)) +
  geom_col(fill='dodgerblue3', alpha=0.8) +
  labs(title=paste('People vaccinated per hundred on', last_date),
       x='Country',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data")+
  coord_flip() +
  theme_bw()

plot_2 <- df_today %>%
  mutate(location = fct_reorder(location, people_fully_vaccinated_per_hundred)) %>%
  ggplot(aes(location, people_fully_vaccinated_per_hundred)) +
  geom_col(fill='dodgerblue3', alpha=0.8) +
  labs(title=paste('People fully vaccinated per hundred on', last_date),
       x='Country',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data")+
  coord_flip() +
  theme_bw()

plot_3 <- df_today %>%
  mutate(location = fct_reorder(location, total_boosters_per_hundred)) %>%
  ggplot(aes(location, total_boosters_per_hundred)) +
  geom_col(fill='dodgerblue3', alpha=0.8) +
  labs(title=paste('Total boosters per hundred on', last_date),
       x='Country',
       y='Number of vaccinations',
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data")+
  coord_flip() +
  theme_bw()

plot(plot_1)
plot(plot_2)
plot(plot_3)


```

I want to try to superimpose the 3 plots. I don't think it's useful from a data visualization point of view, but I believe it is a great opportunity to test how to do this in R. 

```{r}

A <- df_today %>%
  select(location, people_vaccinated_per_hundred) %>%
  rename('n'='people_vaccinated_per_hundred')

A['type'] = rep(c('People Vaccinated'), length(A$location))


B <- df_today %>%
  select(location, people_fully_vaccinated_per_hundred)  %>%
  rename('n'='people_fully_vaccinated_per_hundred')

B['type'] = rep(c('People Fully Vaccinated'), length(B$location))

C <- df_today %>%
  select(location, total_boosters_per_hundred)  %>%
  rename('n'='total_boosters_per_hundred')

C['type'] = rep(c('Total Boosters'), length(C$location))

new <- rbind(A, B, C)

ggplot(new) +
  geom_col(aes(x=location, y=n, fill=type)) +
  labs(title=paste("Number of vaccinations per hundred on", last_date),
      x='Country',
      y='Number of vaccinations per hundred',
      fill = NULL,
      caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data")+
  coord_flip()+
  theme(legend.title=element_blank())+
  theme_bw()

```

Now, I want to do again the same thing in a fancier way, using the library `highcharter`. This library could be very useful in the future.  

```{r, message=FALSE}

library(highcharter)
options(highcharter.theme = hc_theme_google())

new %>%
  hchart(
    'bar', hcaes(x='location', y='n', group='type'),
    stacking='normal'
    ) %>%
  hc_title(text=paste('Number of vaccinations per European country per hundred on', last_date)) %>% 
  hc_xAxis(title = list(text = 'Country')) %>%
  hc_yAxis(title = list(text = 'Number of vaccinations per hundred'))

```

Finally, let's see how these 3 metrics change over time in different continents. To do this, I prefer to define a function, to avoid code replication.

```{r, warning = FALSE}

continents <- c("Europe", "Asia", "Africa", "North America", "South America")
df_vacc_ <- df_vacc[df_vacc$location %in% continents, ]

myplot <- function(y, title_metric){df_vacc_ %>%
  ggplot() +
  geom_line(aes(x=date, y=y, color = location), size=0.6) + 
  scale_color_manual(values = colors)+ 
  labs(title=paste(title_metric, 'per hundred vs time'), 
       x='Date',
       y='Number of vaccinations',
       color='Vaccine', 
       caption = "Source: https://github.com/owid/covid-19-data/tree/master/public/data") +
  scale_x_continuous(
    breaks = seq(as.Date("2021-01-01"), as.Date("2022-04-01"), by = "quarter"),
  ) +
  theme_bw()
}

myplot(df_vacc_$people_vaccinated_per_hundred, 'People vaccinated')
myplot(df_vacc_$people_fully_vaccinated_per_hundred, 'People fully vaccinated')
myplot(df_vacc_$total_boosters_per_hundred, 'Total boosters')

```

From these graphs it is evident that there is a difference in the possibility of accessing vaccines in different geographical areas, especially in Africa.

Data: <https://github.com/owid/covid-19-data/blob/master/public/data/vaccinations/vaccinations.csv>