reproduce-analysis.qmd

---
title: "Reproduce analysis for the paper 'Measuring gig economy in Poland using mobile big data'"
author: Mciej Beręsewicz
format: 
  html:
    self-contained: true
    table-of-contents: true
    number-sections: true
    df-print: kable
editor: visual
execute: 
  eval: true
  warning: false
  message: false
---

# Packages and session info

```{r eval = F}
install.packages(c("tidyverse", "scales", "ggh4x", "ggrepel", "xtable"))
```

```{r}
library(tidyverse)
library(scales)
library(ggh4x)
library(ggrepel)
library(xtable)
theme_set(theme_bw())
```

R version

```
R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 11.7.5
```

# Read the data

Data purchased from Selectivv

+ 1 -- monthly statistics
+ 2 -- characteristics
+ 3 -- place of residence
+ 4 - app activity

```{r}
selectivv <- list()
selectivv[[1]] <- read_csv(file = "data-raw/gig-table1-monthly-counts-stats.csv", show_col_types = F)
selectivv[[2]] <- read_csv(file = "data-raw/gig-table2-halfyear-demo-stats.csv", show_col_types = F)
selectivv[[3]] <- read_csv(file = "data-raw/gig-table3-halfyear-region-stats.csv", show_col_types = F)
selectivv[[4]] <- read_csv(file = "data-raw/gig-table4-halfyear-activity-stats.csv", show_col_types = F)
```

LFS Data
```{r}
bael <- read_csv("data-raw/bael.csv", show_col_types = F) |>
  na.omit() |>
  unite(period, year:yearhalf, sep = ".") |>
  rename(region_name = nazwa_miasta)
```

Indices 

```{r}
selectivv_ind <- read_csv("data-raw/selectivv-indices.csv", show_col_types = F)
selectivv_ind_with_zus <- read_csv("data-raw/selectivv-indices-comp-zus.csv", show_col_types = F)
```

Validation data from registers and Selectivv

```{r}
valid_pesel_polish <- read_csv("data-raw/validation-pesel-register-polish.csv", show_col_types = FALSE)
valid_pesel_ukr <- read_csv("data-raw/validation-pesel-register-ukr.csv", show_col_types = FALSE)
valid_sel_polish <- read_csv("data-raw/validation-selectivv-polish.csv", show_col_types = FALSE)
valid_sel_ukr <- read_csv("data-raw/validation-selectivv-ukr.csv", show_col_types = FALSE)
valid_zus_ukr <- read_csv("data-raw/validation-zus-register-ukr.csv", show_col_types = FALSE)
```

Validation survey 

```{r}
valid_survey <- read_csv("data-raw/validation-survey-1000.csv", show_col_types = F)
head(valid_survey)
```


Self-reported values

+ q1_a -- Poland
+ q1_b -- Ukraine
+ q1_c -- Belarus
+ q2_a -- Female
+ q2_b -- Male
+ q3_a -- 18-24 age
+ q3_b -- 25-29 age
+ q3_c -- 30-39 age
+ q3_d -- 40+ age

Selectivv classifications

+ gender_women -- Female
+ gender_men -- Male
+ age_18_24 -- 18-24 age
+ age_25_29 -- 25-29 age
+ age_30_39 -- 30-39 age
+ age_40 -- 40+ age
+ country_PL -- Poland
+ country_UA -- Ukraine
+ country_BY -- Belarus


# Reproduce results

## Section 'Assessment of the coverage error'

Figure 1 -- Changes in the number of smartphone users based on monthly data collected by Selectivv between 2018 and 2020 (2018.01 = 100)

```{r}
selectivv_ind |>
  mutate(group = factor(group, c("total", "polish", "ukrainian", "other"),
                        c("All", "Polish", "Ukrainian", "Other foreigners")),
         label = if_else(month == max(month), paste0(as.character(group), " (", round(index), "%)"),  NA_character_),
         label = ifelse(group == "Ukrainian" & month == max(month), 
                        paste0(as.character(group), "\n(", round(index), "%)"), label)) |> 
  ggplot(data = _, aes(x = month, y = index, group = group, color = group)) + 
  geom_line(show.legend = FALSE) +
  scale_colour_brewer(type = "qual", palette = "Dark2", name = "Users")  +
  geom_text_repel(aes(label = label),
                  nudge_x = -0.25, 
                  direction = "y", 
                  hjust = "bottom",
                  box.padding = .2,
                  force = 1,
                  na.rm = TRUE,
                  show.legend = FALSE) +
  labs(x = "Month", y = "Index (2018.01=100)") -> p_cov

p_cov

ggsave(plot = p_cov, file = "figures//fig-sel-index-total.png", width = 7, height = 4)

```

