{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Methodology for the Emissions Factors Calculation: Spatial Panel Data Modeling"},{"metadata":{},"cell_type":"markdown","source":"\n*image source*: https://www.economist.com/united-states/2017/10/19/the-story-of-puerto-ricos-power-grid-is-the-story-of-puerto-rico"},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"# Load libraries\nlibrary(tidyverse)\nlibrary(jsonlite)\nlibrary(raster)\nlibrary(sp)\nlibrary(rgdal)\nlibrary(maptools)\nlibrary(rgeos)\nlibrary(ggplot2)\nlibrary(dplyr)\nlibrary(SpatioTemporal)\nlibrary(Matrix)\nlibrary(plotrix)\nlibrary(maps)\nlibrary(shadowtext)\nlibrary(sf)\nlibrary(ggmap)\nlibrary(plyr)\nlibrary(osmar)\nlibrary(geosphere)\nlibrary(osmdata)\nlibrary(repr)\nlibrary(spdep)\nlibrary(tidyr)\nlibrary(splm)\nlibrary(lubridate)\nlibrary(multcomp)\nlibrary(magrittr)\nlibrary(knitr)\nlibrary(kableExtra)\nlibrary(ggpubr)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this notebook we suggest an approach to compute the Average and Marginal Emissions Factors (AEF and MEF) for the electricity generation and demonstrate this on the example of Puerto Rico, using novel data collected through remote sensing techniques. Driven by the goal of boosting the accuracy of the AEF and MEFs calculations, we employ a method called **Spatial Panel Modeling**. This method can be leveraged to record a spatial variability of the data **both across spatial units and over time** (Millo and Piras 2012), hence, tangibly increase the precision of the estimations. Moreover, the approach allows us to easily account for other conditions (like weather) inside the model and make inferences about the factors that are capable of modifying the EFs estimates."},{"metadata":{},"cell_type":"markdown","source":"# Executive summary\n\n* This mini-research uses Spatial Panel Data Modeling to compute the AEF and MEF for Puerto Rico.\n\n* The method is profitable as it allows extracting every piece of information from the variability of data in space and over time and accounts for multiple sources of NO2; this is coupled with simplicity of computations and interpretation.\n\n* The study reports the AEF of 0.5 kg/MWh and MEF of 0.07 kg/MWh.\n\n* It also shows a way to calculate the AEF by months.\n\n* Averaging of the obtained monthly values leads to a yearly AEF of 0.73 kg/MWh - an estimate which is a little greater than the one obtained in the general model; this implies focusing on smaller periods may favor higher precision of the estimates.\n\n* The analysis will benefit from using data on NO2 emission and energy generation sliced by smaller periods, for example, days/hours.\n\n* The suggested methodology is characterized by flexibility of specification, easy scalability, and has an increased potential for defining emissions factors of power stations with a substantial NO2 emission rate.\n\nThe kernel is organized as follows. We start from the data description and proceed with the exploratory analysis. Next, a detailed method explanation is provided, key findings on the emissions factors calculation in Puerto Rico are reported, and computational nuances are presented. We then outline recommendations for the application of this methodology for another geospatial area, elaborate on the important assumptions to be made, the advantages and shortcomings of the approach, provide a comparative perspective with the current methodology, and conclude mentioning some other possibilities in terms of the data to be utilized."},{"metadata":{},"cell_type":"markdown","source":"# Data description\n\n### Dependent Variable\nFor the tasks set, nitrogen dioxide measured in mol/m^2 was converted into kilograms via multiplying the initial values by the number of m^2 in a rectangular raster cell and by 46.006 g/mol -- the molar mass of NO2.\n\n### Predictors\nThe central predictors of our analysis were represented by a set of variables, reflecting the amount of electricity, generated by a respective station, measured in megawatt-hours (MWh). We hence separate the electricity production between each plant since there is a substantive heterogeneity in the ways each plant pollutes the atmosphere relative to the amount of electricity it produces. In other words, some stations might emit harmful substances a lot due to the high amount of energy generated, whereas others can release chemicals just because they are old, ineffective, or exploit a different fuel. Therefore, pulling all the megawatt-hours into one variable might result in misleading findings, and prohibit from capturing some dependencies in the data. We explored only thermal power stations (oil, gas, coal) and excluded 2 plants from the analysis -- Vieques and HIMA San Pablo Caguas -- since their electricity production was negligible compared to the rest of the stations.\n\nImportantly, we used an alternative dataset on electricity production in 2018 provided by U.S. Energy Information Administration (see the references), intending to obtain the data that can change over time differently for each plant. These data are similar to the ones which were offered to us, except for the variability they have over time. Since Puerto Rico was not present in the 2019 data, we obtained the estimated electricity generation for January 2019 - June 2019 by scaling the respective data in 2018 by the ratio between the electricity consumption (see the references) in 2019 and 2018 (see the code for more details).\n\n### Control Variables and Confounding Factors\n\n* Population density *(GPWv411: Basic Demographic Characteristics)*\n\n* Motorways (binary variable, presence or absence in a raster cell) *(OpenStreetMap)*\n\n* Wind velocity calculated using U and V wind components *(Global Forecast System)*\n\n* Temperature, Humidity, and Precipitation *(Global Forecast System)*\n\n* Pressure and Cloudiness *(Sentinel-5P OFFL NO2)*\n\nWe selected only 18:00 rasters from the Global Forecast System since NO2 was measured mostly in the evening."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"setwd('..')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"################################################## Plants Data\n\n# Load plants monthly data\ndf_plants_monthly <- data.frame(readxl::read_excel('input/puerto-rico-plants-monthly-data-eia/PR_mwh_monthly_data_EIA.xlsx'))\ndf_plants_monthly <- df_plants_monthly[,c(\"Plant.Id\", \"Plant.Name\", \"AER..Fuel.Type.Code\",'Netgen..January', 'Netgen..February', 'Netgen..March', 'Netgen..April', \n 'Netgen..May', 'Netgen..June', 'Netgen..July', 'Netgen..August', 'Netgen..September',\n 'Netgen..October', 'Netgen..November', 'Netgen..December')]\ncolnames(df_plants_monthly) <- str_replace(colnames(df_plants_monthly), \"Netgen..\", \"netgen_\")\n\n# Caclulate consumtion mwh difference in 2019 compare to 2018 based on consumption difference\ndf_pr_consumption <- data.frame(\n read.csv('input/ds4g-puerto-rico-electricity-consumption-by-month/EIA_puerto_rico_electric_consumption.csv'))\n\ndf_pr_consumption <- df_pr_consumption[,c('date', 'all_sectors')]\ndf_pr_consumption$date <- as.Date(df_pr_consumption$date)\n\nconsumption_jan_june_2018 <- df_pr_consumption[df_pr_consumption$date >= \"2018-01-01\" & df_pr_consumption$date <= \"2018-06-01\",]\nconsumption_jan_june_2019 <- df_pr_consumption[df_pr_consumption$date >= \"2019-01-01\" & df_pr_consumption$date <= \"2019-06-01\",]\n\nconsumption_difference <- consumption_jan_june_2019$all_sectors / consumption_jan_june_2018$all_sectors\nconsumption_difference <- c(rep(1, 6), consumption_difference)\n\n# Select and normalize monthly netgen\ndf_plants_monthly_netgen <- df_plants_monthly[,c(\"netgen_July\", \"netgen_August\", \"netgen_September\", \"netgen_October\", \"netgen_November\",\n \"netgen_December\", \"netgen_January\", \"netgen_February\", \"netgen_March\",\n \"netgen_April\", \"netgen_May\", \"netgen_June\")]\ndf_plants_monthly_netgen <- data.frame(mapply(`*`, df_plants_monthly_netgen, consumption_difference))\ndf_plants_monthly <- cbind.data.frame(df_plants_monthly[,c(\"Plant.Id\", \"Plant.Name\", \"AER..Fuel.Type.Code\")],\n df_plants_monthly_netgen)\n\n# Get summary GwH for each plant\nnames(df_plants_monthly)[1] <- 'ID'\nplants_id_name_type <- as.data.frame(df_plants_monthly[,c(\"ID\", \"Plant.Name\", \"AER..Fuel.Type.Code\")])\ndf_plants_monthly <- df_plants_monthly %>% \n group_by(ID) %>%\n summarise_if(is.numeric, sum)\ndf_plants_monthly <- as.data.frame(df_plants_monthly)\nplants_id_name_type <- plants_id_name_type[!duplicated(plants_id_name_type$ID),]\ndf_plants_monthly <- df_plants_monthly %>%\n left_join(plants_id_name_type, by = 'ID')\n\n# Read dataframe with coordinates of plants dataframe\ndf_plants_id_coordinates <- data.frame(readxl::read_excel('input/pr-plants-coordinates-eia/PR_plants_coordinates_EIA.xlsx'))\ndf_plants_id_coordinates <- df_plants_id_coordinates[!duplicated(df_plants_id_coordinates$Plant.ID),]\n# Add coordinates to new plants\ndf_plants_monthly <- left_join(df_plants_monthly, df_plants_id_coordinates)\n# Drop \"green\" plants\ndf_plants_monthly <- df_plants_monthly[df_plants_monthly$Technology %in% c('Natural Gas Fired Combined Cycle',\n 'Petroleum Liquids',\n 'Conventional Steam Coal'),]\n\n# Make the spatial points data frame object\nspdf_plants_monthly <- sp::SpatialPointsDataFrame(coords = data.frame(lon=as.numeric(df_plants_monthly$Longitude),\n lat=as.numeric(df_plants_monthly$Latitude)),\n data = df_plants_monthly, \n proj4string = CRS(\"+init=epsg:4326\"))\n\n############################################# S5P NO2 data\n\n# Create list with NO2 file names\nsp5no2_files <- list.files('input/ds4g-environmental-insights-explorer/eie_data/s5p_no2', full.names=T, pattern=\".tif$\")\n\n# Create list with all raster object\n# It takes some time\nsp5no2_list <- list()\nfor(i in 1:length(sp5no2_files)){\n raster_temp <- stack(sp5no2_files[i])\n # Adjust extent for all raster images\n extent(raster_temp) <- extent(-67.32431, -65.19082, 17.90464, 18.5694)\n sp5no2_list <- append(sp5no2_list, raster_temp)\n}\n\n\n# Get rasters from one band\ngetRastersFromBand_SP5NO2 <- function(main_list, band, fact){\n \n band_rasters_list <- list()\n for(i in 1:length(main_list)){\n band_rasters_list <- append(band_rasters_list, main_list[[i]][[band]])\n }\n band_rasters_stack <- stack(band_rasters_list)\n # Horizontal flip\n band_rasters_stack <- flip(band_rasters_stack, direction = \"y\")\n # Aggregate values\n band_rasters_stack <- aggregate(band_rasters_stack, fact=fact)\n return(band_rasters_stack)\n \n}\n\n# Create NO2 rasters\nno2_rasters <- getRastersFromBand_SP5NO2(sp5no2_list, band=1, fact=20)\nnames(no2_rasters) <- substr(sp5no2_files, start = 61, stop= nchar(sp5no2_files[[1]]))\nno2_rasters_spdf <- as(no2_rasters, \"SpatialGridDataFrame\")\n# Create pressure rasters\npressure_rasters <- getRastersFromBand_SP5NO2(sp5no2_list, band=5, fact=20)\nnames(pressure_rasters) <- substr(sp5no2_files, start = 61, stop= nchar(sp5no2_files[[1]]))\n# Create cloud fraction rasters\ncloud_fraction_rasters <- getRastersFromBand_SP5NO2(sp5no2_list, band=7, fact=20)\nnames(cloud_fraction_rasters) <- substr(sp5no2_files, start = 61, stop= nchar(sp5no2_files[[1]]))\n\n\ncreateDataframeFromRasters_SP5NO2 <- function(rasters, weekly=F, monthly=F, no2=F){\n \n ### Remove duplicates from one day\n layer_names <- sapply(names(rasters), FUN = function(x){strsplit(x, '_')[[1]][3]})\n layer_dates <- as.Date(layer_names,format=\"%Y%m%d\")\n # Convert to dataframe\n df <- as.data.frame(t(as.data.frame(rasters)))\n # Remove days which are duplicates and have missings\n df['dates'] <- layer_dates\n df['duplicate_dates'] <- as.logical(duplicated(layer_dates, fromLast=T) + duplicated(layer_dates, fromLast=F))\n df['missings'] <- apply(df, FUN=function(x){any(is.na(x))}, 1)\n df <- df[(as.integer(!df['duplicate_dates']) + as.integer(!df['missings'])) != 0,]\n # Remove other duplicated days\n df <- df[!df['duplicate_dates'],]\n # Remove redundant columns\n df <- df[1:(ncol(df)-2)]\n # Assign dates as row names\n rownames(df) <- df$dates\n \n # Aggregate by week\n if (monthly){\n \n year_month_indices <- format(df$dates, \"%Y-%m\")\n if (no2){\n df <- subset(df, select = -c(dates)) # FOR SUM NO2\n df <- aggregate(df, by=list(year_month_indices), FUN=sum, na.rm=T)\n } else {\n df <- aggregate(df, by=list(year_month_indices), FUN=mean, na.rm=T)\n }\n \n rownames(df) <- df$Group.1\n }\n \n if (weekly){\n \n df <- df[df$dates != \"2018-07-01\",]\n year_week_indices <- strftime(df$dates, format = \"%Y_%V\")\n if (no2){\n df <- subset(df, select = -c(dates)) # FOR SUM NO2\n df <- aggregate(df, by=list(year_week_indices), FUN=sum, na.rm=T)\n } else {\n df <- aggregate(df, by=list(year_week_indices), FUN=mean, na.rm=T)\n }\n \n rownames(df) <- df$Group.1\n }\n \n \n # Remove redundant columns\n df <- df[,startsWith(colnames(df), \"V\")]\n \n return(df)\n \n}\n\n# Create dataframes for no2, pressure and cloud fraction\nno2_df <- createDataframeFromRasters_SP5NO2(no2_rasters, weekly=T, no2=T)\npressure_df <- createDataframeFromRasters_SP5NO2(pressure_rasters, weekly=T)\ncloud_fraction_df <- createDataframeFromRasters_SP5NO2(cloud_fraction_rasters, weekly=T)\n\n###################################### GFS: Global Forecast System\n\n# Select only measurments at 18:00\ngfs_files <- list.files('input/ds4g-environmental-insights-explorer/eie_data/gfs', full.names=T, pattern=\"18.tif$\")\n\n### Create list with all raster objects and with all layers for each raster\ngfs_list<- list()\nfor(i in 1:length(gfs_files)){\n raster_temp <- stack(gfs_files[i])\n # Adjust extent for all raster images\n extent(raster_temp) <- extent(-67.32431, -65.19082, 17.90464, 18.5694)\n gfs_list <- append(gfs_list, raster_temp)\n}\n\n# Get rasters from one band\ngetRastersFromBand_GFS <- function(main_list, band){\n \n band_rasters_list <- list()\n for(i in 1:length(main_list)){\n band_rasters_list <- append(band_rasters_list, main_list[[i]][[band]])\n }\n band_rasters_stack <- stack(band_rasters_list)\n\n return(band_rasters_stack)\n \n}\n\n\n# Create temperature_2m_above_ground rasters\ntemperature_rasters <- getRastersFromBand_GFS(gfs_list, band=1)\nnames(temperature_rasters) <- substr(gfs_files, start = 61, stop=70)\n# Create relative_humidity_2m_above_ground rasters\nhumidity_rasters <- getRastersFromBand_GFS(gfs_list, band=3)\nnames(humidity_rasters) <- substr(gfs_files, start = 61, stop=70)\n# Create u_component_of_wind_10m_above_ground rasters\nu_wind_rasters <- getRastersFromBand_GFS(gfs_list, band=4)\nnames(u_wind_rasters) <- substr(gfs_files, start = 61, stop=70)\n# Create v_component_of_wind_10m_above_ground rasters\nv_wind_rasters <- getRastersFromBand_GFS(gfs_list, band=5)\nnames(v_wind_rasters) <- substr(gfs_files, start = 61, stop=70)\n# Create precipitable_water_entire_atmosphere rasters\nprecipitable_water_rasters <- getRastersFromBand_GFS(gfs_list, band=6)\nnames(precipitable_water_rasters) <- substr(gfs_files, start = 61, stop=70)\n\n\ncreateDataframeFromRasters_GFS <- function(rasters, target_raster, weekly=F, monthly=F){\n \n # Reproject GFS rasters according to NO2\n rasters <- resample(rasters, target_raster)\n # Create dataframe from rastet \n df <- as.data.frame(t(as.data.frame(rasters)))\n # Convert row names to dates\n df$dates <- as.Date(substr(rownames(df), 2, 9), format=\"%Y%m%d\")\n \n # Aggregate by week\n if (weekly){\n \n df <- df[df$dates != \"2018-07-01\",]\n year_week_indices <- strftime(df$dates, format = \"%Y_%V\")\n df <- aggregate(df, by=list(year_week_indices), FUN=mean, na.rm=T)\n df <- df[df$Group.1 != \"2018_01\",]\n rownames(df) <- df$Group.1\n \n }\n \n if (monthly){\n \n # df <- df[df$dates != \"2018-07-01\",]\n year_month_indices <- format(df$dates, \"%Y-%m\")\n df <- aggregate(df, by=list(year_month_indices), FUN=mean, na.rm=T)\n \n rownames(df) <- df$Group.1\n }\n \n # Remove redundant columns\n df <- df[,startsWith(colnames(df), \"V\")]\n \n return(df)\n \n}\n\n# Create dataframes for weather conditions\ntemperature_df <- createDataframeFromRasters_GFS(temperature_rasters, no2_rasters, weekly=T)\nhumidity_df <- createDataframeFromRasters_GFS(humidity_rasters, no2_rasters, weekly=T)\nu_wind_df <- createDataframeFromRasters_GFS(u_wind_rasters, no2_rasters, weekly=T)\nv_wind_df <- createDataframeFromRasters_GFS(v_wind_rasters, no2_rasters, weekly=T)\nprecipitable_water_df <- createDataframeFromRasters_GFS(precipitable_water_rasters, no2_rasters, weekly=T)\n\n################################## Coordinates of cells\ncells_coordinates_df <- as.data.frame(coordinates(no2_rasters))\ncells_coordinates_df$ID <- paste0(\"V\", 1:nrow(cells_coordinates_df), sep = \"\")\n\n\n#################################### Plants: Cell ID's and generation gwh\n\n# Extract NO2 values corresponding to each plant\nno2_plants <- as.data.frame(raster::extract(no2_rasters, spdf_plants_monthly, cellnumbers=T,\n sp = T))\n\n# Add gwh plant names and operating name from plants\nplants_df <- no2_plants[,c(names(no2_plants)[startsWith(names(no2_plants), \"netgen\")],\n \"Plant.Name\", \"Operating.Year\", \"lon\", \"lat\")]\n\n# Add \"V\" to ID\nplants_df$ID <- sapply(no2_plants$cells, function(x){paste0(\"V\", as.character(x))})\n\n# Remove plants with 0 gwh\nplants_df <- plants_df[apply(plants_df[,1:(ncol(plants_df)-5)], 1, sum) > 0,]\n\n# Remove word \"plant\"\nplants_df$Plant.Name <- str_replace(plants_df$Plant.Name, pattern = \" Plant\", replacement = \"\")\n\n# Replace \" \" to \"_\"\nplants_df$Plant.Name <- str_replace_all(plants_df$Plant.Name, pattern = \" \", replacement = \"_\")\nplants_df$Plant.Name <- str_replace_all(plants_df$Plant.Name, pattern = \"_-_\", replacement = \"_\")\n\n# Drop plants outside the map\nplants_monthly_df <- plants_df[plants_df$ID != \"VNA\",]\n\n\n################################################ GPWv411: Basic Demographic Characteristics\n\n# Load demographic data from GEE\ndem_raster <- raster(\"input/pr-population-data-gpwv411/ee_demographic_data.tif\", band=1)\n# Reproject demographic data according to NO2\nextent(dem_raster) <- extent(-67.32431, -65.19082, 17.90464, 18.5694)\ndem_raster <- resample(dem_raster, no2_rasters)\n# Extract values from demographic raster\ndem_df <- as.data.frame(dem_raster)\n# Create ID\ndem_df$ID <- paste0(\"V\", 1:nrow(dem_df), sep = \"\")\ncolnames(dem_df) <- c(\"population\", \"ID\")\n# Fill NA with 0\ndem_df$population[is.na(dem_df$population)] <- 0\n\n\n################################################ OpenStreetMaps: Roads\n\n# osmdata doesn't work correctly in Kaggle notebook\n\n# Parse motorway, trunk and primary roads from Puerto Rico\n# pr_highways <- opq(bbox = 'puerto rico') %>%\n # add_osm_feature(key = 'highway', value = c('motorway')) %>%\n # osmdata_sp()\n# Create SpatialLinesDataFrame\n# highways_spldf <- pr_highways$osm_lines\n\n# Load saved SpatialLinesDataFrame, which can be obtained by running the code above\nhighways_spldf <- readRDS(\"input/puerto-rico-highways/highways_spldf.rds\")\n\n# Get cells coordinates\ncells_coordinates <- data.frame(lon=cells_coordinates_df$y,\n lat=cells_coordinates_df$x)\n\n# Overlap SpatialPoints to SpadialGrid and get road indicator and cell ID\nhighways_df <- over(no2_rasters_spdf, highways_spldf)$highway\nhighways_df <- data.frame('ID'=paste0(\"V\", 1:length(highways_df), sep = \"\"),\n 'highway'=highways_df)\n\nhighways_df$highway <- ifelse(is.na(highways_df$highway), 0, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"################# LOAD THE PREPARED DATA\n\n# cells_coordinates_df <- readRDS('input/saved-dataframes/cells_coordinates_df.rds')\n# cloud_fraction_df <- readRDS('input/saved-dataframes/cloud_fraction_df.rds')\n# dem_df <- readRDS('input/saved-dataframes/dem_df.rds')\n# elevation_df <- readRDS('input/saved-dataframes/elevation_df.rds')\n# highways_df <- readRDS('input/saved-dataframes/highways_df.rds')\n# humidity_df <- readRDS('input/saved-dataframes/humidity_df.rds')\n# no2_df <- readRDS('input/saved-dataframes/no2_df.rds')\n# plants_monthly_df <- readRDS('input/saved-dataframes/plants_monthly_df.rds')\n# precipitable_water_df <- readRDS('input/saved-dataframes/precipitable_water_df.rds')\n# pressure_df <- readRDS('input/saved-dataframes/pressure_df.rds')\n# temperature_df <- readRDS('input/saved-dataframes/temperature_df.rds')\n# u_wind_df <- readRDS('input/saved-dataframes/u_wind_df.rds')\n# v_wind_df <- readRDS('input/saved-dataframes/v_wind_df.rds')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Exploratory Data Analysis\n"},{"metadata":{},"cell_type":"markdown","source":"There are three central factors that can be associated with the NO2 emission: power plants, motorways (intercity roads), and population (within-city traffic and some other activities of a human). The graph below demonstrates the geographical placement of the mentioned causes of the NO2 pollution in Puerto Rico. In order to uncover the contribution of a separate factor (power plants in our case), we need to come up with a model that can effectively account for the existence of the rest."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Adjust plotting window\noptions(repr.plot.width=16, repr.plot.height = 8)\n\n# Create Puerto Rico bbox\nbbox_pr <- make_bbox(lon=c(-67.32431, -65.19082), lat=c(17.89566, 18.5694))\n# Download map of Puerto Rico\nmap_pr <- get_stamenmap(bbox_pr,\n maptype=\"terrain-background\")\n\n# Function for transparent visualization of rasters\ngplot_data <- function(x, maxpixels = 50000) {\n x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)\n coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))\n ## Extract values\n dat <- utils::stack(as.data.frame(raster::getValues(x)))\n names(dat) <- c('value', 'variable')\n \n dat <- dplyr::as.tbl(data.frame(coords, dat))\n \n if (!is.null(levels(x))) {\n dat <- dplyr::left_join(dat, levels(x)[[1]],\n by = c(\"value\" = \"ID\"))\n }\n dat\n}\n\n# Load population raster\ndem_raster_viz <- raster(\"input/pr-population-data-gpwv411/ee_demographic_data.tif\", band=1)\nextent(dem_raster_viz) <- extent(-67.32431, -65.19082, 17.90464, 18.5694)\n# Replace low values with NA\ntemp_values <- values(dem_raster_viz$ee_demographic_data)\ntemp_values[is.na(temp_values)] <- 0\ntemp_values[temp_values < 50] <- NA\nvalues(dem_raster_viz) <- temp_values\n# Create population dataframe for geom_tile\ndem_raster_data <- gplot_data(dem_raster_viz, maxpixels = 10000) # Choose number of pixels\n# Create highways dataframe for geom_tile\nhighways_spldf_f <- fortify(highways_spldf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggmap(map_pr) +\n geom_tile(data = dem_raster_data, aes(x, y, fill = value)) +\n scale_fill_gradient(low = \"gray75\",\n high = \"black\",\n na.value=\"transparent\", name = \"Population\") + \n geom_path(data=highways_spldf_f, aes(x=long, y=lat, group=group), size=1, colour = \"gray25\", alpha=0.3) + \n geom_point(aes(lon, lat),\n data=plants_df, shape=19, size=4, colour= \"brown3\") + \n geom_shadowtext(data = plants_df,\n aes(x = lon, y = lat, label = str_replace_all(Plant.Name, '_', ' ')), \n size = 5, vjust = 0.1, hjust = -0.1, family=\"Helvetica\", fontface = \"bold\",\n bg.colour='white', color=\"black\") + \n annotate(\"text\", label = \"San Juan\", x = -66.1, y = 18.52, size = 8, colour = \"gray50\") + \n ggtitle(\"Population, Motorways, and Power Plants in Puerto Rico\") + \n theme(plot.title = element_text(size = 32, face = \"bold\", hjust = 0.5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, let us glance at how the average yearly amount of NO2 is spread over the sub-national region. From what can be captured by an eye, it is noticeable that NO2 is mostly concentrated around the places with high population density. Some other vivid NO2 locations are barely distinguishable. Let us then investigated a time-related dimension and see what can be extracted from it."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Create NO2 raster for visualization\nno2_rasters_viz <- getRastersFromBand_SP5NO2(sp5no2_list, band=1, fact=5)\nno2_rasters_viz <- calc(no2_rasters_viz, mean, na.rm=T)\nno2_rasters_viz <- disaggregate(no2_rasters_viz, 5, method='bilinear')\n# Create NO2 dataframe for geom_tile\nno2_raster_data <- gplot_data(no2_rasters_viz, maxpixels = 10000) # Choose number of pixels","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggmap(map_pr) +\n geom_tile(data = no2_raster_data, aes(x, y, fill = value), alpha = 0.8) +\n scale_fill_gradient(low =\"white\", high = \"firebrick1\", name = \"NO2 mol/m^2\") + \n geom_path(data=highways_spldf_f, aes(x=long, y=lat, group=group), size=1, colour = \"gray25\", alpha=0.1) + \n geom_point(aes(lon, lat),\n data=plants_df, shape=19, size=3, alpha=0.3, colour= \"gray25\") +\n annotate(\"text\", label = \"San Juan\", x = -66.1, y = 18.52, size = 8, colour = \"black\", fontface =\"bold\") + \n ggtitle(\"Average Yearly NO2 Emission in Puerto Rico\") + \n theme(plot.title = element_text(size = 32, face = \"bold\", hjust = 0.5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The picture below represents the average NO2 in spatial units with power plants (red line) and NO2 in each spatial units without power plants (grey lines) over time. Interestingly, the red line does not behave differently from the grey ones, in fact, the former describes almost an averaging pattern for the latter. Such an inspection combined with a previous one tells us that neither spatial nor time-related dimensions separately are powerful enough to capture some systematic NO2 changes attributable to the electricity production by plants. **Therefore, we need an approach that encompasses both spatial and time-based variability of the data and is capable of properly treating other factors that bring about NO2 pollution.** As has been stated previously, to reach this goal we employ a method called Spatial Panel Data Modeling."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Create raster for NO2 visualization\nno2_df_viz <- createDataframeFromRasters_SP5NO2(no2_rasters, weekly=F, no2=T)\nno2_df_viz[no2_df_viz > quantile(no2_df_viz, 0.99, na.rm=T)] <- NA\nno2_df_viz[no2_df_viz < quantile(no2_df_viz, 0.01, na.rm=T)] <- NA\n\n# Split NO2 data for cells with plants without plants and with population\nno2_df_with_plants <- no2_df_viz[,unique(plants_df$ID)]\nno2_df_no_plants <- no2_df_viz[,colnames(no2_df)[!colnames(no2_df) %in% unique(plants_df$ID)]]\nno2_df_population <- no2_df_viz[,dem_df[dem_df$population > summary(dem_df$population)[5],]$ID]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Plot timeseries for all cells and means for plants\nts.plot(data.frame(no2_df_no_plants, apply(no2_df_with_plants, 1, mean, na.rm=T)),\n gpars=list(col=c(rep(adjustcolor( \"lightgrey\", alpha.f = 0.5), ncol(no2_df_no_plants)),\n adjustcolor(\"brown1\", alpha.f = 0.5)), xaxt=\"n\",\n xlab=\"Date\", ylab=\"NO2 mol/m^2\", main='NO2 in Spatial Units with and without Plants over Time', cex.main=2), \n lwd = c(rep(0.5, ncol(no2_df_no_plants)), 2))\n\nlegend(\"topleft\", bty=\"n\", lty=c(1,1), lwd=c(2,2), col=c(\"brown1\",\"lightgrey\"),\n legend=c(\"Average NO2 in Spatial Units with Plants\", \"NO2 in Spatial Units without Plants\"), cex = 1.5)\n\naxis(1, at=seq(1,length(rownames(no2_df_viz)), 30), labels=rownames(no2_df_viz)[seq(1,length(rownames(no2_df_viz)), 30)])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Methods\n\nThe **spatial element** in spatial panel data modeling can be addressed in three different ways (Millo and Piras 2012).\n\n**1. The first one - the spatial lag panel data model - is denoted by:**\n\n\\begin{equation}\ny=\\lambda\\left(\\mathbf{I}_{T} \\otimes \\mathbf{W}_{N}\\right) y+X \\beta+\\left(l_{T} \\otimes \\mathbf{I}_{N}\\right) \\mu+\\varepsilon\n\\end{equation}\n\n* where $\\otimes$ reflects Kronecker product,\n\n* $y$ is an $N_T \\times 1$ vector of observations on the outcome variable (NO2 in our case),\n\n* $X$ is an $N_T \\times k$ matrix of exogenous independent variables (estimated power generation in our case),\n\n* $\\lambda$ is the spatial parameter,\n\n* $\\beta$ is a $k \\times 1$ vector of the regression coefficients (EFs in our specification),\n\n* $N$ and $T$ are the number of spatial units and periods respectively,\n\n* $I_T$ is an identity matrix of dimension $T$,\n\n* $W_N$ is a $N \\times N$ spatial weight matrix with diagonal elements set to zero,\n\n* $k$ is the number of dependent (explanatory) variables,\n\n* $I_N$ is an $N \\times N$ identity matrix,\n\n* $l_T$ is a $T × 1$ vector of ones,\n\n* $\\mu$ is a vector of time-invariant individual effects that are not spatially\nautocorrelated,\n\n* $\\varepsilon$ is a random disturbance term.\n\n**2. The second - the spatial error panel data model - has the following form:**\n\n\\begin{equation}\n\\begin{array}{l}\ny=X \\beta+\\left(l_{T} \\otimes \\mathbf{I}_{N}\\right) \\mu+\\varepsilon \\\\\n\\varepsilon=\\rho\\left(\\mathbf{I}_{T} \\otimes \\mathbf{W}_{N}\\right) \\varepsilon+v\n\\end{array}\n\\end{equation}\n\n* where $\\rho$ is the spatial autoregressive parameter to be estimated,\n\n* $v$ is an $N_T \\times 1$ vector of random errors with the rest of notations having the above-denoted meaning.\n\n**3. The third one simultaneously incorporates a spatial lag of the dependent variable and a spatially autocorrelated error:**\n\n\\begin{equation}\n\\begin{array}{l}\ny=\\lambda\\left(\\mathbf{I}_{T} \\otimes \\mathbf{W}_{1}\\right) y+\\mathbf{X} \\boldsymbol{\\beta}+\\left(l_{T} \\otimes \\mathbf{I}_{N}\\right) \\boldsymbol{\\mu}+\\boldsymbol{\\varepsilon} \\\\\n\\boldsymbol{\\varepsilon}=\\rho\\left(\\mathbf{I}_{T} \\otimes \\mathbf{W}_{2}\\right) \\boldsymbol{\\varepsilon}+\\boldsymbol{v}\n\\end{array}\n\\end{equation}\n\nwhere $W_1$ and $W_2$ (as $W$ previously) are $N \\times N$ spatial weight matrices.\n\nThe __*panel (or longitudinal) element*__ of spatial panel data models can include the analysis of either random or fixed effects. The current research uses random effects models since our key control variable of focus in the context of NO2 investigation - population density - is assumed to be constant over the available period; otherwise, variables that do not vary over time are wiped out from the model.\n\nTo investigate spatial panel models we exploited a generalized moments (GM) technique due to its relative computational easiness compared to the maximum likelihood (ML) estimator.\n\n\nWe follow Baltagi et al. (2003) to figure out if it is necessary to incorporate spatial error correlation in the model estimation. Lagrange multiplier test (see below) shows a significant LM*-lambda statistic, meaning that spatial error correlation is supported by the data. Moreover, we assume that spatial lag dependence - a dependence of an outcome variable on its values in closest geographical locations - is also essential to account for in the analysis of harmful substances emission.\n\n**Therefore, the current study implements the third model that simultaneously captures a spatial lag of $y$ and a spatially autocorrelated error.**"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"##### Create a spatial matrix for observations\n# Spatial weights matrix\ncoords <- coordinates(cells_coordinates_df[, 1:2])\n# KNN\nnb <- knn2nb(knearneigh(coords, k = 8), row.names = cells_coordinates_df[,3])\nspatial_mat_list <- nb2listw(nb)\n\n# Removing outliers for NO2\nno2_df[no2_df > quantile(as.matrix(no2_df), 0.99)] <- quantile(as.matrix(no2_df), 0.99)\n\n##### Convert to long format\n\n# Define the number of spatial units\nn_spatial_units <- ncol(temperature_df)\n\n# Function for convertion to long format\nShortToLong <- function(df, n, varname, Xlead = FALSE){\n df <- eval(parse(text = paste(\"cbind.data.frame(gather(as.data.frame(df), ID,\",\n varname, \"), 'time' = rep(rownames(df), n))\")))\n df$time <- as.character(df$time)\n if (Xlead == TRUE){\n # if leading X appears in a time variable\n df$time <- as.Date(substr(df$time, 2, nchar(df$time)), format=\"%Y%m%d\")\n }\n else{\n df$time <- df$time\n }\n df$time <- as.character(df$time)\n return(df)\n}\n\n# Applying to the dataframes\ntemperature_df_long <- ShortToLong(temperature_df, n_spatial_units, 'temp', Xlead = F)\nhumidity_df_long <- ShortToLong(humidity_df, n_spatial_units, 'humid', Xlead = F)\nu_wind_df_long <- ShortToLong(u_wind_df, n_spatial_units, 'uwind', Xlead = F)\nv_wind_df_long <- ShortToLong(v_wind_df, n_spatial_units, 'vwind', Xlead = F)\nprecipitable_water_df_long <- ShortToLong(precipitable_water_df, n_spatial_units, 'precip', Xlead = F)\nno2_df_long <- ShortToLong(no2_df, n_spatial_units, 'no2')\npressure_df_long <- ShortToLong(pressure_df, n_spatial_units, 'pressure')\ncloud_fraction_df_long <- ShortToLong(cloud_fraction_df, n_spatial_units, 'cloud')\n\n# Merge plants which were placed in one cell\nplants_monthly_df_ <- plants_monthly_df %>%\n group_by(ID) %>%\n summarise_at(c(\"netgen_July\", \"netgen_August\", \"netgen_September\", \"netgen_October\", \n \"netgen_November\", \"netgen_December\", \"netgen_January\", \"netgen_February\", \n \"netgen_March\", \"netgen_April\", \"netgen_May\", \"netgen_June\" ), sum)\n\n# Renaming 'joint' plants to understand that they actually represent several plants, not one\nplants_monthly_df[plants_monthly_df$ID == 'V158', 'Plant.Name'] <- 'AESPR__Jobos' \nplants_monthly_df[plants_monthly_df$ID == 'V151', 'Plant.Name'] <- 'EcoEl__Costa_Sur'\nplants_monthly_df[plants_monthly_df$ID == 'V38', 'Plant.Name'] <- 'Central_SJ__Palo_Seco'\n\n# Remove duplicates by ID\nplants_monthly_df_ <- plants_monthly_df_[!duplicated(plants_monthly_df_$ID), ]\n#saveRDS(plants_monthly_df_, 'plants_monthly_df_merged.rds')\n\n# Converting to long format the data on plants\nplants_df <- gather(plants_monthly_df_, ID, 'estimated_generation_mwh')\ncolnames(plants_df)[1] <- 'month'\nPlant.Name_ID <- plants_monthly_df[, c('Plant.Name', 'ID')]\nPlant.Name_ID <- Plant.Name_ID[!duplicated(Plant.Name_ID$ID),] %>% arrange(ID)\nplants_df <- cbind.data.frame(plants_df, 'ID' = rep(Plant.Name_ID$ID, 12),\n 'Plant.Name' = rep(Plant.Name_ID$Plant.Name, 12)) \n\n\n# Month as a number\ndates <- cbind.data.frame('month' = unique(plants_df$month), 'inx' = c(7:12, 1:6))\nfor (i in unique(plants_df$month)){\n plants_df[which(plants_df$month == i), 'month'] <- dates[dates$month == i, 'inx']\n}\n\n# Generating a dummy set for plants\nnames(plants_df)[names(plants_df) == 'Plant.Names'] <- ''\nplants_df <- cbind.data.frame(plants_df, model.matrix( ~ Plant.Name -1, data = plants_df))\nnames(plants_df)[5:ncol(plants_df)] <- substr(names(plants_df)[5:ncol(plants_df)],\n nchar('Plant.Name') + 1,\n nchar(names(plants_df)[5:ncol(plants_df)]))\n\n\n# Adding locations WITHOUT electricity stations\ngrid_for_no_plants <- expand.grid('ID' = c(dem_df[!(dem_df$ID %in% plants_df$ID), 'ID']),\n 'month' = unique(plants_df$month), stringsAsFactors = FALSE)\n\n# Merging locations with electricity plants and without\nplants_df_ <- rbind.data.frame(plants_df, cbind.data.frame(grid_for_no_plants,\n 'estimated_generation_mwh' = 0,\n 'Plant.Name' = 'NoPlant',\n 'AESPR__Jobos' = 0, 'Aguirre' = 0,\n 'Cambalache' = 0,\n 'Central_SJ__Palo_Seco' = 0,\n 'Daguao' = 0,\n 'EcoEl__Costa_Sur' = 0,\n 'HIMA_San_Pablo_Caguas' = 0,\n 'Janssen_Ortho_LLC' = 0, 'Mayaguez' = 0, \n 'Vega_Baja' = 0, \n 'Yabucoa' = 0,\n 'Culebra' = 0, 'Vieques' = 0))\n\n# Convert ID to character\nhighways_df$ID <- as.character(highways_df$ID)\n\n# Merging all dataframes\ndf <- no2_df_long %>%\n left_join(dem_df, by = 'ID') %>%\n left_join(highways_df, by = 'ID') %>%\n left_join(temperature_df_long, by = c('ID', 'time')) %>%\n left_join(humidity_df_long, by = c('ID', 'time')) %>%\n left_join(u_wind_df_long, by = c('ID', 'time')) %>%\n left_join(v_wind_df_long, by = c('ID', 'time')) %>%\n left_join(precipitable_water_df_long, by = c('ID', 'time')) %>%\n left_join(pressure_df_long, by = c('ID', 'time')) %>%\n left_join(cloud_fraction_df_long, by = c('ID', 'time'))\n\n\n# Monthly and weekly dummies for further interactions\n\ndf$week <- substr(df$time, 6, nchar(df$time))\ndf$year <- substr(df$time, 1, 4)\ndf$date <- as.Date(paste(df$year, df$week, 1, sep=\"-\"), \"%Y-%U-%u\")\ndf$month <- month(ymd(df$date))\n\n# Format fixing\nplants_df_$month <- as.numeric(plants_df_$month)\nplants_df_$ID <- as.character(plants_df_$ID)\nplants_df_$Plant.Name <- as.character(plants_df_$Plant.Name)\n\n# Merging df with the data on plants\ndf <- df %>%\n left_join(plants_df_, by = c('ID', 'month'))\n\n# Dummy for the presence of an electricity plant\ndf$plant <- ifelse(df$estimated_generation_mwh == 0,0,1)\n# table(df$plant)\n\n# Dropping missing values and hence obtainig a BALANCED PANEL\ntime_to_drop <- df[is.na(df$no2), 'time']\ndf <- df[!df$time %in% time_to_drop, ]\n\n# Arranging\ndf <- df[, c('ID', 'time','no2','estimated_generation_mwh', 'Plant.Name',\n 'population', 'temp', 'humid', 'uwind', 'vwind', 'precip',\n 'pressure', 'cloud',\n 'plant', 'highway',\n 'month', \n 'AESPR__Jobos', 'Aguirre',\n 'Cambalache',\n 'Central_SJ__Palo_Seco',\n 'Daguao',\n 'EcoEl__Costa_Sur',\n 'HIMA_San_Pablo_Caguas',\n 'Janssen_Ortho_LLC', 'Mayaguez', \n 'Vega_Baja', \n 'Yabucoa', 'Culebra', 'Vieques')]\n\n# Genrating overall WEEKLY mwh\ndf$mwh_weekly <- df$estimated_generation_mwh/4\n\n# Genrating weekly mwh fro each plant\ndf$AESPR__Jobos_mwh <- df$mwh_weekly*df$AESPR__Jobos\ndf$Aguirre_mwh <- df$mwh_weekly*df$Aguirre\ndf$Cambalache_mwh <- df$mwh_weekly*df$Cambalache\ndf$Central_SJ__Palo_Seco_mwh <- df$mwh_weekly*df$Central_SJ__Palo_Seco\ndf$Culebra_mwh <- df$mwh_weekly*df$Culebra\ndf$Daguao_mwh <- df$mwh_weekly*df$Daguao\ndf$EcoEl__Costa_Sur_mwh <- df$mwh_weekly*df$EcoEl__Costa_Sur\ndf$HIMA_San_Pablo_Caguas_mwh <- df$mwh_weekly*df$HIMA_San_Pablo_Caguas\ndf$Janssen_Ortho_LLC_mwh <- df$mwh_weekly*df$Janssen_Ortho_LLC\ndf$Mayaguez_mwh <- df$mwh_weekly*df$Mayaguez\ndf$Vega_Baja_mwh <- df$mwh_weekly*df$Vega_Baja\ndf$Vieques_mwh <- df$mwh_weekly*df$Vieques\ndf$Yabucoa_mwh <- df$mwh_weekly*df$Yabucoa\n\n# Wind velocity\ndf$wind_velocity <- sqrt(df$uwind^2 + df$vwind^2)\n\n# NO2 in tons\nm2_per_cell <- 95*1e6\ndf$no2_kg <- df$no2*m2_per_cell*46.00053/1000\n\n# Month as character\ndf$month <- as.factor(df$month)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Testing whether we need an autoregressive error or not\nfm0 <- no2_kg ~ wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\nbsktest_clm <- bsktest(x = fm0, data = df, listw = spatial_mat_list,\n test = \"CLMlambda\")\nbsktest_clm # the test is dignificant -> the error component is necessary","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Findings\n\n### Computational Nuances\n\nTo compute a spatial weight matrix, we used a k-nearest neighbors function, i.e., a function that returns a matrix with point indices, belonging to the set of k nearest neighbors for each available observation. The dataset with all the variables we concentrated on was converted into a long format with the first variable, representing IDs of the spatial units within the considered geographical area (Puerto Rico), and the second variable, responsible for the periods (we selected weeks as time units).\n\n\nThe coordinates that were plugged into the KNN function depicted centers of the respective rectangular raster cells. The number of cells within the whole sub-national region was manually chosen to be equal to 192. Such a decision reflected our attempt to ensure an optimal match between the ability of a cell to 'separate' the power plants as analytical units (without adding much of noise to it) and to grasp a reasonable portion of the territory. Due to the relative closeness between some electricity stations, we did not manage to split them without disrupting the symmetry of the grid. In such cases, 'inseparable' plants were processed together and reflected a single observation. See below a visualized example of how a spatial unit of analysis and its neighbors look like.\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Create data for spatial lag explanation\nspatial_lag_example_raster <- no2_rasters[[1]]\nvalues_spatial_lag_example <- c(rep(NA, 192))\nvalues_spatial_lag_example[57] <- 0.5\nvalues_spatial_lag_example[58] <- 0.5\nvalues_spatial_lag_example[59] <- 0.5\nvalues_spatial_lag_example[81] <- 0.5\nvalues_spatial_lag_example[83] <- 0.5\nvalues_spatial_lag_example[105] <- 0.5\nvalues_spatial_lag_example[106] <- 0.5\nvalues_spatial_lag_example[107] <- 0.5\nvalues_spatial_lag_example[82] <- 1\n\nvalues(spatial_lag_example_raster) <- values_spatial_lag_example\nspatial_lag_example_data <- gplot_data(spatial_lag_example_raster, maxpixels = 10000) # Choose number of pixels\n\n# Create Puerto Rico bbox\nbbox_pr <- make_bbox(lon=c(-67.32431, -65.19082), lat=c(17.84566, 18.5694))\n# Download map of Puerto Rico\nmap_pr <- get_stamenmap(bbox_pr,\n maptype=\"terrain-background\")\n\n### Visualize spatial lag example\nggmap(map_pr) +\n geom_tile(data = spatial_lag_example_data, aes(x, y, fill = value), show.legend = FALSE) +\n scale_fill_gradient(low =adjustcolor(\"lightblue\", alpha.f = 0.7),\n high = adjustcolor(\"blue\", alpha.f = 0.6), na.value=adjustcolor(\"white\", alpha.f = 0.3)) + \n annotate(\"text\", label = \"A spatial unit of interest\", x = -66.12, y = 18.26, size = 9, colour = \"black\") + \n annotate(\"text\", label = \"Nearest neighbors of the spatial unit\", x = -66.75, y = 18.17, size = 7, colour = \"black\") + \n ggtitle(\"Example of the Nearest Neighbors of a Spatial Unit\") + theme(plot.title = element_text(size = 32, face = \"bold\", hjust = 0.5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A simple but critical for our analysis aspect that we account for refers to the fact that wind velocity determines the precision of registering the effect of interest: the highest possible precision can be achieved when wind velocity is equal to zero, meaning that NO2 does not get spread around fast. Contrarily, when the wind is strong, the amount of captured NO2 may be underestimated. Having said this, we then assume that, all things being equal, in periods when the wind velocity is not null, the impact of the generated electricity on NO2 emission is on average the same as in periods with zero wind, although technically we are not able to distinguish it. Nevertheless, we can extrapolate, i.e., contend that the revealed estimate for the windless period is equal to the one in the windy period because the wind does not directly affect the scope of pollution brought on by stations, it just prevents from accurately capturing it.\n\n**Therefore, the Average Emission Factor that we report in this work is described by interaction terms in Spatial Panel Models of the wind speed with the corresponding variables, measuring electricity generation by each plant. The overall AEF for the sub-national region of Puerto Rico is then computed as an average of the regression coefficients of the mentioned interaction terms that appeared to be statistically significant.** The interaction (or moderation) technique is usually employed to reveal the association between a dependent variable and an outcome of interest when this association depends on another variable.\n\nThis logic (imposed by the specificity of the data) implies we can compute the AEF only for a fraction of power stations that emit NO2 to a relatively pronounced extent, i.e., significantly different from the average pollution level in _spatial units lacking power plants_. If the overall pollution is so high everywhere and plants' activity does not stand out, we are not able to disentangle a portion of emissions attributable exclusively to these electricity stations. The latter scenario is only possible (but perhaps not guaranteed) in case of availability of all (accurately measured) spatially detachable factors related to the NO2 emission that we can control for. Regarding the case of Puerto Rico, we account for the population density and motorways/highways.\n\n### Obtaining the AEF\n\nThe table below demonstrates the estimation results of our base spatial panel model with the _**electricity production**_ for each station as predictors, _**wind velocity**_ as a moderator, _**population density**_ and _**motorways/highways**_ as control parameters. Red color highlights statistically significant on 95% confidence level effects for the two spatial units, containing 3 stations in total (given that wind is zero): **Central San Juan and Palo Seco,** and **Vega Baja**. Yellow color accentuates two significant interaction terms for **Vega Baja** and **Central San Juan and Palo Seco**. These terms describe _**the magnitude of reduction in our ability to capture NO2 emission**_ for a particular plant due to the wind. This is logical: NO2 is carried away by the wind and hence with the intensification of wind we end up fixating a lower amount of NO2."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"############## Base model\n# Running random effects with lag and with error\nspgm_base <- spgm(formula = fm0, data = df, listw = spatial_mat_list, lag = TRUE,\n moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n# summary(spgm_base)\n\n# Printing the output of the base model\nspgm_base_table <- round(as.data.frame(summary(spgm_base)[[\"CoefTable\"]]),4)\ncolnames(spgm_base_table) <- c(\"Estimate\", \"SE\", \"t-value\", \"p-value\")\n\n# ! You can get the results bellow by running the following script in RStudio\n# kable(spgm_base_table, \"html\") %>%\n# kable_styling(c(\"striped\"), full_width = F) %>%\n# row_spec(c(7,12), bold = T, background = \"lightcoral\")%>%\n# row_spec(c(20, 22, 25), bold = T,\n# background = \"khaki\")\n\n# Simplier version of the coefficient \n# kable(spgm_base_table, \"rst\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](\"https://sun9-20.userapi.com/c857636/v857636124/1aaac2/vNsyEzUQFok.jpg\")
"},{"metadata":{},"cell_type":"markdown","source":"Overall, the detectable AEF for Puerto Rico is equal to the mean of the significant plant-wise effects of electricity generation (megawatt-hours) on the amount of NO2 (kilograms), controlling for the population and motorways."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"AEF <- round(mean(c(spgm_base_table$Estimate[7],\n spgm_base_table$Estimate[12])),2)\ncat('The AEF for Puerto Rico:', AEF, 'kg/mwh (general model)')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Notice that the coefficient for the population variable is significant and positive, while for motorways/highways -- insignificant. This indicates that growth in population (other things being equal) leads to greater NO2, which seems logical. In turn, an insignificant effect for highways tells us that the intercity traffic in Puerto Rico is not heavy, and therefore areas _without highways_ on average are not different in terms of pollution from those _with highways_.\n\n### Factors Affecting the AEF\n\nFurther, we test if some weather characteristics can modify the amount of NO2 released due to energy generation by the power plants. We do that by introducing additional interaction terms into our base model, i.e., a model where wind velocity decomposes the overall impact of electricity on NO2 concentration. Each weather-related variable (except for wind velocity and cloudiness) was standardized with Z-scores and plugged into separate models for the sake of inferring meaningful and easily interpretable findings (_1 st.dev._ means a feature is highly manifested, _0 st. dev._ means a feature is moderately manifested, _-1 st. dev._ points on the low manifestation of a feature).\n\nThe graphs below display the magnitude of the interaction effects with confidence limits for the power plants that showed some statistically significant estimates. The significance test was performed using a general linear hypotheses function (glht).\n\nOverall, we were able to detect significant findings only for the temperature and cloudiness; moderation effects for such parameters as precipitation, humidity, and pressure were insignificant."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"############## Testing for factors that can potentially change AEF:\n# interactions with weather characteristics\n\n# Temperature formula\nfm1 <- no2_kg ~ scale(temp)*wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\n# Temperature model\nspgm_temp <- spgm(formula = fm1, data = df, listw = spatial_mat_list, lag = TRUE, \n moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n#summary(spgm_temp)\n\n# EcoEl__Costa_Sur_temp\n# linear hypotheses test for each level of temp: -1, 0, 1\nEcoEl__Costa_Sur_temp_1 <- confint(summary(glht(spgm_temp, linfct = \n c(\"EcoEl__Costa_Sur_mwh - scale(temp):EcoEl__Costa_Sur_mwh = 0\"))), level = 0.9)\n\nEcoEl__Costa_Sur_temp0 <- confint(summary(glht(spgm_temp, linfct = \n c(\"EcoEl__Costa_Sur_mwh = 0\"))), level = 0.9)\n\nEcoEl__Costa_Sur_temp1 <- confint(summary(glht(spgm_temp, linfct = \n c(\"EcoEl__Costa_Sur_mwh + scale(temp):EcoEl__Costa_Sur_mwh = 0\"))), level = 0.9)\n# Extracting coefs and confidence limits\nEcoEl__Costa_Sur_temp <- rbind.data.frame(EcoEl__Costa_Sur_temp_1$confint,\n EcoEl__Costa_Sur_temp0$confint,\n EcoEl__Costa_Sur_temp1$confint)\n# A column with names\nrownames(EcoEl__Costa_Sur_temp) <- NULL\nEcoEl__Costa_Sur_temp$names <- factor(c('low (-1 SD)', \n 'middle (0 SD)',\n 'high (1 SD)'),\n levels = c('low (-1 SD)', \n 'middle (0 SD)',\n 'high (1 SD)'))\n\n###############################################################################\n###############################################################################\n###############################################################################\n\n# Humidity --> none of the interactions for zero wind appeared to be significant\nfm2 <- no2_kg ~ scale(humid)*wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\n\n#spgm_humid <- spgm(formula = fm2, data = df, listw = spatial_mat_list, lag = TRUE, \n# moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n#summary(spgm_humid)\n\n###############################################################################\n###############################################################################\n###############################################################################\n\n# Precipitation --> none of the interactions for zero wind appeared to be significant\nfm3 <- no2_kg ~ scale(precip)*wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\n\n\n#spgm_precip <- spgm(formula = fm3, data = df, listw = spatial_mat_list, lag = TRUE, \n# moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n#summary(spgm_precip)\n\n###############################################################################\n###############################################################################\n###############################################################################\n\n# Pressure --> none of the interactions for zero wind appeared to be significant\nfm4 <- no2_kg ~ scale(pressure)*wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\n\n#spgm_pressure <- spgm(formula = fm4, data = df, listw = spatial_mat_list, lag = TRUE, \n# moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n#summary(spgm_pressure)\n\n###############################################################################\n###############################################################################\n###############################################################################\n\n# Clouds\nfm6 <- no2_kg ~ cloud*wind_velocity*(AESPR__Jobos_mwh +\n Aguirre_mwh +\n Cambalache_mwh + \n Central_SJ__Palo_Seco_mwh +\n Daguao_mwh +\n EcoEl__Costa_Sur_mwh +\n Janssen_Ortho_LLC_mwh + \n Mayaguez_mwh +\n Vega_Baja_mwh +\n Yabucoa_mwh +\n Culebra_mwh) + population + highway\n\nspgm_cloud <- spgm(formula = fm6, data = df, listw = spatial_mat_list, lag = TRUE, \n moments = \"fullweights\", model = \"random\", spatial.error = TRUE)\n#summary(spgm_cloud)\n\n################# Cambalache_cloud\n# linear hypotheses test for several levels of cloudiness\n\nCambalache_cloud_0 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Cambalache_mwh = 0\"))), level = 0.9)\n\nCambalache_cloud_0.25 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Cambalache_mwh + 0.25*cloud:Cambalache_mwh = 0\"))), level = 0.9)\n\nCambalache_cloud_0.5 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Cambalache_mwh + 0.5*cloud:Cambalache_mwh = 0\"))), level = 0.9)\n\nCambalache_cloud_0.75 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Cambalache_mwh + 0.75*cloud:Cambalache_mwh = 0\"))), level = 0.9)\n\nCambalache_cloud_1 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Cambalache_mwh + 1*cloud:Cambalache_mwh = 0\"))), level = 0.9)\n# Extracting coefs and confidence limits\nCambalache_cloud <- rbind.data.frame(Cambalache_cloud_0$confint,\n Cambalache_cloud_0.25$confint,\n Cambalache_cloud_0.5$confint,\n Cambalache_cloud_0.75$confint,\n Cambalache_cloud_1$confint)\n# A column with names\nrownames(Cambalache_cloud) <- NULL\nCambalache_cloud$cloudiness <- c(0, 0.25, 0.5, 0.75, 1)\n\n################# Central_SJ__Palo_Seco_cloud\n\nCentral_SJ__Palo_Seco_cloud_0 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Central_SJ__Palo_Seco_mwh = 0\"))), level = 0.9)\n\nCentral_SJ__Palo_Seco_cloud_0.25 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Central_SJ__Palo_Seco_mwh + 0.25*cloud:Central_SJ__Palo_Seco_mwh = 0\"))), level = 0.9)\n\nCentral_SJ__Palo_Seco_cloud_0.5 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Central_SJ__Palo_Seco_mwh + 0.5*cloud:Central_SJ__Palo_Seco_mwh = 0\"))), level = 0.9)\n\nCentral_SJ__Palo_Seco_cloud_0.75 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Central_SJ__Palo_Seco_mwh + 0.75*cloud:Central_SJ__Palo_Seco_mwh = 0\"))), level = 0.9)\n\nCentral_SJ__Palo_Seco_cloud_1 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Central_SJ__Palo_Seco_mwh + 1*cloud:Central_SJ__Palo_Seco_mwh = 0\"))), level = 0.9)\n# Extracting coefs and confidence limits\nCentral_SJ__Palo_Seco_cloud <- rbind.data.frame(Central_SJ__Palo_Seco_cloud_0$confint,\n Central_SJ__Palo_Seco_cloud_0.25$confint,\n Central_SJ__Palo_Seco_cloud_0.5$confint,\n Central_SJ__Palo_Seco_cloud_0.75$confint,\n Central_SJ__Palo_Seco_cloud_1$confint)\n# A column with names\nrownames(Central_SJ__Palo_Seco_cloud) <- NULL\nCentral_SJ__Palo_Seco_cloud$cloudiness <- c(0, 0.25, 0.5, 0.75, 1)\n\n################# Yabucoa_cloud\n# linear hypotheses test for several points of cloudindess\nYabucoa_cloud_0 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Yabucoa_mwh = 0\"))), level = 0.9)\n\nYabucoa_cloud_0.25 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Yabucoa_mwh + 0.25*cloud:Yabucoa_mwh = 0\"))), level = 0.9)\n\nYabucoa_cloud_0.5 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Yabucoa_mwh + 0.5*cloud:Yabucoa_mwh = 0\"))), level = 0.9)\n\nYabucoa_cloud_0.75 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Yabucoa_mwh + 0.75*cloud:Yabucoa_mwh = 0\"))), level = 0.9)\n\nYabucoa_cloud_1 <- confint(summary(glht(spgm_cloud, linfct = \n c(\"Yabucoa_mwh + 1*cloud:Yabucoa_mwh = 0\"))), level = 0.9)\n# Extracting coefs and confidence limits\nYabucoa_cloud <- rbind.data.frame(Yabucoa_cloud_0$confint,\n Yabucoa_cloud_0.25$confint,\n Yabucoa_cloud_0.5$confint,\n Yabucoa_cloud_0.75$confin,\n Yabucoa_cloud_1$confint)\n# A column with names\nrownames(Yabucoa_cloud) <- NULL\nYabucoa_cloud$cloudiness <- c(0, 0.25, 0.5, 0.75, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Plotting\na <- ggplot(EcoEl__Costa_Sur_temp, aes(x = Estimate, y = names)) +\n geom_point(size = 4) +\n geom_text(aes(label = round(Estimate, 3), vjust = -1)) +\n geom_errorbar(aes(xmax = upr, xmin = lwr), width=0.3) +\n geom_vline(xintercept = 0, linetype=\"dashed\", \n color = \"red\", size=1) +\n xlab(\"Plant-specific AEF and its CI\") +\n ylab(\"Temperature\") +\n ggtitle(\"EcoElectrica/Costa Sur Power Plants\") +\n theme(axis.line = element_line(colour = \"black\"),\n panel.grid = element_blank(),\n plot.title = element_text(hjust = 0.5, size=18),\n axis.title = element_text(size = 16),\n axis.text = element_text(color = 'black', size = 12))\n\n# Plotting\nb <- ggplot(Central_SJ__Palo_Seco_cloud, aes(x = Estimate, y = cloudiness)) +\n geom_line() +\n geom_point(size = 4) +\n geom_text(aes(label = round(Estimate, 3), vjust = -1)) +\n geom_errorbar(aes(xmax = upr, xmin = lwr), width = 0.05) +\n geom_vline(xintercept = 0, linetype=\"dashed\", \n color = \"red\", size=1) +\n xlab(\"Plant-specific AEF and its CI\") +\n ylab(\"Cloudiness\") +\n ggtitle(\"Central San Juan and Palo Seco Power Plant\") +\n theme(axis.line = element_line(colour = \"black\"),\n panel.grid = element_blank(),\n plot.title = element_text(hjust = 0.5, size=18),\n axis.title = element_text(size = 16),\n axis.text = element_text(color = 'black', size = 12)) +\n scale_y_continuous(limits = c(-0.1, 1.2))\n\n# Plotting\nc <- ggplot(Cambalache_cloud, aes(x = Estimate, y = cloudiness)) +\n geom_line() +\n geom_point(size = 4) +\n geom_text(aes(label = round(Estimate, 3), vjust = -1)) +\n geom_errorbar(aes(xmax = upr, xmin = lwr), width = 0.05) +\n geom_vline(xintercept = 0, linetype=\"dashed\", \n color = \"red\", size=1) +\n xlab(\"Plant-specific AEF and its CI\") +\n ylab(\"Cloudiness\") +\n ggtitle(\"Cambalache Power Plant\") +\n theme(axis.line = element_line(colour = \"black\"),\n panel.grid = element_blank(),\n plot.title = element_text(hjust = 0.5, size=18),\n axis.title = element_text(size = 16),\n axis.text = element_text(color = 'black', size = 12)) +\n scale_y_continuous(limits = c(-0.1, 1.2))\n\n# Plotting\nd <- ggplot(Yabucoa_cloud, aes(x = Estimate, y = cloudiness)) +\n geom_line() +\n geom_point(size = 4) +\n geom_text(aes(label = round(Estimate, 3), vjust = -1)) +\n geom_errorbar(aes(xmax = upr, xmin = lwr), width = 0.05) +\n geom_vline(xintercept = 0, linetype=\"dashed\", \n color = \"red\", size=1) +\n xlab(\"Plant-specific AEF and its CI\") +\n ylab(\"Cloudiness\") +\n ggtitle(\"Yabucoa Power Plant\") +\n theme(axis.line = element_line(colour = \"black\"),\n panel.grid = element_blank(),\n plot.title = element_text(hjust = 0.5, size=18),\n axis.title = element_text(size = 16),\n axis.text = element_text(color = 'black', size = 12)) +\n scale_y_continuous(limits = c(-0.1, 1.2))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first graph depicts a significant (on 95% confidence interval) moderation effects with the temperature for the **EcoElectrica** and **Costa Sur** electricity stations. When the temperature is high these plants demonstrate a slight increase in the AEF - under high temperatures, the plants increase NO2 emission. However, all negative coefficients indicate that on average growth in the temperature is not enough for these plants to become more pollutant than _an average spatial unit without plants_: the plants seem to be really ecologically friendly (as it actually stems from the name)."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Grouping plots\noptions(repr.plot.width=8, repr.plot.height = 6)\nfig_a <- ggarrange(a, ncol = 1, nrow = 1)\nannotate_figure(fig_a, top = \n text_grob(\"Interaction Effects for Electricity Generation (mwh) and Temperature\", size = 18)) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The next panel displays significant moderation effects for the cloudiness: lower cloudiness increases the AEFs of the **Central San Juan, Palo Seco, Cambalache** and **Yabucoa** power plants so that AEFs become significant for almost all of them. This seems reasonable since cloudiness impedes the collection of precise data through remote sensing techniques, whereas the absence of cloudiness, perhaps, allows capturing more NO2 emissions."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=16, repr.plot.height = 6)\nfig_bcd <- ggarrange(b, c, d, ncol = 3, nrow = 1)\nannotate_figure(fig_bcd, top = \n text_grob(\"Interaction Effects for Electricity Generation (mwh) and Cloudiness\", size = 18)) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Monthly AEF\n\nOur next step is to focus on smaller time slices in the computation of the AEF. We took a month as an example and introduced additional interaction terms (separately for each electricity station) in our base model with _month as a factor moderator_. The significance of the obtained interaction terms was checked as previously by the function glht (General Linear Hypotheses). We estimated the plant-wise effects of the electricity production on NO2 and averaged (across plants and then across months) only those coefficients that happened to be significant.\n\nIt was revealed that if we concentrate on smaller periods, there is actually a greater number of plants that give off a substantial amount of NO2: in addition to **Central San Juan and Palo Seco,** and **Vega Baja**, the **Cambalache** station demonstrated importance in July, August, October, and December. In total, 7 months out of 12 were shown to be significant in our model; the estimates for the rest were not established by our model with certainty (i.e., could be equal to zero). Such a detailed approach seems more advantageous because in the general method shown above the plant-wise estimates are averaged across all months, including those with potentially insignificant parameters; contrariwise when focusing on smaller time slices, we extract and process only the substantial (significant) NO2 contributions by months."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# A vector with plats under focus\n# We drope EcoEl__Costa_Sur_mwh since for this plant a model with an interaction\n# was not able to converge\nmwh <- c('AESPR__Jobos_mwh', 'Cambalache_mwh', \"Central_SJ__Palo_Seco_mwh\",\n 'Daguao_mwh', 'EcoEl__Costa_Sur_mwh', 'Janssen_Ortho_LLC_mwh', 'Mayaguez_mwh',\n 'Vega_Baja_mwh', 'Yabucoa_mwh', 'Culebra_mwh')\n\n# Grid for interaction terms\nmonth_mwh <- expand.grid('month', mwh, stringsAsFactors = FALSE)\nmonth_mwh_interact <- paste(do.call(paste, c(month_mwh, sep=\"*\")))\n\n# Creating formulas for an interction of month with each plants' mwh generation in separate models\nfm7 <- c()\nfor (i in 1:length(mwh)){\n fm7[i] <- paste0(paste0(\"fm7_\", i, \" <- no2_kg ~ wind_velocity*(\", month_mwh_interact[i], sep = \" + \"),\n paste0(mwh[-i], collapse = \" + \"), ') + population + highway')\n}\n\n# Plugging these formulas into the GM estimation of spatial panel data models\nspgm_month <- vector(\"list\", length(mwh))\nfor (i in 1:length(mwh)){\n spgm_month[[i]] <- spgm(formula = formula(fm7[i]),\n data = df,\n listw = spatial_mat_list, lag = TRUE, \n moments = \"fullweights\",\n model = \"random\", spatial.error = TRUE)\n # print(mwh[i])\n}\n\n# Then we have to test for significance of the obtained interactions of months with\n# the plants mwh generation given wind velocity is 0 (that is our assumption that allows\n# to capture the effect that we are interested in).\n\n# We do that with the fuction glht (General Linear Hypotheses) in a loop and plot the results\n\n# Creating a vector - coefs - in which we write the results of glht\nn_months <- 12\nseq <- 2:n_months\ncoefs <- vector(\"list\", length(mwh))\ncoefs2 <- vector(\"list\", n_months)\n\n# List inside a list\nfor (i in 1:length(coefs)){\n coefs[[i]] <- coefs2\n}\n\n# Creating a loop for glht and extracting confidence intervals of the respective point estimates\nfor (i in 1:length(mwh)){\n for (j in seq){\n coefs[[i]][[1]] = eval(parse(text = paste0('confint(summary(glht(spgm_month[[', i, ']], linfct = c(',\n paste0('\"', mwh[i], '= 0\"))), level = 0.9)' ), collapse = \"\")))\n \n coefs[[i]][[j]] = eval(parse(text = paste0('confint(summary(glht(spgm_month[[', i, ']], linfct = c(',\n paste0('\"', mwh[i], ' + month', j, ':', mwh[i], '= 0\"))), level = 0.9)' ), collapse = \"\")))\n }\n}\n\n# Transforming elements in the list into dataframes that represent estimation for\n# a plant within 12 years --> 10 such dataframes\n\ncoefs_fin <- vector(\"list\", length(mwh))\nfor (n in 1:length(coefs_fin)){\n coefs_fin[[n]] <- do.call(rbind.data.frame, lapply(1:n_months, function(i) as.numeric(coefs[[n]][[i]]$confint)))\n}\n\n# Naming the columns\nfor (i in 1:length(mwh)){\n names(coefs_fin[[i]]) <- c(\"est\", \"low\", 'up')\n}\n\n# We need to extract monthly coefficients that are significant for each plant\n\n# Setting all insignificant estmates to zero, leaving intact the significant ones\ncoefs_fin_sign <- coefs_fin\nfor (i in 1:length(coefs_fin)){\n for (j in 1:n_months){\n coefs_fin_sign[[i]]$est[j] <- ifelse(coefs_fin[[i]]$low[j] > 0 & coefs_fin[[i]]$up[j] > 0,\n coefs_fin[[i]]$est[j], 0)\n }\n}\n\n# Keeping only estimates (dropping CIs)\nmonthly_AEF <- list()\nfor (i in 1:length(mwh)){\n monthly_AEF[[i]] <- coefs_fin_sign[[i]]$est\n}\n\n# Transforming our list with estimates to a dataframe\nmonthly_AEF_df <- as.data.frame(monthly_AEF)\ncolnames(monthly_AEF_df) <- mwh\n\n# Computing monthly AEF\nmonthly_AEF_df <- monthly_AEF_df %>% replace(. == 0, NA)\nmonthly_AEF_df$monthly_AEF <- rowMeans(monthly_AEF_df, na.rm = T)\nmonthly_AEF_df$monthly_AEF %<>% replace(is.nan(monthly_AEF_df$monthly_AEF), 0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"AEF_monthly <- round(mean(monthly_AEF_df[is.na(monthly_AEF_df$monthly_AEF) == F, 'monthly_AEF']), 2)\ncat('The AEF for Puerto Rico:', AEF_monthly, 'kg/mwh (monthly model)')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The following picture depicts the AEFs for Puerto Rico by months: brighter colors depict greater AEFs; white space is shown for months with statistically insignificant estimates."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#### Plotting Monthly AEF\noptions(repr.plot.width=16, repr.plot.height = 6)\noptions(warn=-1)\n\n# A vector with month names\nmonth_names<- c('January', 'February', 'March', 'April', \"May\", 'June',\n 'July', 'August', 'September', 'October', 'November', 'December')\n\nmonth_names_levels <- c('July', 'August', 'September', 'October', 'November', 'December',\n 'January', 'February', 'March', 'April', \"May\", 'June')\n\n# Arranging month names\nmonth_names_f <- factor(month_names,\n levels=as.character(month_names_levels))\n\nmonthly_AEF_df <- monthly_AEF_df %>% replace(. == 0, NA)\n\n# Plotting\nggplot(data = monthly_AEF_df, aes(x = month_names_f, y = NA)) +\ngeom_tile(aes(fill = monthly_AEF), colour = \"grey\")+\nscale_fill_distiller(na.value = \"white\", palette = 'Reds', direction = 1, name=\"AEF, kg/MWh\") +\ngeom_text(aes(label = round(monthly_AEF, 3), vjust = -1), size=8) +\ntheme(panel.grid = element_blank(),\npanel.background = element_blank(),\nplot.title = element_text(hjust = 0.5, size=18),\naxis.title = element_text(size = 16),\naxis.title.x = element_blank(),\naxis.text.y = element_blank(),\naxis.text.x = element_text(angle = 30, color = 'black',size = 16),\nplot.margin = unit(c(3,1,3,1), \"cm\")) +\nggtitle(\"Monthly registered AEF in Puerto Rico, 2018 - 2019\") +\nylab(\"AEF, kg/MWh\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Obtaining the MEF\nWe also propose a way to compute the MEF using the same methodology. Following Siler-Evans, Azevedo and Morga (2012), we use *differenced values* of NO2 (kg) and plant-specific energy generation (MWh) to obtain the MEF for Puerto Rico. Differencing in time series refers to the creation of an absolute difference between the current value of a variable and its value in the previous period. Since the energy production variable was available to us only by months, for this task we leveraged the data aggregated by months."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Reading the data\n# no2_df <- readRDS('input/saved-data-month/data_rds_192_month/no2_df.rds')\n# cells_coordinates_df <- readRDS('input/saved-data-month/data_rds_192_month/cells_coordinates_df.rds')\n# plants_monthly_df <- readRDS('input/saved-data-month/data_rds_192_month/plants_monthly_df.rds')\n# u_wind_df <- readRDS('input/saved-data-month/data_rds_192_month/u_wind_df.rds')\n# v_wind_df <- readRDS('input/saved-data-month/data_rds_192_month/v_wind_df.rds')\n# highways_df <- readRDS('input/saved-data-month/data_rds_192_month/highways_df.rds')\n# dem_df <- readRDS('input/saved-data-month/data_rds_192_month/dem_df.rds')","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"# Create dataframes for no2, pressure and cloud fraction\nno2_df <- createDataframeFromRasters_SP5NO2(no2_rasters, monthly=T, no2=T)\npressure_df <- createDataframeFromRasters_SP5NO2(pressure_rasters, monthly=T)\ncloud_fraction_df <- createDataframeFromRasters_SP5NO2(cloud_fraction_rasters, monthly=T)\n\n# Create dataframes for weather conditions\ntemperature_df <- createDataframeFromRasters_GFS(temperature_rasters, no2_rasters, monthly=T)\nhumidity_df <- createDataframeFromRasters_GFS(humidity_rasters, no2_rasters, monthly=T)\nu_wind_df <- createDataframeFromRasters_GFS(u_wind_rasters, no2_rasters, monthly=T)\nv_wind_df <- createDataframeFromRasters_GFS(v_wind_rasters, no2_rasters, monthly=T)\nprecipitable_water_df <- createDataframeFromRasters_GFS(precipitable_water_rasters, no2_rasters, monthly=T)\n\n\n##### Create a spatial matrix for observations\n# Spatial weights matrix\ncoords <- coordinates(cells_coordinates_df[, 1:2])\n# KNN\nnb <- knn2nb(knearneigh(coords, k = 8), row.names = cells_coordinates_df[,3])\nspatial_mat_list <- nb2listw(nb)\n\n# Removing outliers for NO2\nno2_df[no2_df > quantile(as.matrix(no2_df), 0.99)] <- quantile(as.matrix(no2_df), 0.99)\n\n##### Convert to long format\n\n# Define the number of spatial units\nn_spatial_units <- ncol(u_wind_df)\n\n# Function for convertion to long format\nShortToLong <- function(df, n, varname, Xlead = FALSE){\n df <- eval(parse(text = paste(\"cbind.data.frame(gather(as.data.frame(df), ID,\",\n varname, \"), 'time' = rep(rownames(df), n))\")))\n df$time <- as.character(df$time)\n if (Xlead == TRUE){\n # if leading X appears in a time variable\n df$time <- as.Date(substr(df$time, 2, nchar(df$time)), format=\"%Y%m%d\")\n }\n else{\n df$time <- df$time\n }\n df$time <- as.character(df$time)\n return(df)\n}\n\n# Applying to the dataframes\nu_wind_df_long <- ShortToLong(u_wind_df, n_spatial_units, 'uwind', Xlead = F)\nv_wind_df_long <- ShortToLong(v_wind_df, n_spatial_units, 'vwind', Xlead = F)\nno2_df_long <- ShortToLong(no2_df, n_spatial_units, 'no2')\n\n# Merge plants which were placed in one cell\nplants_monthly_df_ <- plants_monthly_df %>%\n group_by(ID) %>%\n summarise_at(c(\"netgen_July\", \"netgen_August\", \"netgen_September\", \"netgen_October\", \n \"netgen_November\", \"netgen_December\", \"netgen_January\", \"netgen_February\", \n \"netgen_March\", \"netgen_April\", \"netgen_May\", \"netgen_June\" ), sum)\n\n# Renaming 'joint' plants to understand that they actually represent several plants, not one\nplants_monthly_df[plants_monthly_df$ID == 'V158', 'Plant.Name'] <- 'AESPR__Jobos' \nplants_monthly_df[plants_monthly_df$ID == 'V151', 'Plant.Name'] <- 'EcoEl__Costa_Sur'\nplants_monthly_df[plants_monthly_df$ID == 'V38', 'Plant.Name'] <- 'Central_SJ__Palo_Seco'\n\n# Remove duplicates by ID\nplants_monthly_df_ <- plants_monthly_df_[!duplicated(plants_monthly_df_$ID), ]\n#saveRDS(plants_monthly_df_, 'plants_monthly_df_merged.rds')\n\n# Converting to long format the data on plants\nplants_df <- gather(plants_monthly_df_, ID, 'estimated_generation_mwh')\ncolnames(plants_df)[1] <- 'month'\nPlant.Name_ID <- plants_monthly_df[, c('Plant.Name', 'ID')]\nPlant.Name_ID <- Plant.Name_ID[!duplicated(Plant.Name_ID$ID),] %>% arrange(ID)\nplants_df <- cbind.data.frame(plants_df, 'ID' = rep(Plant.Name_ID$ID, 12),\n 'Plant.Name' = rep(Plant.Name_ID$Plant.Name, 12)) \n\n\n# Month as a number\ndates <- cbind.data.frame('month' = unique(plants_df$month), 'inx' = c(7:12, 1:6))\nfor (i in unique(plants_df$month)){\n plants_df[which(plants_df$month == i), 'month'] <- dates[dates$month == i, 'inx']\n}\n\n# Generating a dummy set for plants\nnames(plants_df)[names(plants_df) == 'Plant.Names'] <- ''\nplants_df <- cbind.data.frame(plants_df, model.matrix( ~ Plant.Name -1, data = plants_df))\nnames(plants_df)[5:ncol(plants_df)] <- substr(names(plants_df)[5:ncol(plants_df)],\n nchar('Plant.Name') + 1,\n nchar(names(plants_df)[5:ncol(plants_df)]))\n\n\n# Adding locations WITHOUT electricity stations\ngrid_for_no_plants <- expand.grid('ID' = c(dem_df[!(dem_df$ID %in% plants_df$ID), 'ID']),\n 'month' = unique(plants_df$month), stringsAsFactors = FALSE)\n\n# Merging locations with electricity plants and without\nplants_df_ <- rbind.data.frame(plants_df, cbind.data.frame(grid_for_no_plants,\n 'estimated_generation_mwh' = 0,\n 'Plant.Name' = 'NoPlant',\n 'AESPR__Jobos' = 0, 'Aguirre' = 0,\n 'Cambalache' = 0,\n 'Central_SJ__Palo_Seco' = 0,\n 'Daguao' = 0,\n 'EcoEl__Costa_Sur' = 0,\n 'HIMA_San_Pablo_Caguas' = 0,\n 'Janssen_Ortho_LLC' = 0, 'Mayaguez' = 0, \n 'Vega_Baja' = 0, \n 'Yabucoa' = 0,\n 'Culebra' = 0, 'Vieques' = 0))\n\n# Convert ID to character\nhighways_df$ID <- as.character(highways_df$ID)\n\n# Merging all dataframes\ndf <- no2_df_long %>%\n left_join(dem_df, by = 'ID') %>%\n left_join(highways_df, by = 'ID') %>%\n left_join(u_wind_df_long, by = c('ID', 'time')) %>%\n left_join(v_wind_df_long, by = c('ID', 'time'))\n \n\n# Months\ndf$month <- as.numeric(substr(df$time, 6, nchar(df$time)))\n\n# Format fixing\nplants_df_$month <- as.numeric(plants_df_$month)\nplants_df_$ID <- as.character(plants_df_$ID)\nplants_df_$Plant.Name <- as.character(plants_df_$Plant.Name)\n\n# Merging df with the data on plants\ndf <- df %>%\n left_join(plants_df_, by = c('ID', 'month'))\n\n# Dummy for the presence of an electricity plant\ndf$plant <- ifelse(df$estimated_generation_mwh == 0,0,1)\n# table(df$plant)\n\n# Dropping missing values and hence obtainig a BALANCED PANEL\ntime_to_drop <- df[is.na(df$no2), 'time']\ndf <- df[!df$time %in% time_to_drop, ]\n\n# Arranging\ndf <- df[, c('ID', 'time','no2',\n 'estimated_generation_mwh', 'Plant.Name',\n 'uwind', 'vwind', 'population', \n 'plant', 'highway', \n 'month', \n 'AESPR__Jobos', 'Aguirre',\n 'Cambalache',\n 'Central_SJ__Palo_Seco',\n 'Daguao',\n 'EcoEl__Costa_Sur',\n 'HIMA_San_Pablo_Caguas',\n 'Janssen_Ortho_LLC', 'Mayaguez', \n 'Vega_Baja', \n 'Yabucoa', 'Culebra', 'Vieques')]\n\n# renaming\ndf$mwh_monthly <- df$estimated_generation_mwh\n\n# Genrating weekly mwh fro each plant\ndf$AESPR__Jobos_mwh <- df$mwh_monthly*df$AESPR__Jobos\ndf$Aguirre_mwh <- df$mwh_monthly*df$Aguirre\ndf$Cambalache_mwh <- df$mwh_monthly*df$Cambalache\ndf$Central_SJ__Palo_Seco_mwh <- df$mwh_monthly*df$Central_SJ__Palo_Seco\ndf$Culebra_mwh <- df$mwh_monthly*df$Culebra\ndf$Daguao_mwh <- df$mwh_monthly*df$Daguao\ndf$EcoEl__Costa_Sur_mwh <- df$mwh_monthly*df$EcoEl__Costa_Sur\ndf$HIMA_San_Pablo_Caguas_mwh <- df$mwh_monthly*df$HIMA_San_Pablo_Caguas\ndf$Janssen_Ortho_LLC_mwh <- df$mwh_monthly*df$Janssen_Ortho_LLC\ndf$Mayaguez_mwh <- df$mwh_monthly*df$Mayaguez\ndf$Vega_Baja_mwh <- df$mwh_monthly*df$Vega_Baja\ndf$Vieques_mwh <- df$mwh_monthly*df$Vieques\ndf$Yabucoa_mwh <- df$mwh_monthly*df$Yabucoa\n\n# Wind velocity\ndf$wind_velocity <- sqrt(df$uwind^2 + df$vwind^2)\n\n# NO2 in tons\nm2_per_cell <- 95*1e6\ndf$no2_kg <- df$no2*m2_per_cell*46.00053/1000\n\n# Testing the correctness of merging\n# df[df$month == 6 & df$Plant.Name == \"Daguao\", 'estimated_generation_mwh']\n\n# Generating differenced data\nvars_to_diff <- df[, c('ID', 'time',\n 'wind_velocity', 'population', 'highway',\n 'no2_kg', 'AESPR__Jobos_mwh', 'Aguirre_mwh',\n 'Cambalache_mwh',\n 'Central_SJ__Palo_Seco_mwh',\n 'Daguao_mwh',\n 'EcoEl__Costa_Sur_mwh',\n 'HIMA_San_Pablo_Caguas_mwh',\n 'Janssen_Ortho_LLC_mwh', 'Mayaguez_mwh', \n 'Vega_Baja_mwh', \n 'Yabucoa_mwh', 'Culebra_mwh', 'Vieques_mwh')]\n# Computing differenced variables\nDiff <- function(varname){\n varname - lag(varname)\n}\ndf_diff <- vars_to_diff %>% \n group_by(ID) %>%\n mutate_at(vars(-c(group_cols(), time, \n wind_velocity,\n population, highway)), Diff) %>% drop_na()\n\n# Changing names\ncolnames(df_diff)[6:ncol(df_diff)] <- paste0(colnames(df_diff)[6:ncol(df_diff)], '_diff')\n\n##############################################################################################\n######################################## MEF ############################################\n##############################################################################################\n\n# Running random effects with lag and with error\nfm_diff <- no2_kg_diff ~ wind_velocity*(AESPR__Jobos_mwh_diff +\n Aguirre_mwh_diff +\n Cambalache_mwh_diff + \n Central_SJ__Palo_Seco_mwh_diff +\n Daguao_mwh_diff +\n EcoEl__Costa_Sur_mwh_diff +\n Janssen_Ortho_LLC_mwh_diff + \n Mayaguez_mwh_diff +\n Vega_Baja_mwh_diff +\n Yabucoa_mwh_diff +\n Culebra_mwh_diff) + population + highway\n\nspgm_diff <- spgm(formula = fm_diff, data = df_diff, listw = spatial_mat_list, lag = TRUE, \n moments = \"fullweights\", model = \"random\", spatial.error = TRUE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Table below shows the results of the MEF calculation. It can be seen that there is only one statistically significant slope, reflecting the MEF for the Cambalache power station, we report this as the MEF detectable for Puerto Rico.\n\nOverall, the devised specification seems adequate for the MEF estimation. However, to boost the efficiency of the MEF calculation one can use a more nuanced data of plants' activity and NO2 emission by, for example, hours (we did not have that at our disposal)."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# Printing the output for the MEF\nspgm_MEF_table <- round(as.data.frame(summary(spgm_diff)[[\"CoefTable\"]]),4)\ncolnames(spgm_MEF_table) <- c(\"Estimate\", \"SE\", \"t-value\", \"p-value\")\n\n# You can get the results bellow by running the following script in RStudio\n# kable(spgm_MEF_table, \"html\") %>%\n# kable_styling(c(\"striped\"), full_width = F) %>%\n# row_spec(c(7,12), bold = T, background = \"lightcoral\")%>%\n# row_spec(c(20, 22, 25), bold = T,\n# background = \"khaki\")\n\n# Simplier version of the coefficient \n# kable(spgm_MEF_table, \"rst\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"![](\"https://sun9-70.userapi.com/c857236/v857236853/12f889/0OMFi4pT34I.jpg\")
"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"MEF <- round(spgm_MEF_table$Estimate[6],2)\ncat('The MEF for Puerto Rico:', MEF, 'kg/mwh')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Discussion of the Method"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"**Advantages of the model and its improvements compared to the existing methodology:**\n\n* the crucial advantage is that the model accounts for both spatial and time-related dimensions of the data in a single system, therefore, we can effectively decompose emission factors between plants and over time;\n\n* the possibility to control for any factor that we can operationalize and measure such as population, traffic, factories, weather, etc., hence, any types of area (not only islands) are suitable for the exploration;\n\n* rather high computational speed for reasonable data sizes;\n\n* the ability to define the emission factors for smaller time slices (through interaction effects);\n\n* easiness of interpretation (like in a regular liner model);\n\n* the possibility to account for the heterogeneity of plants (as we plug their energy generation separately) coupled with coherence and integrity (as it is a single model);\n\n* easy scalability of the model.\n\n\n**Disadvantages of the model:**\n\n* the emission factors for electricity plants that do not give off a tangible amount of NO2 might not be captured, nevertheless, there is no guarantee any other approach can ever do that accurately;\n\n* the model requires relatively large RAM for producing parameter estimates.\n\n**Spatial Panel Data Modeling that we suggested for the calculation of the emissions factor requires several important assumptions/recommendations to consider in order to enlarge the precision of computations:**\n\n* electricity plants should be powerful and considerably contribute to NO2 emission, otherwise, the respective effect might be hard to register;\n\n* the method will benefit if the analyzed plants are geographically separable from each other;\n\n* it is highly recommended to utilize smaller time units, for example, hours; in this case, only a few samples of hourly data by, say, months might be enough to sensibly capture the variation of the data over the year;\n\n* NO2 generation should be independent of the wind since we extrapolate the registered effect from windless to windy periods (the ideal situation is when the wind is low or moderate).\n\n**Additional useful datasets/factors/assumptions:**\n\n* the data on factories and manufacturing capacity as another source of chemical emissions can be useful to control for;\n\n* the data on electricity consumption can be leveraged to add a time-related variability to the electricity production data;\n\n* in case of the absence of information about some plants in a dataset it is possible to extrapolate the estimated effects between similar plants in terms of a fuel type, generation capacity, etc. "},{"metadata":{},"cell_type":"markdown","source":"### References\n\nEIA-923 Puerto Rico Monthly Generation and Fuel Consumption Time Series File, 2018 Final Revision\n\nEIA Puerto Rico Sales of Electricity to Ultimate Customers: Total by End-Use Sector, by Month, 2016 - December 2019\n\nBaltagi, B. H., Song, S. H., & Koh, W. (2003). Testing panel data regression models with spatial error correlation. Journal of econometrics, 117(1), 123-150.\n\nMillo, G., & Piras, G. (2012). splm: Spatial panel data models in R. Journal of Statistical Software, 47(1), 1-38.\n\nSiler-Evans, K., Azevedo, I. L., & Morgan, M. G. (2012). Marginal emissions factors for the US electricity system. Environmental science & technology, 46(9), 4742-4748."}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":4}