CTIS_Analysis_Mental_Health.Rmd

---
title: "CTIS Mental Health"
author: "Christian Hobelsberger, Max Lang, Elisa Shima"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# devtools::install_github("CaroHaensch/CTIS")
library(CTIS)
library(tidyverse)
library(timetk)
library(lubridate)
library(reshape2)
library(corrplot)
library(ozmaps)
library(RColorBrewer)
library(countrycode)
library(gganimate)
library(gifski)
library(rgdal)
library(tmap)
library(maps)
library(sf)
library(tigris)
library(rworldmap)
library(xts)
library(lubridate)
library(lares)
library(highcharter)
library(dygraphs)
library(GGally)


```


## Data Import & Wrangling

We are using the R Package CTIS to get access to the API Interface, more information on <https://github.com/CaroHaensch/CTIS>.

```{r Data Import & Wrangling}
indicator_list = list("anxious_7d", "depressed_7d", "worried_become_ill", "finance", "food_security")

data <- lapply(X = indicator_list, 
               FUN = function(x){
                 CTIS_open_data_country(indicator = x, type = "daily", country = "all", daterange = "")
                 })
data <- as.data.frame(do.call(rbind, data))
data_trimmed <- data[,c(1:5, 7)]
data_aggregated <- dcast(data = data_trimmed, 
                         formula = data.country+data.iso_code+data.gid_0+data.survey_date~data.indicator, 
                         value.var = "data.pct")

data_aggregated$data.country <- as.factor(data_aggregated$data.country)
data_aggregated$data.iso_code <- as.factor(data_aggregated$data.iso_code)
data_aggregated$data.gid_0 <- as.factor(data_aggregated$data.gid_0)

data_aggregated$data.survey_date <- as.Date(x = data_aggregated$data.survey_date, "%Y%m%d")



write.csv(data_aggregated,"data/data_aggregated.csv", row.names = FALSE)
data_CTIS = read.csv(file = "data/data_aggregated.csv")
data_CTIS$data.country <- as.factor(data_CTIS$data.country)
data_CTIS$data.iso_code <- as.factor(data_CTIS$data.iso_code)
data_CTIS$data.gid_0 <- as.factor(data_CTIS$data.gid_0)

data_CTIS$data.survey_date <- as.Date(x = data_CTIS$data.survey_date)

```

# Read our world in data (<https://ourworldindata.org/policy-responses-covid>)

```{r Policy Data Import}
policy_income_support <- read.csv("data/ourworldindata-policy-responses/income-support-covid.csv")
policy_stay_at_home <- read.csv("data/ourworldindata-policy-responses/stay-at-home-covid.csv")
policy_changes_visitors <- read.csv("data/ourworldindata-policy-responses/changes-visitors-covid.csv")
policy_school_closures <- read.csv("data/ourworldindata-policy-responses/school-closures-covid.csv")

names(policy_income_support)[1:3] <- c("data.country", "data.iso_code", "data.date")
names(policy_stay_at_home)[1:3] <- c("data.country", "data.iso_code", "data.date")
names(policy_changes_visitors)[1:3] <-  c("data.country", "data.iso_code", "data.date")
names(policy_school_closures)[1:3] <- c("data.country", "data.iso_code", "data.date")

policy_data <- policy_stay_at_home %>% 
  left_join(policy_income_support, by = c("data.country", "data.iso_code", "data.date")) %>% 
  left_join(policy_changes_visitors, by = c("data.country", "data.iso_code", "data.date")) %>% 
  left_join(policy_school_closures, by = c("data.country", "data.iso_code", "data.date"))

policy_data$data.date <- as.Date(policy_data$data.date)
# Scale from 0 to 3 restrictions
policy_data$stay_home_requirements <- factor(policy_data$stay_home_requirements, ordered = TRUE) 
policy_data$income_support <- factor(policy_data$income_support, ordered = TRUE)
policy_data$school_closures<- factor(policy_data$school_closures, ordered = TRUE)
policy_data$data.country<- factor(policy_data$data.country)
policy_data$data.iso_code<- factor(policy_data$data.iso_code)

saveRDS(object = policy_data, file = "data/policy_data.RDS")
```

# Large Data with policy

```{r}
data_CTIS_policy <- data_CTIS %>% 
  left_join(policy_data, by = c("data.country", "data.iso_code", "data.survey_date" = "data.date"))
# Adding continent
data_CTIS_policy$continent <- countrycode(data_CTIS_policy$data.country, origin = "country.name", destination = "continent")

saveRDS(object = data_CTIS_policy, file = "data/data_CTIS_policy.RDS")
```


# Create mapping data
```{r}
data("World")
world <- World[c("iso_a3", "sovereignt", "geometry")]
colnames(world) <- c("data.iso_code", "data.country", "geometry")

# WGS 84 standard
data_CTIS_map <- world %>%
  full_join(data_CTIS_policy, by = c("data.iso_code", "data.country"))

#data_CTIS_map <- data_CTIS_map[!st_is_empty(data_CTIS_map),,drop=FALSE]

data_CTIS_map <-  data_CTIS_map %>% 
filter( is.na(st_dimension(.)) == FALSE )

# Match Countrycode again because now new countries that were not in the survey but now joined
data_CTIS_map$continent <- countrycode(data_CTIS_map$data.country, origin = "country.name", destination = "continent")