Figure 2 -- A~comparison of changes in the number of users according to Selectivv and the number of insured foreigners according to ZUS (2019.01 = 100). Note: ZUS (M) -- monthly statistics; ZUS (Q) -- quarterly statistics


```{r}
selectivv_ind_with_zus  |>
  mutate(group = factor(group, c("ukrainian", "other"), c("Ukraine", "Other")),
         source = factor(source, c("selectivv", "zus", "zus (kw)"), c("Selectivv", "ZUS (M)", "ZUS (Q)")),
         label = if_else(month == max(month), as.character(source), NA_character_)) |>
  ggplot(data = _, aes(x = month, y = index_new, color = source, group = source)) +
  geom_line(show.legend = FALSE) + 
  geom_text_repel(aes(label = label),
                  nudge_x = -0.25, 
                  direction = "y", 
                  hjust = "left",
                  box.padding = .1,
                  force = .8,
                  na.rm = TRUE,
                  show.legend = FALSE) + 
  facet_wrap(~group)  +
  theme_bw() +
  scale_color_brewer(type = "qual", palette = "Dark2") + 
  labs(x = "Month", y = "Index", linetype = "Source") +
  scale_x_date(limits = c(as.Date("2018-01-01"), as.Date("2021-09-01"))) + 
  theme(strip.text = element_text(size = 15)) -> p1

p1

ggsave(plot = p1, file = "figures/fig-sel-zus.png", width = 7, height = 4)

```

## Section 'Assessment of the measurement error'

Figure 3 -- Comparison of gender and age of Polish citizens according to the Population Register (PESEL) and Selectivv

```{r}
valid_pesel_polish |> 
  add_count(period, wt = pesel, name = "total") |>
  mutate(p_pesel = pesel / total) |>
  select(-pesel, -total) |>
  left_join(
    valid_sel_polish
  ) |>
  rename(pesel=p_pesel, selectivv=p_sel) |> 
  pivot_longer(cols = pesel:selectivv, names_to = "source", values_to = "share") |>
  mutate(gender = factor(gender, c("Females", "Males")),
         source = factor(source, c("pesel", "selectivv"), c("Population\nregister", "Selectivv"))) |>
  ggplot(aes(x = period, y = share, 
             linetype = source,  
             group = interaction(source, gender),  
             color = gender)) +
  geom_line() +
  #geom_point() + 
  facet_wrap(~age, ncol = 4) + 
  scale_y_continuous(limits = c(0, .31), labels = scales::percent, n.breaks = 10) +
  scale_color_brewer(type = "qual", palette = "Set1") + 
  labs(x = "Period", y = "Share", linetype = "Source", color = "Gender") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) -> p3

p3

ggsave(plot = p3, filename = "figures/comp-polish-gender-age.png", width = 6, height = 4)
```


Figure 4 -- Comparison of gender and age of Ukrainian citizens according to the Population Register (PESEL), Social Insurance Register (ZUS) and Selectivv

```{r}
valid_pesel_ukr |> 
  add_count(period, wt = pesel, name = "total") |> 
  mutate(p_pesel = pesel / total) |>
  select(-pesel, -total) |>
  left_join(
    valid_sel_ukr
  ) |>
  left_join(
    valid_zus_ukr |>
      add_count(period, wt = zus, name = "total") |>
      mutate(p_zus = zus/total) |>
      select(-zus, -total) 
  ) |>
  rename(pesel=p_pesel, selectivv=p_sel, zus=p_zus) |> 
  pivot_longer(cols = pesel:zus, names_to = "source", values_to = "share") |>
  mutate(gender = factor(gender, c("Females", "Males")),
         source = factor(source, c("pesel", "zus", "selectivv"), 
                         c("Population\nregister", "ZUS", "Selectivv"))) |>
  ggplot(aes(x = period, y = share, 
             linetype = source,  
             group = interaction(source, gender),  
             color = gender)) +
  geom_line() +
  #geom_point() + 
  facet_wrap(~age, ncol = 4) + 
  scale_y_continuous(limits = c(0, .31), labels = scales::percent, n.breaks = 10) +
  scale_linetype_manual(values = c(1,3,2)) + 
  scale_color_brewer(type = "qual", palette = "Set1") + 
  labs(x = "Period", y = "Share", linetype = "Source", color = "Gender") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) -> p4

p4
ggsave(plot = p4, filename = "figures/comp-ukr-gender-age.png", width = 6, height = 4)

```

