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

knitr::opts_chunk$set(dev = "ragg_png")

# 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

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

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

Loch	Volume.km.3.	Volume.mi.3.	Area.km.2.	Area.mi.2.	Length.km.	Length.mi.	Max..depth.m.	Max..depth.ft.	Mean.depth.m.	Mean.depth.ft.
Loch Ness	7.45	1.790	56.0	22.0	36.2	22.5	227	745	132.0	433
Loch Lomond	2.60	0.620	71.0	27.0	36.0	22.0	190	620	37.0	121
Loch Morar	2.30	0.550	26.7	10.3	18.8	11.7	310	1020	87.0	285
Loch Tay	1.60	0.380	26.4	10.2	23.0	14.0	150	490	60.6	199
Loch Awe	1.20	0.290	39.0	15.0	41.0	25.0	94	308	32.0	105
Loch Maree	1.09	0.260	28.6	11.0	20.0	12.0	114	374	38.0	125
Loch Ericht	1.08	0.260	18.6	7.2	23.0	14.0	156	512	57.6	189
Loch Lochy	1.07	0.260	16.0	6.2	16.0	9.9	162	531	70.0	230
Loch Rannoch	0.97	0.230	19.0	7.3	15.7	9.8	134	440	51.0	167
Loch Shiel	0.79	0.190	19.5	7.5	28.0	17.0	128	420	40.0	130
Loch Katrine	0.77	0.180	12.4	4.8	12.9	8.0	151	495	43.4	142
Loch Arkaig	0.75	0.180	16.0	6.2	19.3	12.0	109	358	46.5	153
Loch Shin	0.35	0.084	22.5	8.7	27.8	17.3	49	161	15.5	51

# 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()

Loch	Volume	Area
Loch Ness	7.45	56.0
Loch Lomond	2.60	71.0
Loch Morar	2.30	26.7
Loch Tay	1.60	26.4
Loch Awe	1.20	39.0
Loch Maree	1.09	28.6
Loch Ericht	1.08	18.6
Loch Lochy	1.07	16.0
Loch Rannoch	0.97	19.0
Loch Shiel	0.79	19.5
Loch Katrine	0.77	12.4
Loch Arkaig	0.75	16.0
Loch Shin	0.35	22.5

1. Evaluate the highest and lowest volume and area lake

# 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'
  
  )

## Highest volume lake: Loch Ness , with volume  7.45 km^3
##  Lowest volume lake: Loch Shin , with volume  0.35 km^3
##  Highest area lake: Loch Lomond , with area  71 km^2
##  Lowest area lake: Loch Katrine , with area  12.4 km^2

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

Ordered data frame:

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

	Loch	Volume	Area
2	Loch Lomond	2.60	71.0
1	Loch Ness	7.45	56.0
5	Loch Awe	1.20	39.0
6	Loch Maree	1.09	28.6
3	Loch Morar	2.30	26.7
4	Loch Tay	1.60	26.4
13	Loch Shin	0.35	22.5
10	Loch Shiel	0.79	19.5
9	Loch Rannoch	0.97	19.0
7	Loch Ericht	1.08	18.6
8	Loch Lochy	1.07	16.0
12	Loch Arkaig	0.75	16.0
11	Loch Katrine	0.77	12.4

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

## 2 largest area lakes: Loch Lomond and Loch Ness

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

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

## Area of Scotland covered by water: 371.7 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.

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)

## 'data.frame':    160 obs. of  4 variables:
##  $ Entity: chr  "World" "World" "World" "World" ...
##  $ Code  : chr  "OWID_WRL" "OWID_WRL" "OWID_WRL" "OWID_WRL" ...
##  $ Year  : int  1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 ...
##  $ Price : num  0.49 1.05 3.15 8.06 6.59 3.74 2.41 3.63 3.64 3.86 ...

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

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 ?

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}.')

## Highest price in hystory: 111.67 $/barrel.
## 
## It occured in 2012.

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

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

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)

## 'data.frame':    159 obs. of  2 variables:
##  $ Year      : int  1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 ...
##  $ Derivative: num  0.56 2.1 4.91 -1.47 -2.85 ...

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

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