saveRDS(data_CTIS_map, file = "shiny/ctis-app/app-data/data_CTIS_map.RDS")
```

## Chris

### First look at data

Brief first look at the CTIS data.

```{r Chris: Data glimpse}
# Read data
data_CTIS = read.csv(file = "data/data_aggregated.csv")
data_CTIS$data.country <- as.factor(data_CTIS$data.country)
data_CTIS$data.iso_code <- as.factor(data_CTIS$data.iso_code)
data_CTIS$data.gid_0 <- as.factor(data_CTIS$data.gid_0)

data_CTIS$data.survey_date <- as.Date(x = data_CTIS$data.survey_date)

# Number of NAs for each row
colSums(is.na(data_CTIS))

# Number of countries with equal or more than 600 observations
sum(table(data_CTIS$data.country) >= 600)

# Names of the countries with equal or more than 600 observations
names_sufficient_data = names(which(table(data_CTIS$data.country) >= 600))
names_sufficient_data

# Get the 10 countries with the highest amount of observations
names(sort(table(data_CTIS$data.country), decreasing = TRUE)[1:10])

head(data_CTIS)

summary(data_CTIS)

# Get mean anxious_7d for each country
data_CTIS %>% group_by(data.country) %>% summarise(mean(anxious_7d, na.rm = TRUE))

# Arrange by mean_anxious_7d of each country
mean_anxious_7d_arranged <- data_CTIS %>% 
  filter(data.country %in% names_sufficient_data) %>%
  group_by(data.country) %>% 
  summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  arrange(mean_anxious_7d)
mean_anxious_7d_arranged

# Arrange by median_anxious_7d of each country
median_anxious_7d_arranged <- data_CTIS %>% 
  filter(data.country %in% names_sufficient_data) %>%
  group_by(data.country) %>% 
  summarise(median_anxious_7d = median(anxious_7d, na.rm = TRUE)) %>% 
  arrange(median_anxious_7d)

# Calculating pearson correlation between anxious_7d and depressed_7d: 0.8064084
drop_na(data_CTIS[, 5:6]) |> cor()

# Plot: Denisity of mean_anxious_7d
ggplot(data = mean_anxious_7d_arranged, mapping = aes(y = mean_anxious_7d)) + geom_density() + 
  coord_flip() +
  labs(title = "Denisity of mean_anxious_7d")
ggsave("graphics/density_mean_anxious_7d.png", width = 12, height = 9)

# First set of plot
country_subset1 <- c("Germany", "Mexico", "South Africa")
data_CTIS_country_subset1 = data_CTIS %>%
  filter(data.country %in% country_subset1)

ggplot(data = data_CTIS_country_subset1, 
       mapping = aes(x = data.survey_date, y = anxious_7d, color = data.country)) +
  geom_line() +
  labs(title = "anxious_7d vs. date (1)")
ggsave("graphics/anxious_7d_set1.png", width = 12, height = 9) 

ggplot(data = data_CTIS_country_subset1, 
       mapping = aes(x = data.survey_date, y = depressed_7d, color = data.country)) +
  geom_line() +
  labs(title = "depressed_7d vs. date (1)")
ggsave("graphics/depressed_7d_set1.png", width = 12, height = 9)

# Second set of plot
country_subset2 = c("Denmark", "India", "Turkey")
data_CTIS_country_subset2 = data_CTIS %>%
  filter(data.country %in% country_subset2)

ggplot(data = data_CTIS_country_subset2, 
       mapping = aes(x = data.survey_date, y = anxious_7d, color = data.country)) +
  geom_line() +
  labs(title = "anxious_7d vs. date (2)")
ggsave("graphics/anxious_7d_set2.png", width = 12, height = 9)

ggplot(data = data_CTIS_country_subset2, 
       mapping = aes(x = data.survey_date, y = depressed_7d, color = data.country)) +
  geom_line() +
  labs(title = "depressed_7d vs. date (2)")
ggsave("graphics/depressed_7d_set2.png", width = 12, height = 9)

# Third set of plot
country_subset3 = c("Norway", "Hungary", "Egypt")
data_CTIS_country_subset3 = data_CTIS %>%
  filter(data.country %in% country_subset3)

ggplot(data = data_CTIS_country_subset3, 
       mapping = aes(x = data.survey_date, y = anxious_7d, color = data.country)) +
  geom_line() +
  labs(title = "anxious_7d vs. date (3)")
ggsave("graphics/anxious_7d_set3.png", width = 12, height = 9)

ggplot(data = data_CTIS_country_subset3, 
       mapping = aes(x = data.survey_date, y = depressed_7d, color = data.country)) +
  geom_line() +
  labs(title = "depressed_7d vs. date (3)")
ggsave("graphics/depressed_7d_set3.png", width = 12, height = 9)

# Boxplots: Total subset
total_subset <- c(country_subset1, country_subset2, country_subset3)
data_CTIS_total_subset = data_CTIS %>%
  filter(data.country %in% total_subset)

ggplot(data = data_CTIS_total_subset, 
       mapping = aes(y = anxious_7d, color = data.country)) +
  geom_boxplot() +
  labs(title = "anxious_7d Boxplot")
ggsave("graphics/anxious_7d_boxplot.png", width = 12, height = 9)

ggplot(data = data_CTIS_total_subset, 
       mapping = aes(y = depressed_7d, color = data.country)) +
  geom_boxplot() +
  labs(title = "depressed_7d Boxplot")