Tables: 

+ Classification table for nationality

```{r}
valid_survey |>
  mutate(country_stated = case_when(q1_a == 1 ~ "Poland",
                                    q1_b == 1 ~ "Ukraine",
                                    q1_c == 1 ~ "Belarus"),
         country_class = case_when(country_PL == 1 ~ "Poland",
                                   country_UA == 1 ~ "Ukraine",
                                   country_BY == 1 ~ "Belarus"),
         country_stated = factor(country_stated, c("Poland", "Ukraine", "Belarus")),
         country_class = factor(country_class, c("Poland", "Ukraine", "Belarus"))) |>
  xtabs(formula=~country_class + country_stated) 
```

+ Classification table for gender

```{r}
valid_survey |>
  mutate(gender_stated = case_when(q2_a == 1 ~ "Female",
                                    q2_b == 1 ~ "Male"),
         gender_class = case_when(gender_women == 1 ~ "Female",
                                   gender_men == 1 ~ "Male")) |>
  xtabs(formula=~gender_class + gender_stated) 
```

+ Classification table for age groups


```{r}
valid_survey |>
  mutate(age_stated = case_when(q3_a == 1 ~ "18-24",
                                q3_b == 1 ~ "25-29",
                                q3_c == 1 ~ "30-39",
                                q3_d == 1 ~ "40+"),
         age_class = case_when(age_18_24 == 1 ~ "18-24",
                               age_25_29 == 1 ~ "25-29",
                               age_30_39 == 1 ~ "30-39",
                               age_40 == 1 ~ "40+")) |>
  xtabs(formula=~age_class + age_stated) 
```
## Section 'Characteristics of Polish gig workers'

Figure -- The monthly number of active users in Poland for selected apps by category between 2018 and 2020

```{r}
selectivv[[1]] |>
  mutate(label = if_else(month == max(month), as.character(apps), NA_character_),
         category = ifelse(apps == "Uber", "Uber & Uber Eats", category),
         apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = factor(category, c("Transportation", "Delivery", "Uber & Uber Eats"))) |>
  ggplot(data = _, aes(x = month, y = number_of_users/1000, color = apps, group = apps)) +
  geom_line(show.legend = FALSE) + 
  geom_text_repel(aes(label = label),
                  nudge_x = -0.25, 
                  direction = "y", 
                  hjust = "left",
                  box.padding = .1,
                  force = .8,
                  na.rm = TRUE,
                  show.legend = FALSE) + 
  scale_color_brewer(type = "qual", palette = "Dark2") + 
  labs(x = "The reference period", y = 'Number of users (thousands)', color = "App") + 
  facet_wrap(~category, scales = "free_y") +
  scale_x_date(limits = c(as.Date("2018-01-01"), as.Date("2021-09-01"))) + 
  theme(strip.text = element_text(size = 15)) ->  p2

p2

ggsave(plot = p2, file = "figures/fig-sel-month.png", width = 6, height = 4)

```

Figure -- Change in the number of active users in Poland for selected apps by type between 2018 and 2020. Base 2018HY1