coal_prod <- read_csv(filepath, col_names=TRUE)

## Rows: 11528 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (2): Entity, Code
## dbl (2): Year, Coal production (TWh)
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

colnames(coal_prod) <- c(names(coal_prod[1:3]), 'Production')

glimpse(coal_prod)

## Rows: 11,528
## Columns: 4
## $ Entity     <chr> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan",~
## $ Code       <chr> "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AF~
## $ Year       <dbl> 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909,~
## $ Production <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~

is_tibble(coal_prod)

## [1] TRUE

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

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

## Number of countries: 200

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

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

## Number of countries using 'Code' instead of 'Entity': 177

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.

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

# 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()

	Country	IntProd
176	World	1260113.2
9	Asia Pacific	690240.5
8	Asia and Oceania	682343.9
32	China	459564.6
120	OECD	428754.2

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

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()

	Country	IntProd
32	China	459564.65
166	United States	226068.84
135	Russia	98251.65
10	Australia	83077.07
77	India	77043.45

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

# 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)

## tibble [200 x 4] (S3: tbl_df/tbl/data.frame)
##  $ Entity    : chr [1:200] "Australia" "Australia" "Australia" "Australia" ...
##  $ Code      : chr [1:200] "AUS" "AUS" "AUS" "AUS" ...
##  $ Year      : num [1:200] 1981 1982 1983 1984 1985 ...
##  $ Production: num [1:200] 789 805 839 881 1074 ...

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

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

# 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)

## Rows: 36389 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (2): location, vaccine
## dbl  (1): total_vaccinations
## date (1): date
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(vacc_bymanu)

## Rows: 36,389
## Columns: 4
## $ location           <chr> "Argentina", "Argentina", "Argentina", "Argentina",~
## $ date               <date> 2020-12-29, 2020-12-29, 2020-12-29, 2020-12-29, 20~
## $ vaccine            <chr> "Moderna", "Oxford/AstraZeneca", "Sinopharm/Beijing~
## $ total_vaccinations <dbl> 2, 5, 1, 20481, 2, 5, 1, 40583, 2, 5, 1, 43388, 2, ~

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

location	date	vaccine	total_vaccinations
Italy	2020-12-27	Moderna	2
Italy	2020-12-27	Pfizer/BioNTech	7347
Italy	2020-12-28	Moderna	6
Italy	2020-12-28	Pfizer/BioNTech	8848
Italy	2020-12-29	Moderna	11
Italy	2020-12-29	Pfizer/BioNTech	9942

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

# inspect the manufacturers
unique(vacc_Italy$vaccine)

## [1] "Moderna"            "Pfizer/BioNTech"    "Johnson&Johnson"   
## [4] "Oxford/AstraZeneca" "Novavax"

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

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.

# 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()

Problematic Dates
2021-11-01
2021-11-03
2021-11-05
2021-11-21
2021-11-28
2021-12-01
2021-12-06
2021-12-15
2021-12-18
2021-12-21
2021-12-31
2022-01-02
2022-01-12
2022-03-01
2022-03-10
2022-03-12
2022-03-16
2022-03-20
2022-03-24

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

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))
}

## [1] "2021-11-01 : is it NOT empty? FALSE"
## [1] "2021-11-03 : is it NOT empty? FALSE"
## [1] "2021-11-05 : is it NOT empty? FALSE"
## [1] "2021-11-21 : is it NOT empty? FALSE"
## [1] "2021-11-28 : is it NOT empty? FALSE"
## [1] "2021-12-01 : is it NOT empty? FALSE"
## [1] "2021-12-06 : is it NOT empty? FALSE"
## [1] "2021-12-15 : is it NOT empty? FALSE"
## [1] "2021-12-18 : is it NOT empty? FALSE"
## [1] "2021-12-21 : is it NOT empty? FALSE"
## [1] "2021-12-31 : is it NOT empty? FALSE"
## [1] "2022-01-02 : is it NOT empty? FALSE"
## [1] "2022-01-12 : is it NOT empty? FALSE"
## [1] "2022-03-01 : is it NOT empty? FALSE"
## [1] "2022-03-10 : is it NOT empty? FALSE"
## [1] "2022-03-12 : is it NOT empty? FALSE"
## [1] "2022-03-16 : is it NOT empty? FALSE"
## [1] "2022-03-20 : is it NOT empty? FALSE"
## [1] "2022-03-24 : is it NOT empty? FALSE"

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

