---
title: "01_Analysis_SpatialScale"
author: "Camille"
date: "3/1/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Generates code for spatial scale robustness check

```{r Load data and package}
#load packages
library(tidyverse)
library(lfe)

#load data 
master_full <- read.csv("DataKeep/MasterSheet_full_points.csv")
master_full$Tourism <- master_full$Tourism/1000
master_full$PercentDeforest <- master_full$PercentDeforest*100
master_full <- master_full %>% filter(Year >2000)
master_full$OBJECTID2 <- master_full$OBJECTID

master_full_buffer <- read.csv("DataKeep/MasterSheet_full_points_buffer.csv")
master_full_buffer$Tourism <- master_full_buffer$Tourism/1000
master_full_buffer$PercentDeforest <- master_full_buffer$PercentDeforest*100
master_full_buffer <- master_full_buffer %>% filter(Year >2000)
master_full_buffer$OBJECTID2 <- paste0(master_full_buffer$OBJECTID, "New")

master_full_combined <- rbind(master_full_buffer, master_full)
master_full_combined$OBJECTID <- master_full_combined$OBJECTID2


mean(na.omit(master_full$PercentDeforest))
mean(na.omit(master_full_buffer$PercentDeforest))

```

```{r Full Analysis}
#FULL ANALYSIS
fem_full_formula <- as.formula("PercentDeforest ~ Tourism + PopDensity + precipitation + OilPriceAdj | Year + OBJECTID | 0 | Name")
dd_reg_full <- felm(fem_full_formula, data = master_full)
summary_dd_reg_full <- summary(dd_reg_full)

main_models <- as.data.frame(summary_dd_reg_full$coefficients)
main_models$Variable <- rownames(main_models)

main_models$`Pr(>|t|)` <- round(main_models$`Pr(>|t|)`, digits = 3)
main_models$Estimate <- round(main_models$Estimate, digits = 3)
main_models$`Cluster s.e.` <- round(main_models$`Cluster s.e.`, digits = 3)
main_models$`t value` <- round(main_models$`t value`, digits = 3)
main_models$Model <- "90m Spatial Resolution"
main_models <- main_models %>% select(Estimate:`Pr(>|t|)`, Model, Variable)

100*main_models[1,1]/mean(na.omit(master_full$PercentDeforest))

write.csv(main_models, "Outputs/points_summary_output.csv")


#BUFFER ANALYSIS
fem_full_formula_buffer <- as.formula("PercentDeforest ~ Tourism + PopDensity + precipitation + OilPriceAdj | Year + OBJECTID | 0 | Name")
dd_reg_full_buffer <- felm(fem_full_formula_buffer, data = master_full_buffer)
summary_dd_reg_full_buffer <- summary(dd_reg_full_buffer)

main_models_buffer <- as.data.frame(summary_dd_reg_full_buffer$coefficients)
main_models_buffer$Variable <- rownames(main_models_buffer)

main_models_buffer$`Pr(>|t|)` <- round(main_models_buffer$`Pr(>|t|)`, digits = 3)
main_models_buffer$Estimate <- round(main_models_buffer$Estimate, digits = 3)
main_models_buffer$`Cluster s.e.` <- round(main_models_buffer$`Cluster s.e.`, digits = 3)
main_models_buffer$`t value` <- round(main_models_buffer$`t value`, digits = 3)
main_models_buffer$Model <- "90m Spatial Resolution"
main_models_buffer <- main_models_buffer %>% select(Estimate:`Pr(>|t|)`, Model, Variable)

100*main_models_buffer[1,1]/mean(na.omit(master_commune$PercentDeforest))

write.csv(main_models_buffer, "Outputs/points_summary_output_buffer.csv")


#COMBINED ANALYSIS
fem_full_formula_combined <- as.formula("PercentDeforest ~ Tourism + PopDensity + precipitation + OilPriceAdj | Year + OBJECTID | 0 | Name")
dd_reg_full_combined <- felm(fem_full_formula_combined, data = master_full_combined)
summary_dd_reg_full_combined <- summary(dd_reg_full_combined)

main_models_combined <- as.data.frame(summary_dd_reg_full_combined$coefficients)
main_models_combined$Variable <- rownames(main_models_combined)

main_models_combined$`Pr(>|t|)` <- round(main_models_combined$`Pr(>|t|)`, digits = 3)
main_models_combined$Estimate <- round(main_models_combined$Estimate, digits = 3)
main_models_combined$`Cluster s.e.` <- round(main_models_combined$`Cluster s.e.`, digits = 3)
main_models_combined$`t value` <- round(main_models_combined$`t value`, digits = 3)
main_models_combined$Model <- "90m Spatial Resolution"
main_models_combined <- main_models_combined %>% select(Estimate:`Pr(>|t|)`, Model, Variable)

100*main_models_combined[1,1]/mean(na.omit(master_commune$PercentDeforest))

write.csv(main_models_combined, "Outputs/points_summary_output_combined.csv")

```

```{r Sensitivity Table}
sensitivity_table <- read.csv("Outputs/sensitivity_table.csv")
sensitivity_table <- sensitivity_table[,-1]
colnames(sensitivity_table) <- colnames(main_models)
sensitivity_table2 <- rbind(sensitivity_table, main_models)

write.csv(sensitivity_table2, "Outputs/sensitivity_table2.csv")

```