ggsave("graphics/depressed_7d_boxplot.png", width = 12, height = 9)

```




# Max

## Data preparation

```{r Max: Data preparation}
country_subset_large <- c("Germany", "Mexico", "South Africa", "Norway", "Hungary",
                          "Egypt", "Denmark", "India", "Turkey")
data_CTIS_country_subset_large = data_CTIS %>%
  filter(data.country %in% country_subset_large)


data_policy_subset_large = policy_data %>% 
  filter(data.country %in% country_subset_large)

data_CTIS_policy_large <- data_CTIS_country_subset_large %>% 
  left_join(data_policy_subset_large, by = c("data.country",
                                             "data.iso_code",
                                              "data.survey_date" = "data.date"))
```

-   Only turkey and india have level 3 stay home requirements

```{r Check for hard restricting countries}
hard_countries <- unique(data_CTIS_policy_large %>% 
  filter(stay_home_requirements %in% 2:3 & school_closures %in% 2:3) %>% dplyr::select(data.country))[["data.country"]]
data_CTIS_policy_large %>% 
  filter(stay_home_requirements %in% 3) %>% dplyr::select(data.country) %>%  unique()
```

## Anxiety

```{r Anxiety over time for country selection}
anx_hard_countries <- data_CTIS_policy_large %>% filter(data.country %in% hard_countries) %>% dplyr::select(anxious_7d)

data_CTIS_policy_large %>% 
  filter(stay_home_requirements %in% 2:3 & school_closures %in% 2:3) %>% ggplot(aes(x= data.survey_date, y = anxious_7d, color = data.country)) + 
  geom_line() + 
  geom_hline(yintercept = mean(anx_hard_countries$anxious_7d, na.rm = T)) + 
  geom_hline(yintercept = mean(data_CTIS_policy_large$anxious_7d), color = "blue")+
  scale_color_viridis_d(option = "A")
```




```{r school closure in 3}
hard_countries <- unique(data_CTIS_policy_large %>% 
  filter(school_closures %in% 3) %>% dplyr::select(data.country))[["data.country"]]

anx_hard_countries <- data_CTIS_policy_large %>% filter(data.country %in% hard_countries) %>% dplyr::select(anxious_7d)

data_CTIS_policy_large %>% 
  filter(school_closures %in% 3) %>% ggplot(aes(x= data.survey_date, y = anxious_7d, color = data.country)) + 
  geom_line() + 
  geom_hline(yintercept = mean(anx_hard_countries$anxious_7d, na.rm = T)) + 
  geom_hline(yintercept = mean(data_CTIS_policy_large$anxious_7d), color = "blue")+
  scale_color_viridis_d(option = "A")
```

## Depression

```{r Depression over time by country}
hard_countries <- unique(data_CTIS_policy_large %>% 
  filter(stay_home_requirements %in% 2:3 & school_closures %in% 2:3) %>% dplyr::select(data.country))[["data.country"]]

anx_hard_countries <- data_CTIS_policy_large %>% filter(data.country %in% hard_countries) %>% dplyr::select(depressed_7d)

data_CTIS_policy_large %>% 
  filter(stay_home_requirements %in% 2:3 & school_closures %in% 2:3) %>% ggplot(aes(x= data.survey_date, y = depressed_7d, color = data.country)) + 
  geom_line() + 
  geom_hline(yintercept = mean(anx_hard_countries$depressed_7d, na.rm = T)) + 
  geom_hline(yintercept = mean(data_CTIS_policy_large$depressed_7d), color = "blue")+
  scale_color_viridis_d(option = "A")
```

## Correlation Matrix

```{r corrplot for subset}
numeric_data <- data_CTIS_policy_large %>% dplyr::select(anxious_7d, depressed_7d, finance, food_security, worried_become_ill, parks, workplaces, residential, transit_stations)
cormat <- cor(numeric_data, use = "complete.obs")

corrplot(cormat, method = "color")
```

# Scatterplot stay home park movement and anxiety

```{r}
data_CTIS_policy_large %>% 
  #filter(data.survey_date %in% as.Date("2021-03-01")) %>% 
  ggplot(aes(x = parks, y = anxious_7d, color = stay_home_requirements))+
  facet_wrap(facets = vars(data.country))+
  geom_point(alpha = 0.2)
  
```

# Boxplots Stay Home requirements / school closure influence on mental health

```{r}
data_CTIS_policy_large %>% 
  ggplot(aes(x = stay_home_requirements, y = anxious_7d))+
  facet_wrap(vars(data.country))+
  geom_boxplot()
```

```{r}
data_CTIS_policy_large %>% 
  ggplot(aes(x = school_closures, y = anxious_7d))+
  facet_wrap(vars(data.country))+
  geom_boxplot()+
  geom_hline(yintercept = mean(data_CTIS$anxious_7d, na.rm = TRUE))
```

```{r}
data_CTIS_policy_large %>% 
  ggplot(aes(x = stay_home_requirements, y = finance))+
  facet_wrap(vars(data.country))+
  geom_boxplot()+
  geom_hline(yintercept = mean(data_CTIS$finance, na.rm = TRUE), linetype = "dotted")+
  ggtitle("Influence of different levels of stay home requirements on household’s finances")+
  labs(x = "Level of stay home requirements",
       y = "Relative amount of users worried about finance",
       subtitle = "Dotted line shows mean of 'finance' variable over the complete data set")