## [1] "Last information about Oxford/AstraZeneca vaccine in  2022-01-11"

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:

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.
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.

# 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))
}

## [1] "2021-11-01 : is it NOT empty? TRUE"
## [1] "2021-11-03 : is it NOT empty? TRUE"
## [1] "2021-11-05 : is it NOT empty? TRUE"
## [1] "2021-11-21 : is it NOT empty? TRUE"
## [1] "2021-11-28 : is it NOT empty? TRUE"
## [1] "2021-12-01 : is it NOT empty? TRUE"
## [1] "2021-12-06 : is it NOT empty? TRUE"
## [1] "2021-12-15 : is it NOT empty? TRUE"
## [1] "2021-12-18 : is it NOT empty? TRUE"
## [1] "2021-12-21 : is it NOT empty? TRUE"
## [1] "2021-12-31 : is it NOT empty? TRUE"

Let’s see how the cumulative distribution changed.

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.

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

# 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)

## Warning: Transformation introduced infinite values in continuous y-axis

plot(G_perday)

United States of America

# 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.

# 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))
}

## [1] "There is only one corrupted sample: 2022-03-16"

# 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

# 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)

## Rows: 93774 Columns: 16
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr   (2): location, iso_code
## dbl  (13): total_vaccinations, people_vaccinated, people_fully_vaccinated, t...
## date  (1): date
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

# 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.

# import iso-codes of Europian countries from csvfile

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

## Rows: 28 Columns: 3
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (3): Country, Alpha-2, Alpha-3
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

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()

European Countries in the vaccinations tibble
Austria
Belgium
Bulgaria
Croatia
Cyprus
Czechia
Denmark
Estonia
Finland
France
Germany
Greece
Hungary
Ireland
Italy
Latvia
Lithuania
Luxembourg
Malta
Netherlands
Poland
Portugal
Romania
Slovakia
Slovenia
Spain
Sweden
United Kingdom

# 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)

## Warning: Removed 30 row(s) containing missing values (geom_path).

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

cat('Data structure:\n\n')

## Data structure:

glimpse(df_vacc)

## Rows: 93,774
## Columns: 16
## $ location                            <chr> "Afghanistan", "Afghanistan", "Afg~
## $ iso_code                            <chr> "AFG", "AFG", "AFG", "AFG", "AFG",~
## $ date                                <date> 2021-02-22, 2021-02-23, 2021-02-2~
## $ total_vaccinations                  <dbl> 0, NA, NA, NA, NA, NA, 8200, NA, N~
## $ people_vaccinated                   <dbl> 0, NA, NA, NA, NA, NA, 8200, NA, N~
## $ people_fully_vaccinated             <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA~
## $ total_boosters                      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA~
## $ daily_vaccinations_raw              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA~
## $ daily_vaccinations                  <dbl> NA, 1367, 1367, 1367, 1367, 1367, ~
## $ total_vaccinations_per_hundred      <dbl> 0.00, NA, NA, NA, NA, NA, 0.02, NA~
## $ people_vaccinated_per_hundred       <dbl> 0.00, NA, NA, NA, NA, NA, 0.02, NA~
## $ people_fully_vaccinated_per_hundred <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA~
## $ total_boosters_per_hundred          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA~
## $ daily_vaccinations_per_million      <dbl> NA, 34, 34, 34, 34, 34, 34, 40, 45~
## $ daily_people_vaccinated             <dbl> NA, 1367, 1367, 1367, 1367, 1367, ~
## $ daily_people_vaccinated_per_hundred <dbl> NA, 0.003, 0.003, 0.003, 0.003, 0.~

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

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.

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.

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.

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

Exercises Laboratory Session 01

Nicola Zomer

08 April, 2022