```{r}
selectivv[[3]] |>
  count(category, apps, group, period, wt = number_of_users) |>
  spread(group, n) |>
  mutate(fua = `Functional area` - Cities, 
         woj = Voivodeship - `Functional area`) |>
  select(category, apps, period, cities = Cities, fua, woj) |>
  gather(geo, val, -apps, -period, -category) |>
  mutate(geo = factor(geo, c("cities", "fua", "woj"),
                      c("Cities", "Functional UAs", "Provinces")),
         category = ifelse(apps == "Uber", "Uber & Uber Eats", category),
         category  = factor(category, c("Transportation", "Delivery", "Uber & Uber Eats"))) |>
  group_by(apps, geo) |>
  mutate(index = val/val[1],
         diff = val-val[1],
         apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier"))) |>
  ggplot(data = _, aes(x = as.factor(period), y = diff/1000, group = apps, color = apps)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_grid(category ~ geo, scales = "free_y") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "The reference period (t0=2018.1)", 
       y = "Change in the number of users\n(in thousands)", color = "Apps") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  guides(color = guide_legend(nrow = 2, byrow = TRUE)) ->  p1b

p1b

ggsave(plot = p1b, file = "figures/fig-sel-diff-hy.png", width = 8, height = 5)

```

Figure -- The monthly number of active users of selected apps in Poland by country, age and gender between 2018 and 2020

```{r}
selectivv[[2]] |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, country, period, wt = n_unique_users) |>  
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = country, color = country)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Nationality") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2a

p2a

ggsave(plot = p2a, file = "figures/fig-sel-half-country.png", width = 8, height = 4)
```

```{r}
selectivv[[2]]  |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, age, period, wt = n_unique_users) |>
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = age, color = age)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Age group") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2b

p2b

ggsave(plot = p2b, file = "figures/fig-sel-half-age.png", width = 8, height = 4)


```


```{r}
selectivv[[2]]  |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, gender, period, wt = n_unique_users)  |>
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = gender, color = gender)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Sex") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2c

p2c

ggsave(plot = p2c, file = "figures/fig-sel-half-sex.png", width = 8, height = 4)

```


Figure -- Point and 95\% interval of ratio of the number of active users and the LFS working 18--65 population by app for 2020HY1. Note that Y axis differs between the apps

```{r}
selectivv[[3]]  |>
    filter(group == "Cities")  |>
    mutate(period = as.character(period)) |>
    select(category, period, region_name, apps, number_of_users) |>
    left_join(
      bael |>
        mutate(region_name = str_remove(region_name, "M. "))
    ) |>
    filter(number_of_users > 0 ) -> cities_bael
```


```{r}
cities_bael |>
  filter(region_name != "Sopot", period == "2020.1") |>
  filter(!is.na(y))  |>
  mutate(share = number_of_users/y*100,
         share_hi = number_of_users / (y - 1.96*rmse)*100,
         share_lo = number_of_users / (y + 1.96*rmse)*100) |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier"))) |>
  ggplot(data = _, aes(x = reorder(region_name, -y), y = share, ymin = share_lo, ymax = share_hi)) +
  geom_errorbar() +
  geom_point()  +
  facet_wrap(~apps, ncol = 4, scale = "free_y") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 5))  +
  #scale_y_continuous(labels = scales::label_percent(accuracy = 2)) +
  labs(x = "Cities (ordered by the size of the working population according to the LFS)", 
       y = "Number of users / LFS working population [%]") -> p3

p3

ggsave(plot = p3, file = "figures/fig-sel-cities-shares.png", width = 8, height = 4)

```
Figure -- Average and standard deviation (Std. Dev.) of activity time by app type, app, weekdays and hours


```{r}
selectivv[[4]] |>
  mutate(category = ifelse(apps == "Uber", "Uber & Uber Eats", category),
         category  = factor(category, c( "Uber & Uber Eats", "Transportation", "Delivery")),
         apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         day = factor(day, c("Mondays-Thursdays", "Fridays-Sundays")),
         hour = str_replace(hour, "day", "daytime"),
         hour = str_replace(hour, "night", "nightime"))  |>
  ggplot(data = _, aes(x = as.character(period), y = activity_time, color = apps, group = apps)) +
  geom_line() +
  geom_point() + 
  facet_nested(day ~ category +  hour + statistic)  +
  scale_color_brewer(type = "qual", palette = "Dark2", name = "Apps") +  
  theme(axis.text.x=element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") +
  labs(x = "Half-year", y = "Time spent in the app") + 
  scale_y_continuous(breaks = seq(0,7,1)) +
  guides(color = guide_legend(nrow = 2, byrow = TRUE))  -> p4

p4

ggsave(plot = p4, file = "figures/fig-sel-half-time-dev.png", width = 8, height = 5)

```