ggsave("plots/Max/box_sub_stayhome_finance.png", width = 12, height = 9)
```

```{r}
data_CTIS_policy_large %>% 
  ggplot(aes(x = stay_home_requirements, y = finance))+
  facet_wrap(vars(data.country))+
  geom_boxplot()+
  geom_hline(yintercept = mean(data_CTIS$finance, na.rm = TRUE))
```



# Using complete data set

```{r}
numeric_data <- data_CTIS_policy %>% dplyr::select(anxious_7d, depressed_7d, finance, food_security, worried_become_ill, parks, workplaces, residential, transit_stations)
cormat <- cor(numeric_data, use = "complete.obs")

corrplot(cormat, method = "color")
```

# Continent grouped

```{r Anxiety continent level}
data_CTIS_policy %>% 
  group_by(continent, data.survey_date) %>% 
  summarise(continent_mean_anxious7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  ggplot(aes(x = data.survey_date, y = continent_mean_anxious7d, color = continent))+
    geom_line()+
  scale_color_viridis_d()
```

```{r Anxiety continent level}
data_CTIS_policy %>% 
  group_by(continent, data.survey_date) %>% 
  summarise(continent_mean_anxious7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  ggplot(aes(x = data.survey_date, y = continent_mean_anxious7d))+
  facet_wrap(vars(continent))+
    geom_line()+
  labs(x = "Survey date",
       y = "Continent mean rel. anxious users over last 7 days",
       title = "Anxiety levels by continent over time")

ggsave("plots/Max/anx_cont_level.png", width = 12, height = 9)
```

# Just for Europe

-   Peak in Europe in on February 27 is probably bc of Ukraine war

```{r Anxiety continent level}
data_CTIS_policy %>% 
  group_by(continent, data.survey_date) %>% 
  filter(continent %in% "Europe") %>% 
  summarise(continent_mean_anxious7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  ggplot(aes(x = data.survey_date, y = continent_mean_anxious7d))+
  facet_wrap(vars(continent))+
    geom_line()
```

```{r}
data_CTIS_policy %>% filter(continent %in% "Europe") %>% filter(data.survey_date %in% as.Date("2022-02-27")) %>% dplyr::select(anxious_7d, data.country) %>% arrange(desc(anxious_7d))
```

```{r Anxiety continent level}
data_CTIS_policy %>% 
  group_by(continent, data.survey_date) %>% 
  summarise(continent_mean_depressd7d = mean(depressed_7d, na.rm = TRUE)) %>% 
  ggplot(aes(x = data.survey_date, y = continent_mean_depressd7d))+
  facet_wrap(vars(continent))+
    geom_line()
```



```{r}
# data: A dataframe containing the CTIS Survey data
# date: A string setting the date of the data (Format: 2022-03-22)
# variable: A string containing the variable name

plot_survey_world <- function(data, date, variable) {
  map_data <- data %>% filter(data.survey_date %in% as.Date(date))
  map_data <- joinCountryData2Map(map_data, nameJoinColumn = "data.iso_code")
  mapCountryData(map_data, nameColumnToPlot = variable, 
                 colourPalette = brewer.pal(7,"OrRd"), mapTitle = paste(variable, date))
}
head(data_CTIS_policy)
plot_survey_world(data_CTIS_policy, "2020-04-28", "anxious_7d")
plot_survey_world(data_CTIS_policy, "2021-01-28", "anxious_7d")
plot_survey_world(data_CTIS_policy, "2021-04-28", "anxious_7d")




```

# Creating map for shiny app
## Raw function for plotting for debugging
```{r}
# Continent: "Asia" "Africa" "Europe" "Americas" "Oceania" | valid inputs
area <- "Europe"
# Variable: "anxious_7d" "finance" "depressed_7d" 
variable_choices <- c("anxious_7d", "finance","depressed_7d", "food_security", "worried_become_ill")
variable <- "anxious_7d"
# survey_date
survey_date <- "2021-01-28"
# Mode: "plot" "view
mode <- "plot"


mapdata <- data_CTIS_map %>% 
  dplyr::select(continent, variable, geometry, data.country, data.iso_code, data.survey_date) %>% 
  filter(data.survey_date %in% as.Date(survey_date))

# If we filter for specific data we lose geometries for countries with no obs at that data therefore we need to add the geo data + NA for survey data
missing_countries <- setdiff(world$data.country, mapdata$data.country)
missing_data <- world %>% filter(data.country %in% missing_countries)
missing_columns <- setdiff(colnames(mapdata), colnames(missing_data))
for (column in 1:length(missing_columns)) {
  missing_data[missing_columns[column]] <- NA
}

mapdata <- rbind(mapdata, missing_data)

tmap_mode("plot")
tm_shape(mapdata)+
  tm_fill("anxious_7d")+
  tm_borders()
```

# plot_data_map
```{r}
plot_global_map <- function(data, variable, survey_date, mode = "plot", ...){
  mapdata <- data %>% 
  dplyr::select(continent, variable, geometry, data.country, data.iso_code, data.survey_date, school_closures, stay_home_requirements) %>% 
  filter(data.survey_date %in% as.Date(survey_date))
# If we filter for specific data we lose geometries for countries with no obs at that data therefore   we need to add the geo data + NA for survey data
  missing_countries <- setdiff(world$data.country, mapdata$data.country)
  missing_data <- world %>% filter(data.country %in% missing_countries)
  missing_columns <- setdiff(colnames(mapdata), colnames(missing_data))
  for (column in 1:length(missing_columns)) {
    missing_data[missing_columns[column]] <- NA
  }
  mapdata <- rbind(mapdata, missing_data)
  tmap_mode(mode)
  tm_shape(mapdata)+
    tm_fill(variable, ...)+
    tm_borders()
}
```


```{r}
plot_continent_map <- function(data, variable, area, survey_date, mode = "plot", ...){
  mapdata <- data %>% 
  dplyr::select(continent, variable, geometry, data.country, data.iso_code, data.survey_date, school_closures, stay_home_requirements) %>% 
  filter(data.survey_date %in% as.Date(survey_date)) %>% 
  filter(continent %in% area)

  tmap_mode(mode)
  tm_shape(mapdata)+
    tm_borders()+
    tm_polygons(variable, ...)
}
```

```{r}
plot_global_map(data_CTIS_map, variable = "anxious_7d",survey_date = "2021-01-28", mode = "view", id = "data.country", popup.vars = c("Anxiety last 7 days:" = "anxious_7d", "Stay home req. level" = "stay_home_requirements", "School closure level" = "school_closures"))

plot_global_map(data_CTIS_map, variable = "anxious_7d",survey_date = "2021-01-28", mode = "plot",  title='Anxiety last 7 days')+
   tm_layout(legend.position = c("left", "bottom"), 
            main.title= 'Global anxiety levels on 2021-01-28',
            legend.outside = TRUE,
            main.title.position = c('center', 'top'))


plot_continent_map(data_CTIS_map, variable = "anxious_7d", area = "Europe", survey_date = "2021-01-28", mode = "view", id = "data.country", popup.vars = c("Anxiety last 7 days:" = "anxious_7d", "Stay home req. level" = "stay_home_requirements", "School closure level" = "school_closures"))
```



=======
```{r Elisa}
########################### anxious_7d ########################################

# Sorted mean anxious_7d for each country
data_CTIS %>% group_by(data.country) %>% summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% arrange(desc(mean_anxious_7d)) 

# first 10 countries with highest anxious_7d
top10_mean_anxious_7d <- data_CTIS %>% group_by(data.country) %>% summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% arrange(desc(mean_anxious_7d)) %>% top_n(n = 10)

# vector of first 10 countries with highest anxious_7d
top10_countires_highest_mean_anxious <- top10_mean_anxious_7d$data.country

# Boxplot for first 10 countries with highest anxious_7d
data_CTIS %>% filter(data.country %in% top10_countires_highest_mean_anxious) %>% 
ggplot(aes(data.country, anxious_7d)) + 
  geom_boxplot() +
  theme(legend.position = "bottom") +
  geom_hline(yintercept = mean(data_CTIS$anxious_7d, na.rm = TRUE))


## Bar plot for first 10 countries with highest anxious_7d
data_CTIS %>% group_by(data.country) %>% filter(data.country %in% top10_countires_highest_mean_anxious) %>% summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% ggplot(aes(y = mean_anxious_7d, x = data.country)) + 
  geom_bar(stat = "identity") 

# Timeline plot for first 10 countries with highest anxious_7d
options(dplyr.summarise.inform = FALSE)
data_CTIS %>% group_by(data.country, data.survey_date) %>% filter(data.country %in% top10_countires_highest_mean_anxious) %>% summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>%
  ggplot(aes(x = data.survey_date, y = mean_anxious_7d)) + 
  facet_wrap(vars(data.country)) + geom_line()

## Timeline plot for Turkey + Tunisia
data_CTIS %>% group_by(data.country, data.survey_date) %>% filter(data.country %in% c("Tunisia", "Turkey")) %>% summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>%
  ggplot(aes(x = data.survey_date, y = mean_anxious_7d, color = data.country)) + geom_line()

```

```{r Elisa}
########################### depressed_7d #######################################

data_CTIS %>% group_by(data.country) %>% summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>% arrange(desc(mean_depressed_7d)) 

# first 10 countries with highest depressed_7d
top10_mean_depressed_7d <- data_CTIS %>% group_by(data.country) %>% summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>% arrange(desc(mean_depressed_7d)) %>% top_n(n = 10)

# vector of first 10 countries with highest depressed_7d
top10_countries_highest_mean_depressed <- top10_mean_depressed_7d$data.country

## countries in the top 10 for anxiety and depression -> all ten except for Tanzania
intersect(top10_countires_highest_mean_anxious, top10_countries_highest_mean_depressed)

# Boxplot for first 10 countries with highest depressed_7d
data_CTIS %>% filter(data.country %in% top10_countries_highest_mean_depressed) %>% 
ggplot(aes(data.country, depressed_7d)) + 
  geom_boxplot() +
  theme(legend.position = "bottom") +
  geom_hline(yintercept = mean(data_CTIS$depressed_7d, na.rm = TRUE))


## Bar plot for first 10 countries with highest depressed_7d
data_CTIS %>% group_by(data.country) %>% filter(data.country %in% top10_countries_highest_mean_depressed) %>% summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>% ggplot(aes(y = mean_depressed_7d, x = data.country)) + 
  geom_bar(stat = "identity") 

# Timeline plot for first 10 countries with highest depressed_7d
options(dplyr.summarise.inform = FALSE)
data_CTIS %>% group_by(data.country, data.survey_date) %>% filter(data.country %in% top10_countries_highest_mean_depressed) %>% summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>% ggplot(aes(x = data.survey_date, y = mean_depressed_7d)) + 
  facet_wrap(vars(data.country)) + geom_line()

## Timeline plot for Egypt + Iraq
data_CTIS %>% group_by(data.country, data.survey_date) %>% filter(data.country %in% c("Egypt", "Iraq")) %>% summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>%
  ggplot(aes(x = data.survey_date, y = mean_depressed_7d, color = data.country)) + geom_line()


```


```{r}
############## Code from Chris
# Number of NAs for each row
colSums(is.na(data_CTIS))

# Number of countries with equal or more than 600 observations
sum(table(data_CTIS$data.country) >= 600)

# Names of the countries with equal or more than 600 observations
names_sufficient_data <- names(which(table(data_CTIS$data.country) >= 600))
names_sufficient_data

# Get the 10 countries with the highest amount of observations
names(sort(table(data_CTIS$data.country), decreasing = TRUE)[1:10])

head(data_CTIS)

summary(data_CTIS)
############# End code from Chris
```


```{r Only analysis for countries with more than 600 obs.}

############## Countries with > 600 obs. #####################################


table(data_CTIS$data.country)

# Subset of countries >= 600 observations
countries_600_subset <- data_CTIS %>% filter(data.country %in% names_sufficient_data)
summary(countries_600_subset)

# top 10 countries with more than 600 observations
top10 <- names(sort(table(data_CTIS$data.country), decreasing = TRUE)[1:10])
top10_subset <- data_CTIS %>% filter(data.country %in% top10)

# Argentina dataset
Argentina_dataset <- top10_subset %>% 
  filter(data.country == "Argentina")

# Argentina plot
Argentina_dataset %>% ggplot(mapping = aes(x = anxious_7d, 
                         y = depressed_7d)) + geom_jitter() + geom_abline()


# Arrange by mean_anxious_7d for countries with > 600 obs. 
mean_anxious_7d_arranged <- data_CTIS %>% 
  filter(data.country %in% names_sufficient_data) %>%
  group_by(data.country) %>% 
  summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  arrange(mean_anxious_7d)
mean_anxious_7d_arranged

## first 10 countries with highest anxious and > 600 obs.
mean_anxious_7d_top10 <- top10_subset %>%
  group_by(data.country) %>% 
  summarise(mean_anxious_7d = mean(anxious_7d, na.rm = TRUE)) %>% 
  arrange(mean_anxious_7d)
mean_anxious_7d_top10 

# Bar plots for first 10 countries with > 600 obs. 
# anxious_7d for top10 countries
ggplot(data = mean_anxious_7d_top10, 
       aes(y = mean_anxious_7d, x = data.country)) + 
  geom_bar(stat = "identity") 

## Arrange by mean_depressed_7d for top 10 countries (> 600 obs.)
mean_depressed_7d_top10 <- top10_subset %>%
  group_by(data.country) %>% 
  summarise(mean_depressed_7d = mean(depressed_7d, na.rm = TRUE)) %>% 
  arrange(mean_depressed_7d)
mean_depressed_7d_top10 

# depressed_7d for top10 countries
ggplot(data = mean_depressed_7d_top10, 
       aes(y = mean_depressed_7d, x = data.country)) + 
  geom_bar(stat = "identity")

## Arrange by food_security, finance, worried_become_ill for top 10 countries
mean_food_finance_ill_top10 <- top10_subset %>%
  group_by(data.country) %>% 
      summarise(mean_ill_7d = mean(worried_become_ill, na.rm = TRUE),
                mean_depressed_7d = mean(depressed_7d, na.rm = TRUE), 
                mean_anxious_7d = mean(anxious_7d, na.rm = TRUE), 
                mean_food_7d = mean(food_security, na.rm = TRUE), 
                mean_finance_7d = mean(finance, na.rm = TRUE))

mean_food_finance_ill_top10


# food security for top10 countries
ggplot(data = mean_food_finance_ill_top10, 
       aes(y = mean_food_7d, x = data.country)) + 
  geom_bar(stat = "identity")

# finance for top10 countries
ggplot(data = mean_food_finance_ill_top10, 
       aes(y = mean_finance_7d, x = data.country)) + 
  geom_bar(stat = "identity")

# worried_become_ill for top10 countries
ggplot(data = mean_food_finance_ill_top10, 
       aes(y = mean_ill_7d, x = data.country)) + 
  geom_bar(stat = "identity")



```

```{r}
# Create a new data.frame with the newly formatted date field
data_CTIS_new <- data_CTIS %>%
  mutate(DATE = as.Date(data.survey_date, format = "%m/%d/%y"))


data_CTIS_new <- data_CTIS_new %>% drop_na() %>% mutate(month_new = as.Date(data.survey_date, 
                          format = "%m/%y")) %>% group_by(month_new)

```

## Correlations
```{r}
## daily correlation
data_CTIS %>%  drop_na() %>% group_by(data.survey_date) %>%
  summarize(COR=cor(depressed_7d,anxious_7d)) %>% ggplot(aes(x = data.survey_date, y = COR)) + geom_line() + 
  xlab("Date")

cor_daily <- data_CTIS %>%  drop_na() %>% group_by(data.survey_date) %>%
  summarize(COR=cor(depressed_7d,anxious_7d))

highchart() %>% 
  hc_add_series(cor_daily$COR, type = "line")
######################################## CORRECT
# month
data_CTIS_policy_large$month <- floor_date(data_CTIS_policy_large$data.survey_date, "1 month")
month_cor <- data_CTIS_policy_large %>% group_by(month) %>% summarise(correlation = cor(anxious_7d, depressed_7d))
# week
data_CTIS_policy_large$week <- floor_date(data_CTIS_policy_large$data.survey_date, "1 week")
week_cor <- data_CTIS_policy_large %>% group_by(week) %>% summarise(correlation = cor(anxious_7d, depressed_7d))
######################################### END
highchart() %>% 
  hc_add_series(month_cor$correlation, type = "line")
```


```{r}
# correlation table 
res <- data_CTIS %>% select(anxious_7d, depressed_7d, worried_become_ill, 
                            finance, food_security) %>% drop_na()
correlation <- cor(res)
correlation

## add week numbers to the dataset
data_CTIS$week_num <- floor_date(data_CTIS$data.survey_date, "1 week")

data_CTIS$week_num <- strftime(data_CTIS$data.survey_date, format = "%V")
## add month to the dataset
data_CTIS$month <- floor_date(data_CTIS$data.survey_date, "1 month")

# Ukraine has the highest correlation between depressed and anxious
data_CTIS %>%
  group_by(data.country) %>%
  summarize(COR=cor(depressed_7d,anxious_7d)) %>% drop_na() %>% arrange(desc(COR))
Ukraine_data <- data_CTIS %>%  filter(data.country == "Ukraine")

## correlations
# remove not needed variables
ukraine <- Ukraine_data %>% select(-c(1:4), -10, -11)
summary(ukraine)


corr_cross(ukraine, rm.na = T, max_pvalue = 0.05, top = 15, grid = T)
```

```{r Turkey Correlation anxious and depressed}
## with dygraphs
# Turkey has highest anxious_7d
data_CTIS_Turkey <- data_CTIS %>% filter(data.country == "Turkey")
don <- xts(x = data_CTIS_Turkey$anxious_7d, order.by = data_CTIS_Turkey$data.survey_date)

p <- dygraph(don) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
p
```


```{r Daily correlation anxious and depressed}
## corr daily
daily <- xts(x = cor_daily$COR, order.by = cor_daily$data.survey_date)
d <- dygraph(daily) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
d
```


```{r}
## monthly correlation
monthly_correlation <- xts(x = month_cor$correlation, order.by = month_cor$month)
m <- dygraph(monthly_correlation) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
m
```


```{r}
## weekly correlation, continent level
data_CTIS_policy$week <- floor_date(data_CTIS_policy$data.survey_date, "1 week")
week_cor <- data_CTIS_policy %>% drop_na() %>% group_by(week, continent) %>% summarise(correlation = cor(anxious_7d, depressed_7d))
# plot
ggplot(week_cor, aes(x = week, y = correlation))+ facet_wrap(~ continent, ncol = 2) + geom_line() +
  theme(axis.text = element_text(size = 8))
ggsave("plots/Elisa/weekly_correlation_anxious_dep_continent.png", width = 12, height = 9)
```


```{r}
## weekly
week_correlation <- xts(x = week_cor$correlation, order.by = week_cor$week)
w <- dygraph(week_correlation) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
w
```


```{r}
## weekly in Ukraine
week_cor_Ukraine <- data_CTIS_policy %>% filter(data.country == "Ukraine") %>% 
  group_by(week) %>% summarise(correlation = cor(anxious_7d, depressed_7d))


data_CTIS_policy$week <- floor_date(data_CTIS_policy$data.survey_date, "1 week")

week_cor <- data_CTIS_policy %>% filter(data.country == "Ukraine") %>%
  group_by(week) %>% summarise(correlation = cor(anxious_7d, depressed_7d))
# plot ukraine
week_cor_ukraine <- xts(x = week_cor$correlation, order.by = week_cor$week)
w_ukraine <- dygraph(week_cor_ukraine) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
w_ukraine

```

```{r Correlation}
## correlation with more plot types
res <- data_CTIS %>% select(anxious_7d, depressed_7d, worried_become_ill, 
                            finance, food_security) %>% drop_na()
ggpairs(res) 
ggsave("plots/Elisa/correlation_policy_data.png", width = 12, height = 9)

## correlation of each variable as a square
data_CTIS_map_ohnegeom <- st_set_geometry(data_CTIS_map, NULL)
cor_map <- data_CTIS_map_ohnegeom %>% select(anxious_7d, depressed_7d, finance, 
                                    food_security, worried_become_ill,
                                    retail_and_recreation, grocery_and_pharmacy,
                                    residential, transit_stations, parks, 
                                    workplaces) %>% drop_na()

ggcorr(cor_map, method = c("everything", "pearson")) 


```

```{r Correlation}
## Correlation in the data set with policy restrictions
ctis_map_correlation <- data_CTIS_map_ohnegeom %>% select(-c(1:4), -19) %>% drop_na()
summary(ctis_map_correlation)
# transform variables as numeric
ctis_map_correlation[1:14] <- sapply(ctis_map_correlation, as.numeric)
cor_map2 <- cor(ctis_map_correlation)

# 
ggcorr(cor_map2, method = c("everything", "pearson")) 

```


```{r}
## countries with high correlation (> 0.9) between food security and finance

cor_food_finance <- data_CTIS_map_ohnegeom %>%
  group_by(data.country) %>% drop_na() %>%
  summarize(COR = cor(food_security,finance)) %>% arrange(desc(COR)) %>% top_n(15)

top10 <- cor_food_finance $data.country
top10
data <- head(cor_food_finance, 15)

data_CTIS_map_ohnegeom %>% filter(data.country %in% top10) %>% ggplot(aes(x=finance, y=food_security)) +
  geom_point()


# "Brazil"  "Mexico" "Taiwan""Indonesia", "Thailand", "Colombia", "Argentina", 
# "Portugal"   
# "Philippines", "Chile", "India", "Peru", "France", "Spain", "Italy" between finance and food security

```


```{r}
## monthly correlation
# month
data_CTIS_policy_large_wo_na <- data_CTIS_policy_large %>% drop_na()
data_CTIS_policy_large_wo_na$month <- floor_date(data_CTIS_policy_large_wo_na$data.survey_date, "1 month") 
month_cor_food_finance <- data_CTIS_policy_large_wo_na %>% group_by(month) %>% summarise(correlation = cor(finance, food_security))

## plot
monthly_cor_food_finance <- xts(x = month_cor_food_finance$correlation, order.by = month_cor_food_finance$month)
m_food_finance <- dygraph(monthly_cor_food_finance) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)
m_food_finance

```



```{r}
## work with this dataset: data_CTIS_map_ohnegeom
head(data_CTIS_map_ohnegeom)

# table
data_CTIS_map_ohnegeom %>% group_by(data.country) %>% drop_na() %>% summarise(cor = cor(finance, 
                          food_security)) %>% arrange(desc(cor))

# bar plot for first 10 countries with high correlation between finance and food security
cor_data <- data_CTIS_map_ohnegeom %>% group_by(data.country) %>% drop_na() %>% summarise(cor = cor(finance, 
                          food_security)) %>% arrange(desc(cor)) %>% top_n(10)

ggplot(cor_data) + geom_col(aes(x = data.country, y = cor))


```

```{r Correlation Plots with global map}
## Correlation Plot
# new dataset with correlation between anxious and depressed
data_cor_plot <- data_CTIS_map %>% drop_na() %>% group_by(data.country, data.iso_code) %>% summarise(cor_anx_dep = cor(anxious_7d, depressed_7d))
data_cor_plot <- data_cor_plot %>% drop_na()
# join data
joinData <- joinCountryData2Map(data_cor_plot,
            nameJoinColumn = "data.iso_code")
# map plot
mapCountryData(joinData, nameColumnToPlot = "cor_anx_dep", catMethod = "FixedWidth", colourPalette = brewer.pal(7,"OrRd"))
```

```{r Correlation Plots with global map}
## Correlation Plot

# new dataset with correlation between finance and food security
data_cor_fin_food <- data_CTIS_map %>% drop_na() %>% group_by(data.country, data.iso_code) %>% summarise(cor_fin_food = cor(finance, food_security))
data_cor_fin_food <- data_cor_fin_food %>% drop_na()
# join data
joinData2 <- joinCountryData2Map(data_cor_fin_food,
            nameJoinColumn = "data.iso_code")
# map plot
mapCountryData(joinData2, nameColumnToPlot = "cor_fin_food", catMethod = "FixedWidth", colourPalette = brewer.pal(7,"OrRd"))
```


## Code from Max -> using it for correlation interactive plots (Elisa) ## 

# plot_data_map
```{r}

plot_global_map <- function(data, variable, mode = "plot", ...){
  mapdata <- data_cor_fin_food %>% 
  dplyr::select(variable, geometry, data.country, data.iso_code) 
# If we filter for specific data we lose geometries for countries with no obs at that data therefore   we need to add the geo data + NA for survey data
  missing_countries <- setdiff(world$data.country, mapdata$data.country)
  missing_data <- world %>% filter(data.country %in% missing_countries)
  missing_columns <- setdiff(colnames(mapdata), colnames(missing_data))
  for (column in 1:length(missing_columns)) {
    missing_data[missing_columns[column]] <- NA
  }
  mapdata <- rbind(mapdata, missing_data)
  tmap_mode(mode)
  tm_shape(mapdata)+
    tm_fill(variable, ...)+
    tm_borders()
}
```


```{r}
plot_continent_map <- function(data, variable, area, survey_date, mode = "plot", ...){
  mapdata <- data_cor_fin_food %>% 
  dplyr::select(variable, geometry, data.country, data.iso_code) 

  tmap_mode(mode)
  tm_shape(mapdata)+
    tm_borders()+
    tm_polygons(variable, ...)
}
```

```{r}
plot_global_map(data_cor_fin_food, variable = "cor_fin_food", mode = "view", id = "data.country", popup.vars = c("Correlation between finance and food security:" = "cor_fin_food"))

plot_global_map(data_cor_fin_food, variable = "cor_fin_food", mode = "plot", 
                title ='Finance and food security correlation')+
   tm_layout(legend.position = c("left", "bottom"), 
            main.title= 'Correlation between finance and food security',
            legend.outside = TRUE,
            main.title.position = c('center', 'top'))


plot_continent_map(data_cor_fin_food, variable = "cor_fin_food",
                   area = "Europe", mode = "view", id = "data.country", 
                   popup.vars = c("Correlation between finance and food security:" = "cor_fin_food"))
```