TourismDeforestationSubmission / 01b_Analysis_BufferMain.Rmd
01b_Analysis_BufferMain.Rmd
Raw
---
title: "01_Analysis_BufferMain"
output: html_document
date: "2024-02-13"
---

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

Code and output for the main analysis in the buffer zone for the 3km, 6km, and 9km buffer zones

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

master_buffer_analysis_3km <- read.csv( "DataKeep/master_buffer_analysis_3km.csv")
master_buffer_analysis_commune_3km <- read.csv( "DataKeep/master_buffer_analysis_commune_3km.csv")
master_buffer_analysis_6km <- read.csv("DataKeep/master_buffer_analysis_6km.csv")
master_buffer_analysis_commune_6km <- read.csv("DataKeep/master_buffer_analysis_commune_6km.csv")
master_buffer_analysis_9km <- read.csv("DataKeep/master_buffer_analysis_9km.csv")
master_buffer_analysis_commune_9km <- read.csv("DataKeep/master_buffer_analysis_commune_9km.csv")
```

```{r 3Km Buffer}
#Full analysis
fem_full_formula_3km <- as.formula("PercentDeforest ~ Tourism + PopDensity + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_3km <- felm(fem_full_formula_3km, data = master_buffer_analysis_3km)
summary_dd_reg_full_3km <- summary(dd_reg_full_3km)

#Main analysis no precipitation

fem_full_formula_3km_precip <- as.formula("PercentDeforest ~ Tourism + PopDensity + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_3km_precip <- felm(fem_full_formula_3km_precip, data = master_buffer_analysis_3km)
summary_dd_reg_full_3km_precip <- summary(dd_reg_full_3km_precip)

#Commune analysis
fem_full_formula_commune_3km <- as.formula("PercentDeforest ~ Tourism + PopDensityCommune + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_commune_3km <- felm(fem_full_formula_commune_3km, data = master_buffer_analysis_commune_3km)
summary_dd_reg_full_commune_3km <- summary(dd_reg_full_commune_3km)

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

#Output compilation
models_3km <- as.data.frame(summary_dd_reg_full_3km$coefficients)
models_3km$Analysis <- "PA"
models_3km_commune <- as.data.frame(summary_dd_reg_full_commune_3km$coefficients)
models_3km_commune$Analysis <- "Commune"
models_3km  <- rbind(models_3km, models_3km_commune)

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

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

100*models_3km[1,1]/mean(master_buffer_analysis_commune_3km$PercentDeforest[master_buffer_analysis_commune_3km$Year>2000])

write.csv(models_3km, "Outputs/regression_table_3km.csv")


models_3km_precip <- as.data.frame(summary_dd_reg_full_3km_precip$coefficients)
models_3km_precip$Analysis <- "PA"

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

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

write.csv(models_3km_precip, "Outputs/regression_table_3km_precip.csv")


```

```{r 6Km Buffer}
#Main analysis

fem_full_formula_6km <- as.formula("PercentDeforest ~ Tourism + PopDensity + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_6km <- felm(fem_full_formula_6km, data = master_buffer_analysis_6km)
summary_dd_reg_full_6km <- summary(dd_reg_full_6km)


#Commune analysis

fem_full_formula_commune_6km <- as.formula("PercentDeforest ~ Tourism + PopDensityCommune + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_commune_6km <- felm(fem_full_formula_commune_6km, data = master_buffer_analysis_commune_6km)
summary_dd_reg_full_commune_6km <- summary(dd_reg_full_commune_6km)

#Output compilation
models_6km <- as.data.frame(summary_dd_reg_full_6km$coefficients)
models_6km$Analysis <- "PA"
models_6km_commune <- as.data.frame(summary_dd_reg_full_commune_6km$coefficients) 
models_6km_commune$Analysis <- "Commune"
models_6km  <- rbind(models_6km, models_6km_commune)

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

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

write.csv(models_6km, "Outputs/regression_table_6km.csv")

```

```{r 9km Buffer}
#Main analysis

fem_full_formula_9km <- as.formula("PercentDeforest ~ Tourism + PopDensity + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_9km <- felm(fem_full_formula_9km, data = master_buffer_analysis_9km)
summary_dd_reg_full_9km <- summary(dd_reg_full_9km)

#Commune analysis
fem_full_formula_commune_9km <- as.formula("PercentDeforest ~ Tourism + PopDensityCommune + Precipitation + OilPriceAdj | Year + Name | 0 | Regions")
dd_reg_full_commune_9km <- felm(fem_full_formula_commune_9km, data = master_buffer_analysis_commune_9km)
summary_dd_reg_full_commune_9km <- summary(dd_reg_full_commune_9km)

#Output Compilation
models_9km <- as.data.frame(summary_dd_reg_full_9km$coefficients)
models_9km$Analysis <- "PA"
models_9km_commune <- as.data.frame(summary_dd_reg_full_commune_9km$coefficients) 
models_9km_commune$Analysis <- "Commune"
models_9km  <- rbind(models_9km, models_9km_commune)

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

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

write.csv(models_9km, "Outputs/regression_table_9km.csv")


```