# Supplementary materials


Table  Number of active users in cities for 2020HY2

```{r}
selectivv[[3]] |>
    filter(group == "Cities") |>
    mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                                 "Takeaway", "Wolt", "Glover",  "Bolt Courier"))) |>
    filter(number_of_users > 0, period == "2020.2")  |> 
    select(period, apps, region_name, number_of_users) |>
    spread(apps, number_of_users) |> 
    select(-period) |> 
    xtable(digits = 0) |> 
    print.xtable(include.rownames = F, NA.string = "--")
```


Additional plots for characteristics

```{r}
selectivv[[2]]  |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, gender, period, wt = n_unique_users*parents_of_children_0_4_years)  |>
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = gender, color = gender)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  scale_y_continuous(limits = c(0,2)) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Sex") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2d

p2d

ggsave(plot = p2d, file = "figures/fig-sel-half-parents04.png", width = 8, height = 4)

```
```{r}
selectivv[[2]]  |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, gender, period, wt = n_unique_users*parents_of_children_5_10_years)  |>
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = gender, color = gender)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  scale_y_continuous(limits = c(0,2)) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Sex") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2e

p2e

ggsave(plot = p2e, file = "figures/fig-sel-half-parents510.png", width = 8, height = 4)

```

```{r}
selectivv[[2]]  |>
  mutate(apps = factor(apps, c("Uber", "Bolt Driver", "iTaxi", "FREE NOW",
                               "Takeaway", "Wolt", "Glover",  "Bolt Courier")),
         category = ifelse(apps == "Uber", "Uber &\n Uber Eats", category),
         category  = factor(category, c("Uber &\n Uber Eats", "Transportation", "Delivery"))) |>
  count(category, apps, gender, period, wt = n_unique_users*students)  |>
  ggplot(data = _, aes(x = as.character(period), y = n/1000, group = gender, color = gender)) +
  geom_line() +
  geom_point() + 
  scale_color_brewer(type = "qual", palette = "Dark2") +  
  facet_nested(~category + apps) + 
  scale_y_continuous(limits = c(0,2)) + 
  labs(x = "Half-year", y = "Number of users (in thousands)", color = "Sex") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "bottom") -> p2f

p2f

ggsave(plot = p2f, file = "figures/fig-sel-half-students.png", width = 8, height = 4)

```



Table  -- Direct estimates, its standard errors (in thousands) and sample sizes of the working population size (18--64) based on the Labour Force Survey between 2018HY1 and 2020HY

```{r}
bael |>
  filter(period == "2018.1")  |> 
  mutate(region_name = str_remove(region_name, "M. "),
         y = y / 1000,
         rmse = rmse / 1000,
         cv = rmse/y*100) |>
  select(region_name, size = y, std = rmse, sample = n) |>
  left_join(
    bael |>
    filter(period == "2018.2")  |> 
    mutate(region_name = str_remove(region_name, "M. "),
         y = y / 1000,
         rmse = rmse / 1000,
         cv = rmse/y*100) |>
    select(region_name, size = y, std = rmse, sample = n),
    by = "region_name"
  ) |>
  left_join(
    bael |>
    filter(period == "2019.1")  |> 
    mutate(region_name = str_remove(region_name, "M. "),
         y = y / 1000,
         rmse = rmse / 1000,
         cv = rmse/y*100) |>
    select(region_name, size = y, std = rmse, sample = n),
    by = "region_name"
  )  |>
  left_join(
    bael |>
    filter(period == "2019.2")  |> 
    mutate(region_name = str_remove(region_name, "M. "),
         y = y / 1000,
         rmse = rmse / 1000,
         cv = rmse/y*100) |>
    select(region_name, size = y, std = rmse, sample = n),
    by = "region_name"
  )  |>
  left_join(
    bael |>
    filter(period == "2020.1")  |> 
    mutate(region_name = str_remove(region_name, "M. "),
         y = y / 1000,
         rmse = rmse / 1000,
         cv = rmse/y*100) |>
    select(region_name, size = y, std = rmse, sample = n),
    by = "region_name"
  )  |>
  mutate_at(vars(starts_with("sample")), as.integer) |>
  xtable(digits = 1) |>
  print.xtable(include.rownames = F)
  
```