---
title: "01_Analysis_CombinedMain"
output: html_document
date: "2024-02-14"
---

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

Code and outputs for the main analysis for the buffer and PA area combined

```{r Load data and packages}
library(tidyverse);library(lfe)

master_buffer_combined <- read.csv("DataKeep/master_buffer_combined.csv")
master_buffer_combined_commune <- read.csv("DataKeep/master_buffer_combined_commune.csv")

```

```{r Main analysis}

master_buffer_combined$Tourism <- master_buffer_combined$Tourism/1000

#Model
fem_full_formula <- as.formula("PercentDeforest ~ Tourism + PopDensity + combined_precip + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full <- felm(fem_full_formula, data = master_buffer_combined)
summary_dd_reg_full <- summary(dd_reg_full)

#Model without precipitation
fem_full_formula_precip <- as.formula("PercentDeforest ~ Tourism + PopDensity + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_precip <- felm(fem_full_formula_precip, data = master_buffer_combined)
summary_dd_reg_full_precip <- summary(dd_reg_full_precip)

100*summary_dd_reg_full$coefficients[1,1]/mean(master_buffer_combined$PercentDeforest[master_buffer_combined$Year>2000])

#Output
models_total <- as.data.frame(summary_dd_reg_full$coefficients)
models_total$Analysis <- "PA"

models_total$Variable <- rownames(models_total)
models_total$Location <- "Inside PA + Buffer"
models_total$`Pr(>|t|)` <- round(models_total$`Pr(>|t|)`, digits = 3)
models_total$Estimate <- round(models_total$Estimate, digits = 3)
models_total$`Cluster s.e.` <- round(models_total$`Cluster s.e.`, digits = 3)
models_total$`t value` <- round(models_total$`t value`, digits = 3)

models_total$Variable[models_total$Variable=="PopDensity"] <- "Pop. Density"
models_total$Variable[models_total$Variable=="OilPriceAdj"] <- "Oil Price"
models_total$Variable[models_total$Variable=="Tourism1"] <- "Tourism"
models_total$Variable[models_total$Variable=="combined_precip"] <- "Pop. Density"
models_total$Variable[models_total$Variable=="Precipitation1"] <- "Precipitation"
models_total$Variable[models_total$Variable=="OilPriceAdj1"] <- "Oil Price"

write.csv(models_total, "Outputs/regression_table_total.csv")

#Output no precipitation
summary_dd_reg_full_buffer_precip <- summary_dd_reg_full_precip
combined_precip <- as.data.frame(summary_dd_reg_full_buffer_precip$coefficients)
combined_precip$Analysis <- "PA"
combined_precip$Variable <- rownames(combined_precip)
combined_precip$Location <- "Inside PA + Buffer 3km"
combined_precip$`Pr(>|t|)` <- round(combined_precip$`Pr(>|t|)`, digits = 5)
combined_precip$Estimate <- round(combined_precip$Estimate, digits = 5)
combined_precip$`Cluster s.e.` <- round(combined_precip$`Cluster s.e.`, digits = 5)
combined_precip$`t value` <- round(combined_precip$`t value`, digits = 5)
combined_precip$Variable[combined_precip$Variable=="PopDensity"] <- "Pop. Density"
combined_precip$Variable[combined_precip$Variable=="OilPriceAdj"] <- "Oil Price"
combined_precip$Variable[combined_precip$Variable=="Tourism1"] <- "Tourism"
combined_precip$Variable[combined_precip$Variable=="PopDensityCommune"] <- "Pop. Density"

write.csv(combined_precip, "Outputs/regression_table_total_precip.csv")

#with funding
master_combined_fund <- master_buffer_combined %>% filter(Year <=2012& Year >= 2005, Name != "Anjanaharibe-Sud", Name != "Pic d'Ivohibe")
fem_combined_formula_fund <- as.formula("PercentDeforest ~ Tourism + PopDensity + combined_precip + OilPriceAdj+ funding| Year + Name | 0 | Regions")
dd_reg_combined_fund <- felm(fem_combined_formula_fund, data = master_combined_fund)
summary_dd_reg_combined_fund  <- summary(dd_reg_combined_fund )

#no funding 
fem_combined_formula_nofund <- as.formula("PercentDeforest ~ Tourism + PopDensity + combined_precip + OilPriceAdj| Year + Name | 0 | Regions")
dd_reg_combined_nofund <- felm(fem_combined_formula_nofund, data = master_combined_fund)
summary_dd_reg_combined_nofund  <- summary(dd_reg_combined_nofund )

```

```{r Commune analysis}
#Model
fem_full_formula <- as.formula("PercentDeforest ~ Tourism + PopDensityCommune + Combined_commune_precip + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full <- felm(fem_full_formula, data = master_buffer_combined_commune)
summary_dd_reg_full <- summary(dd_reg_full)


100*summary_dd_reg_full$coefficients[1,1]/mean(master_buffer_combined_commune$PercentDeforest[master_buffer_combined_commune$Year>2000])

#Output
models_total <- as.data.frame(summary_dd_reg_full$coefficients)
models_total$Analysis <- "Commune"

models_total$Variable <- rownames(models_total)
models_total$Location <- "Inside PA + Buffer 3km"
models_total$`Pr(>|t|)` <- round(models_total$`Pr(>|t|)`, digits = 3)
models_total$Estimate <- round(models_total$Estimate, digits = 3)
models_total$`Cluster s.e.` <- round(models_total$`Cluster s.e.`, digits = 3)
models_total$`t value` <- round(models_total$`t value`, digits = 3)

models_total$Variable[models_total$Variable=="PopDensity"] <- "Pop. Density"
models_total$Variable[models_total$Variable=="OilPriceAdj"] <- "Oil Price"
models_total$Variable[models_total$Variable=="Tourism1"] <- "Tourism"
models_total$Variable[models_total$Variable=="PopDensityCommune"] <- "Pop. Density"
models_total$Variable[models_total$Variable=="Combined_commune_precip"] <- "Precipitation"
models_total$Variable[models_total$Variable=="OilPriceAdj1"] <- "Oil Price"

write.csv(models_total, "Outputs/regression_table_total_commune.csv")

```