--- title: "Covid-19 lockdown effects on LBV hedgehog observation numbers" author: "Fabio" date: "13/11/2021" output: pdf_document: default html_document: default editor_options: chunk_output_type: console --- ```{r, include = FALSE, echo=FALSE} options(RStudioGD.backend = "cairo") knitr::opts_chunk$set(echo = FALSE) ``` # Setup Firstly, the necessary packages need to be loaded. ```{r load packages, include=FALSE} library(ggplot2) library(tidyverse) library(lubridate) library(mgcv) library(visreg) library(itsadug) library(patchwork) ``` Then the data will be loaded in. ```{r load data, include=FALSE} DATA <- read.csv(file = "ENTER_DATA_LOCATION") DATA$igel_sicht <- as.Date(DATA$igel_sicht, format="%d/%m/%Y") ``` The dates will be summarized per day per month per year ```{r summarize per day per month per year, include=FALSE} DATA_MONTH <- DATA %>% dplyr::mutate(Year = lubridate::year(igel_sicht), Month = lubridate::month(igel_sicht), Day = lubridate::day(igel_sicht), Week = lubridate::week(igel_sicht)) ``` The years 2016-2019 will be grouped together, and the year 2020 will be its own group in order to compare this time period for the year 2020 with its expected values based on previous years. ```{r selection on year and weeks, include=FALSE} DATA_MONTH2 <- DATA_MONTH %>% mutate(Year_Groups = case_when( between(Year, 2016, 2019) ~ "2016-2019", between(Year, 2020, 2020) ~ "2020", TRUE ~ NA_character_ )) %>% filter(Year != 2015 & Year != 2021) DATA_CALC_WEEK <- DATA_MONTH2 %>% group_by(Year, Week, Year_Groups) %>% summarise(count = n()) DATA_CALC_WEEK$Year <- as.factor(DATA_CALC_WEEK$Year) DATA_CALC_WEEK$Year_Groups <- as.factor(DATA_CALC_WEEK$Year_Groups) ``` # Are there more sightings of hedgehogs in 2020 than normal? ```{r GAM setup and results, include=FALSE} set.seed(1) GAM <- gam(data = DATA_CALC_WEEK , count ~ s(Week, by = Year_Groups, k=53) , method = "GCV.Cp") summary(GAM) visreg(GAM, "Week",by = "Year_Groups", overlay = FALSE, band = TRUE) visreg(GAM, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE) visreg(GAM, "Week",by = "Year_Groups", overlay = TRUE, band = FALSE) par(mfrow = c(2,1)) plot(GAM) par(mfrow = c(1,1)) par(mfrow = c(2,2)) gam.check(GAM) par(mfrow = c(1,1)) plot_diff(GAM, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE) ``` ```{r Visualization preparation, echo=FALSE, warning=FALSE, results='hide'} P1 <- visreg(GAM, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE, gg=TRUE)+ geom_segment(aes(x = 13, xend = 19, y = -200, yend = -200), color = "#4BC34B", size = 5)+ theme_classic()+ theme(legend.title = element_blank(), legend.position = "top")+ ylab("n hedgehog \n observations") P2 <- plot_diff(GAM, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE, plot = FALSE) P2$diff_to_0 <- P2$est-P2$sim.CI ``` ```{r Visializations made in ggplot and ready for one graph, echo=FALSE, warning=FALSE, include=FALSE} P2 <- P2 %>% mutate(Coloration = case_when( between(diff_to_0, -2000, 0) ~ "C_1", between(diff_to_0, 0, 2000) ~ "C_2" )) P2_P <- ggplot(data = P2, aes(x = Week, y = est))+ geom_segment(aes(x = 13, xend = 19, y = -550, yend = -550), color = "#4BC34B", size = 5)+ geom_line(aes(colour=(Coloration == "C_1"), group=1))+ geom_ribbon(aes(ymin = est - sim.CI, ymax = est + sim.CI), alpha = 0.1)+ geom_hline(yintercept = 0, linetype = "dashed")+ theme_classic()+ theme(legend.position = "none")+ geom_vline(xintercept = 15.110553 , linetype = "dotted", color = "red")+ geom_vline(xintercept = 17.723618, linetype = "dotted", color = "red")+ scale_color_manual(values=c("red", "black"))+ ylab("Est. diff. \n [2020] - [2016-2019]") P1 / P2_P + plot_annotation(tag_levels = 'A') summary(GAM) ``` # Is this difference because there are more hedgehogs, just because there are more people? ```{r data preparation people observing, include=FALSE} DATA_WEEK3 <- DATA_MONTH %>% filter(Year != 2015 & Year != 2021) %>% filter(Email_ID != 0) %>% mutate(Year_Groups = case_when( between(Year, 2016, 2019) ~ "2016-2019", between(Year, 2020, 2020) ~ "2020" )) DATA_CALC_WEEK_3_1 <- DATA_WEEK3 %>% group_by(Year, Week, Year_Groups) %>% summarise(NUM = n_distinct(Email_ID)) DATA_CALC_WEEK_3_1$Year_Groups <- as.factor(DATA_CALC_WEEK_3_1$Year_Groups) DATA_CALC_WEEK_3_2 <- DATA_WEEK3 %>% group_by(Year, Week, Year_Groups) %>% summarise(NUM = n()/ n_distinct(Email_ID)) DATA_CALC_WEEK_3_2$Year_Groups <- as.factor(DATA_CALC_WEEK_3_2$Year_Groups) ``` ## The increase in the amount of people that did hedgehog observations within the first Covid-19 lockdown in Germany ```{r, analyses people observing, include=FALSE} GAM_2<- gam(data = DATA_CALC_WEEK_3_1 , NUM ~ s(Week, by = Year_Groups, k = 53), method = "GCV.Cp") summary(GAM_2) visreg(GAM_2, "Week",by = "Year_Groups", overlay = FALSE, band = TRUE) visreg(GAM_2, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE) visreg(GAM_2, "Week",by = "Year_Groups", overlay = TRUE, band = FALSE) par(mfrow = c(2,1)) plot(GAM_2) par(mfrow = c(1,1)) par(mfrow = c(2,2)) gam.check(GAM_2) par(mfrow = c(1,1)) plot_diff(GAM_2, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE) ``` ```{r Visualization preparation pt 2, echo=FALSE, warning=FALSE, results='hide'} P3 <- visreg(GAM_2, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE, gg=TRUE)+ geom_segment(aes(x = 13, xend = 19, y = -150, yend = -150), color = "#4BC34B", size = 5)+ theme_classic()+ theme(legend.title = element_blank(), legend.position = "top")+ ylab("n registered \n observers") P4 <- plot_diff(GAM_2, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE, plot = FALSE) P4$diff_to_0 <- P4$est-P4$sim.CI ``` ```{r Visializations made in ggplot and ready for one graph pt 2, echo=FALSE, warning=FALSE, include=FALSE} P4 <- P4 %>% mutate(Coloration = case_when( between(diff_to_0, -2000, 0) ~ "C_1", between(diff_to_0, 0, 2000) ~ "C_2" )) P4_P <- ggplot(data = P4, aes(x = Week, y = est))+ geom_segment(aes(x = 13, xend = 19, y = -500, yend = -500), color = "#4BC34B", size = 5)+ geom_line(aes(colour=(Coloration == "C_1"), group=1))+ geom_ribbon(aes(ymin = est - sim.CI, ymax = est + sim.CI), alpha = 0.1)+ geom_hline(yintercept = 0, linetype = "dashed")+ theme_classic()+ theme(legend.position = "none")+ geom_vline(xintercept = 15.110553 , linetype = "dotted", color = "red")+ geom_vline(xintercept = 17.201005, linetype = "dotted", color = "red")+ scale_color_manual(values=c("red", "black"))+ ylab("Est. diff. \n [2020] - [2016-2019]") P3 / P4_P + plot_annotation(tag_levels = 'A') summary(GAM_2) ``` ## The amount of hedgehog sightings per person, however, stayed the same. ```{r analyses number of hedgehog sightings per person, include=FALSE} GAM_3 <- gam(data = DATA_CALC_WEEK_3_2 , NUM ~ s(Week, by = Year_Groups, k=53) , method = "GCV.Cp") summary(GAM_3) visreg(GAM_3, "Week",by = "Year_Groups", overlay = FALSE, band = TRUE) visreg(GAM_3, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE) visreg(GAM_3, "Week",by = "Year_Groups", overlay = TRUE, band = FALSE) par(mfrow = c(2,1)) plot(GAM_3) par(mfrow = c(1,1)) par(mfrow = c(2,2)) gam.check(GAM_3) par(mfrow = c(1,1)) plot_diff(GAM_3, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE) ``` ```{r Visualization preparation pt 3, echo=FALSE, warning=FALSE, results='hide'} P5 <- visreg(GAM_3, "Week",by = "Year_Groups", overlay = TRUE, band = TRUE, gg=TRUE)+ geom_segment(aes(x = 13, xend = 19, y = 0.9, yend = 0.9), color = "#4BC34B", size = 5)+ theme_classic()+ theme(legend.title = element_blank(), legend.position = "top")+ ylab("Observations per \n registered observer") P6 <- plot_diff(GAM_3, view = "Week", comp=list(Year_Groups=c("2020", "2016-2019")), rm.ranef= FALSE, sim.ci = TRUE, plot = FALSE) P6$diff_to_0 <- P6$est+P6$sim.CI ``` ```{r Visializations made in ggplot and ready for one graph pt 3, echo=FALSE, warning=FALSE, include=FALSE} P6 <- P6 %>% mutate(Coloration = case_when( between(diff_to_0, -2000, 0) ~ "C_1", between(diff_to_0, 0, 2000) ~ "C_2" )) P6_P <- ggplot(data = P6, aes(x = Week, y = est))+ geom_segment(aes(x = 13, xend = 19, y = -0.2, yend = -0.2), color = "#4BC34B", size = 5)+ geom_line(aes(colour=(Coloration == "C_1"), group=1))+ geom_ribbon(aes(ymin = est - sim.CI, ymax = est + sim.CI), alpha = 0.1)+ geom_hline(yintercept = 0, linetype = "dashed")+ theme_classic()+ theme(legend.position = "none")+ geom_vline(xintercept = 22.688442 , linetype = "dotted", color = "red")+ geom_vline(xintercept = 24.256281, linetype = "dotted", color = "red")+ scale_color_manual(values=c("black", "red"))+ ylab("Est. diff. \n [2020] - [2016-2019]") P5 / P6_P + plot_annotation(tag_levels = 'A') summary(GAM_3) ``` ```{r Final combined plot , echo=FALSE, warning=FALSE} PLOT <- (((P1 | P3 | P5) + plot_layout(guides = "collect") & theme(legend.position='none')) / ((P2_P | P4_P | P6_P ) + plot_layout(guides = "collect") & theme(legend.position='none'))) + plot_annotation(tag_levels = c('1')) PLOT[[1]] <- PLOT[[1]] + plot_layout(tag_level = 'new') PLOT[[2]] <- PLOT[[2]] + plot_layout(tag_level = 'new') PLOT + plot_annotation(tag_levels = c('1', 'A'), tag_sep = '.', tag_suffix = ':') & theme(plot.tag.position = c(0, 1), plot.tag = element_text(size = 12, hjust = 0, vjust = 0)) summary(GAM) summary(GAM_2) summary(GAM_3) ``` #Differences in number of observations during lockdown based on imperviousness levels. ```{r} DATA_CALC_WEEK_IMP <- DATA_MONTH2 %>% filter(Week > 12 & Week < 20) %>% mutate(Imperviousness_proportion = case_when( between(Mean_imperviousness, 0, 19.99999999999) ~ "0-20", between(Mean_imperviousness, 20, 39.9999999999) ~ "20-40", between(Mean_imperviousness, 40, 59.9999999999) ~ "40-60", between(Mean_imperviousness, 60, 79.9999999999) ~ "60-80", between(Mean_imperviousness, 80, 100) ~ "80-100", TRUE ~ NA_character_ )) %>% group_by(Year, Year_Groups, Imperviousness_proportion) %>% summarise(count = n()) %>% group_by(Year) %>% mutate(percent = prop.table(count)*100) DATA_CALC_WEEK_IMP$Year_Groups <- as.factor(DATA_CALC_WEEK_IMP$Year_Groups) P7 <- ggplot(DATA_CALC_WEEK_IMP,aes(x=Year,y=percent, fill = Imperviousness_proportion))+ geom_bar(position="stack", stat="identity")+ theme_classic()+ ggsci::scale_fill_npg()+ ylab("Percent observations \n per impervious \n surface level") + scale_y_continuous(expand=c(0,0))+ guides(fill=guide_legend(title="Percentage \n impervious surface")) DATA_WEEK3_IMP <- DATA_MONTH %>% filter(Year != 2015 & Year != 2021) %>% filter(Week > 12 & Week < 20) %>% mutate(Imperviousness_proportion = case_when( between(Mean_imperviousness, 0, 19.99999999999) ~ "0-20", between(Mean_imperviousness, 20, 39.9999999999) ~ "20-40", between(Mean_imperviousness, 40, 59.9999999999) ~ "40-60", between(Mean_imperviousness, 60, 79.9999999999) ~ "60-80", between(Mean_imperviousness, 80, 100) ~ "80-100", TRUE ~ NA_character_ )) %>% filter(Email_ID != 0) %>% mutate(Year_Groups = case_when( between(Year, 2016, 2019) ~ "2016-2019", between(Year, 2020, 2020) ~ "2020" )) DATA_CALC_WEEK_3_1_IMP <- DATA_WEEK3_IMP %>% group_by(Year, Year_Groups, Imperviousness_proportion) %>% summarise(NUM = n_distinct(Email_ID)) %>% group_by(Year) %>% mutate(percent = prop.table(NUM)*100) P8 <- ggplot(DATA_CALC_WEEK_3_1_IMP,aes(x=Year,y=percent, fill = Imperviousness_proportion))+ geom_bar(position="stack", stat="identity")+ theme_classic()+ ggsci::scale_fill_npg()+ ylab("Percent registered \n users per impervious \n surface level")+ scale_y_continuous(expand=c(0,0))+ guides(fill=guide_legend(title="Percentage \n impervious surface")) DATA_CALC_WEEK_3_2_IMP <- DATA_WEEK3_IMP %>% group_by(Year, Year_Groups, Imperviousness_proportion) %>% summarise(NUM = n()/ n_distinct(Email_ID)) %>% group_by(Year) %>% mutate(percent = prop.table(NUM)*100) P9 <- ggplot(DATA_CALC_WEEK_3_2_IMP,aes(x=Year,y=percent, fill = Imperviousness_proportion))+ geom_bar(position="dodge", stat="identity")+ theme_classic()+ ggsci::scale_fill_npg()+ ylab("n observations \n per registered user")+ scale_y_continuous(expand=c(0,0))+ guides(fill=guide_legend(title="Percentage \n impervious surface")) P7 + P8 # Try 2020 proportion prediction method similar to Crimmins et al 2021 DATA_CALC_WEEK_IMP_DIFF <- DATA_CALC_WEEK_IMP %>% filter( Year != 2020) TEST <- lm(data = DATA_CALC_WEEK_IMP_DIFF, percent ~ Year * Imperviousness_proportion ) summary(TEST) par(mfrow = c(2,2)) plot(TEST) par(mfrow = c(1,1)) TEST_predict <- predict(TEST, newdata = DATA_CALC_WEEK_IMP , interval = "prediction") TEST_predict <- data.frame(TEST_predict) DATA_CALC_WEEK_IMP_TEST <- cbind(DATA_CALC_WEEK_IMP, TEST_predict) ggplot(data = DATA_CALC_WEEK_IMP_TEST,aes( x = Year, color = Imperviousness_proportion, fill = Imperviousness_proportion)) + geom_ribbon(aes(ymin=lwr, ymax=upr), alpha = 0.3)+ geom_line(aes(y= fit))+ geom_point(aes(y = percent))+ facet_wrap(~Imperviousness_proportion)+ theme_classic()+ ylab("Percent observations in imperviousness class")+ theme(legend.position = "none") DATA_CALC_WEEK_3_1_IMP_DIFF <- DATA_CALC_WEEK_3_1_IMP %>% filter( Year != 2020) TEST_2 <- lm(data = DATA_CALC_WEEK_3_1_IMP_DIFF, percent ~ Year * Imperviousness_proportion ) summary(TEST_2) par(mfrow = c(2,2)) plot(TEST_2) par(mfrow = c(1,1)) TEST_2_predict <- predict(TEST_2, newdata = DATA_CALC_WEEK_3_1_IMP , interval = "prediction") TEST_2_predict <- data.frame(TEST_2_predict) DATA_CALC_WEEK_IMP_TEST_2 <- cbind(DATA_CALC_WEEK_3_1_IMP, TEST_2_predict) ggplot(data = DATA_CALC_WEEK_IMP_TEST_2,aes( x = Year, color = Imperviousness_proportion, fill = Imperviousness_proportion)) + geom_ribbon(aes(ymin=lwr, ymax=upr), alpha = 0.3)+ geom_line(aes(y= fit))+ geom_point(aes(y = percent))+ facet_wrap(~Imperviousness_proportion)+ theme_classic()+ ylab("Percent observers in imperviousness class")+ theme(legend.position = "none") DATA_CALC_WEEK_IMP_TEST_X <- DATA_CALC_WEEK_IMP_TEST %>% filter( Year == 2020) DATA_CALC_WEEK_IMP_TEST_2_X <- DATA_CALC_WEEK_IMP_TEST_2 %>% filter( Year == 2020) P10 <- ggplot() + geom_col(data = DATA_CALC_WEEK_IMP_TEST_X, aes(x = Imperviousness_proportion, y = fit, fill = "#0C9BBD"), position = position_nudge(x = -0.12) , color = "black", width = 0.2) + geom_errorbar(data = DATA_CALC_WEEK_IMP_TEST_X, aes(x = Imperviousness_proportion, ymax = upr, ymin = lwr), width = 0.15, position = position_nudge(x = -0.12)) + geom_point(data = DATA_CALC_WEEK_IMP_TEST_X, aes(x = Imperviousness_proportion, y = percent), position = position_nudge(x = -0.12), shape = 8)+ geom_col(data = DATA_CALC_WEEK_IMP_TEST_2_X, aes(x = Imperviousness_proportion, y = fit , fill = "#FFB700"), position = position_nudge(x = 0.12), color = "black", width = 0.2) + geom_errorbar(data = DATA_CALC_WEEK_IMP_TEST_2_X, aes(x = Imperviousness_proportion, ymax = upr, ymin = lwr), width = 0.15, position = position_nudge(x = 0.12)) + geom_point(data = DATA_CALC_WEEK_IMP_TEST_2_X, aes(x = Imperviousness_proportion, y = percent), position = position_nudge(x = 0.12), shape = 8)+ ylab("Percentage")+ xlab("Imperviousness percentage class")+ coord_cartesian(ylim = c(0, 40), expand = FALSE)+ theme_classic()+ scale_fill_identity(name = '', guide = 'legend',labels = c('Observations', 'Observers')) P10 # There is no deviation from the prediction that you would expect for the percentage of observations anywhere in 2020. There might however be a small slightly lower percentage of people doing observations on the lowest level of imperviousness in 2020 compared to what you would expect. ```