MCDA-EVSE.Rmd

---
title: "Advancing Sustainable Mobility: A GIS-Based Multi-Criteria Decision Analysis for Equitable Electric Vehicle Supply Equipment Deployment in Philadelphia"
author: "Emily Zhou, Junyi Yang"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: simplex
    toc: yes
    toc_float: yes
    code_folding: hide
    code_download: yes
editor_options:
  markdown:
    wrap: sentence
---

Version 5.0 \| First Created Nov 17, 2023 \| Updated Mar 9, 2025

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

# Abstract 

The transition from traditional fossil fuel vehicles to electric vehicles (EVs) stands out as a pivotal solution to decarbonize the transportation systems and combat the climate crisis. However, the widespread adoption of EVs faces challenges, with one significant obstacle being the expansion of electric vehicle supply equipment (EVSE). The City of Philadelphia’s Office of Innovation Technology (OIT) is actively engaged in developing a network of extensive, equitable, and accessible EVSEs as part of its Smart City initiatives. Despite these efforts, the initiative faces a common challenge related to the selection of suitable sites.

We share a geographic information system (GIS)-based multi-criteria decision analysis (MCDA) method that can be used to evaluate the suitability of potential EVSE sites to support the sustainable and equitable deployment of EVSE in Philadelphia. Our MCDA approach considers key criteria ranging from socio-demographic indicators (e.g. driving-age population) to site-specific characteristics (e.g. spatial accessibility of existing EVSEs, availability of public parking garages, and city-wide power grid) and is based on the Analytic Hierarchy Process (AHP). To identify the optimal sites, three MCDA methods – WSM (weighted sum method), PROMETHEE (Method for Organizing Ranking of Preferences for Enrichment), as well as TOPSIS (Technique for Order Preference by Similarity to Ideal Solutions) – are applied and compared. 

We implemented our study with open-source R markdown and acquired data from the American Community Survey (ACS), OpenStreetMap, and the US Department of Energy. Specifically, a five-step solution approach is developed for the problem: 1) creating a fishnet for Philadelphia as the spatial unit for analysis and removing any water features, 2) determining and preprocessing criteria, among which we calculated the spatial accessibility of EVSEs using 2SFCA and distance to parking lots using k-nearest neighbor, 3) aggregating criteria into the fishnet, 4) prioritizing the criteria using AHP and finally 5) ranking the potential sites using WSM, PROMETHEE, and TOPSIS. 

The results of our MCDAs highlighted several areas in Philadelphia with a pronounced demand for new EVSE. We selected one of those sites in South Philadelphia and proposed a public-private partnership model between Philadelphia’s OIT and local grocery stores to install and maintain the EVSEs. We subsequently conducted financial analyses for the cost and revenue of breakdowns and designed a phased implementation of EVSE infrastructure given the current site conditions. Our study highlighted the challenges of agreeing on the input criteria and weighting schemes in MCDA but provided a spatially informed starting point to identify sites for new EVSE that contribute to Philadelphia’s sustainable transportation goals while addressing social disparities. It also provides a scalable and replicable model for other cities facing similar challenges in deploying EVSE and advancing smart city initiatives. More importantly, it emphasized the potential of geospatial analysis in shaping a more climate-smart and equitable future. 

The study is conducted collaboratively by graduate students from the Weitzman School of Design and Wharton Business School at the University of Pennsylvania and is available as a reproducible repository on [GitHub](https://github.com/emilyzhou112/CAGIS-UCGIS-2024).

## Keywords

Multi-Criteria Decision Analysis, Analytic Hierarchy Process, Electric Vehicle Supply Equipment, Sustainable Transportation, Smart City, GIS

# Setup

```{r packages and processing environment, message=FALSE, warning=FALSE}

# list of packages required 
packages = c("tidycensus", "tidyverse", "viridis", "FNN",  "dplyr", "sf", "classInt", "readr", "ggplot2", "here", "tmap", "SpatialAcc", "svDialogs", "MCDA", "nngeo", "leaflet", "gghalves", "NbClust", "kableExtra", "factoextra", "irlba", "devtools", "cluster", "NbClust")

# load and install required packages
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE, quietly=TRUE)
      library(x, character.only = TRUE)
    }
  }
)

# load functions
source(here("code", "KNN.R"))

# save the R processing environment
writeLines(
  capture.output(sessionInfo()),
  here("environment", paste0("r-environment-", Sys.Date(), ".txt"))
)
```

# Study Area

```{r process study area, warning=FALSE, message=FALSE}

philly <- st_read(here("data", "raw", "city_boundary.geojson")) %>% st_transform('ESRI:102728')
water <- st_read(here("data", "raw", "water_features.geojson")) %>% st_transform('ESRI:102728')

# generate fishnet for Philadelphia
fishnet <- st_make_grid(philly,
                        cellsize = 1640, # 500 meters
                        square = TRUE) %>% 
  .[philly] %>%           
  st_sf() %>%
  mutate(uniqueID = 1:n())

# remove grids in water with 330 ft buffer
fishnet_nowater <- fishnet %>%
  filter(!(uniqueID %in% (water %>%
                            st_buffer(dist = 330) %>%
                            st_intersection(st_centroid(fishnet), .) %>%
                            st_drop_geometry() %>%
                            dplyr::select(uniqueID) %>%
                            pull(uniqueID))))

```


```{r map study grids}
ggplot() +
  geom_sf(data=fishnet, color="black", fill="#726DA8") +
  geom_sf(data=fishnet_nowater, color="black", fill="white") +
  labs(title = "Fishnet of Philly") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth=0.8)
        )
```

```{r save the intermediaries}

st_write(fishnet, here("data", "derived", "philly_fishnet.geojson"), driver = "GeoJSON")
st_write(fishnet_nowater, here("data", "derived", "fishnet_nowater.geojson"), driver = "GeoJSON")

```

# Criteria Processing 

## Demographic Information and EV Ownership

```{r load demographic info, message=FALSE, warning=FALSE}

options(timeout=1000) 

# get API Key
census_api_key(dlgInput("Enter a Census API Key", 
   Sys.getenv("CENSUS_API_KEY"))$res,
   overwrite = TRUE)
acs_vars <- load_variables(year = 2022, dataset = "acs5", cache = TRUE)

# query 2022 acs data
philly22 <- get_acs(geography = "tract", 
          variables = c(
            "B01001_001E", # total population
            "B01001_010E", # male 22-24
            "B01001_011E",
            "B01001_012E",
            "B01001_013E",
            "B01001_014E",
            "B01001_015E",
            "B01001_016E",
            "B01001_017E",
            "B01001_018E",
            "B01001_019E", # male 62-64
            "B01001_034E",
            "B01001_035E",
            "B01001_036E",
            "B01001_037E",
            "B01001_038E",
            "B01001_039E",
            "B01001_040E",
            "B01001_041E",
            "B01001_042E",
            "B01001_043E" #female 62-64
            ), 
          year=2022, state="PA", county="Philadelphia", 
          geometry=TRUE, output="wide") %>%
          st_transform('ESRI:102728')

philly22 <- philly22 %>%
  mutate(popden = B01001_001E / (st_area(geometry) * 9.2903e-8),
         density = as.integer(gsub("\\[.*\\]", "", popden)),
         pop22_64 = B01001_010E + B01001_011E + B01001_012E + B01001_013E + B01001_014E +
                    B01001_015E + B01001_016E + B01001_017E + B01001_018E + B01001_019E +
                    B01001_034E + B01001_035E + B01001_036E + B01001_037E + B01001_038E + 
                    B01001_039E + B01001_040E + B01001_041E + B01001_042E + B01001_043E,
         driving_pop = pop22_64 / B01001_001E,
         totalpop = B01001_001E) %>%
  dplyr::select(GEOID, totalpop, density, driving_pop) %>%
  mutate(driving_pop = ifelse(is.nan(driving_pop), 0, driving_pop))


```


```{r save the queried acs information}

st_write(philly20, here("data", "derived", "acs2020.geojson"), driver = "GeoJSON")

```


```{r load registered ev, message=FALSE, warning=FALSE}

reg_ev <- st_read(here("data", "raw", "reg_ev.geojson"))
Zip_Code <- as.character(c(
  19120, 19124, 19111, 19143, 19149, 19134, 19140, 19148, 19104, 19144,
  19145, 19131, 19139, 19146, 19147, 19115, 19136, 19128, 19135, 19121,
  19154, 19141, 19132, 19152, 19114, 19116, 19151, 19138, 19142, 19119,
  19130, 19125, 19133, 19103, 19150, 19122, 19126, 19123, 19107, 19106,
  19153, 19129, 19118, 19137, 19127, 19102, 19108, 19109, 19176, 19112,
  19160, 19155, 19162, 19161, 19171, 19170, 19173, 19172, 19175, 19178,
  19177, 19181, 19179, 19183, 19182, 19185, 19184, 19188, 19187, 19192,
  19191, 19194, 19193, 19196, 19195, 19244, 19197, 19255, 19019, 19093,
  19092, 19101, 19099, 19105, 19110, 19190
))
reg_ev <- reg_ev %>% 
  filter(zip %in% Zip_Code) %>% 
  dplyr::select(elec_cnt) %>% 
  st_transform('ESRI:102728')

```


## Existing EVSE Distribution and Access

```{r load existing evse, message=FALSE, warning=FALSE}

evse <- read.csv(here::here("data", "raw", "existing_evse.csv"))
evse <- evse %>% 
  filter(City == "Philadelphia") %>% 
  dplyr::select(Latitude, Longitude, EV.Level2.EVSE.Num) %>% 
  filter(is.na(Latitude) == FALSE & is.na(Longitude) == FALSE) %>% 
  st_as_sf(., coords = c("Longitude", "Latitude"), crs = 4326) %>% 
  st_transform('ESRI:102728') %>% 
  mutate(num = ifelse(is.na(EV.Level2.EVSE.Num), 1, EV.Level2.EVSE.Num))

```


```{r compute spatial accessibility, message=FALSE, warning=FALSE}

tracts.coords <- st_coordinates(st_centroid(philly22))
evse.coords <- st_coordinates(evse)
dist.matrix <- distance(tracts.coords, evse.coords, type = "euclidean")

TSFCA <- ac(p = philly22$totalpop, 
            n = evse$num, 
            D = dist.matrix, d0 = 2000, family = "2SFCA")
philly22 <- philly22 %>% 
  mutate(TSFCA = TSFCA)
```


## Local Environment

```{r load site environmental metrics, message=FALSE, warning=FALSE}

parking <- st_read(here("data", "raw", "osm_features", "parking.geojson"))
parking <- parking %>% 
  filter(!access %in% c("no", "permissive", "permit", "private")) %>% 
  st_transform('ESRI:102728')

police_fire <- st_read(here("data", "raw", "osm_features", "police_fire.geojson")) 
police_fire <- police_fire %>% 
  st_transform('ESRI:102728') %>% 
  st_centroid()

zoning <- st_read(here("data", "raw", "zoning", "zoning.shp"))
zoning <- zoning %>% 
  st_transform('ESRI:102728') %>% 
  mutate(zoning = case_when( ZONINGGROU == "Commercial/Commercial Mixed-Use" ~ "Commercial",
                             ZONINGGROU == "Industrial/Industrial Mixed-Use" ~ "Industrial",
                             ZONINGGROU == "Residential/Residential Mixed-Use" ~ "Residential",
                             ZONINGGROU == "Special Purpose" ~ "Special",
                             TRUE ~ "Undefined")) %>% 
  mutate(zone_score = case_when( zoning == "Commercial" ~ 0.4,
                             zoning == "Industrial" ~ 0.15,
                             zoning == "Residential" ~ 0.3,
                             TRUE ~ 0.15))
```

## Energy

```{r load property electricity data, message=FALSE, warning=FALSE}

electricity <- read.csv(here("data", "raw", "property-electricity-data.csv")) %>% 
  dplyr::select(Y = y_coord, X = x_coord, portfolio_bldg_id, postal_code, primary_property_type, sector, electricity_2021) %>%
  na.omit() %>%
  st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
  st_transform('ESRI:102728') %>% 
  st_intersection(., fishnet) # clip to fishnet extent
```

# Aggregate Criteria

```{r aggregate all criteria into fishnet, message=FALSE, warning=FALSE}

net_centroid <- st_centroid(fishnet_nowater)
fishnet_nowater <- fishnet_nowater %>%
  left_join(net_centroid %>% 
              mutate(evse.nn = nn_function(st_coordinates(net_centroid), 
                                           st_coordinates(evse), 2)*0.3048) %>% # distance to evse
              st_drop_geometry(), by = "uniqueID") %>% 
  left_join(net_centroid %>% 
              st_intersection(philly22) %>% 
              st_drop_geometry(), by = "uniqueID") %>% # demographic info and spatial acc
  left_join(net_centroid %>% 
              st_intersection(reg_ev) %>% 
              st_drop_geometry(), by = "uniqueID") %>% 
  left_join(net_centroid %>% 
              mutate(parking.nn = nn_function(st_coordinates(net_centroid), 
                                           st_coordinates(parking), 2)*0.3048) %>% # dist to parking
              st_drop_geometry(), by = "uniqueID") %>% 
  left_join(net_centroid %>% 
              mutate(policefire.nn = nn_function(st_coordinates(net_centroid), 
                                           st_coordinates(police_fire), 1)*0.3048) %>% # police and fire
              st_drop_geometry(),by = "uniqueID") %>% 
  left_join(net_centroid %>% 
              st_buffer(dist = 200) %>% 
              st_intersection(zoning) %>% 
              st_drop_geometry(), by = "uniqueID") %>% 
  mutate(TSFCA = ifelse(is.infinite(TSFCA), 0, TSFCA)) %>% 
  mutate(elec_cnt = ifelse(is.na(elec_cnt), 0, elec_cnt)) %>% 
  filter(is.na(GEOID) == FALSE) %>% 
  filter(is.na(zone_score) == FALSE) %>% 
  dplyr::select(-c(GEOID, totalpop, featCount, ZONINGGROU, zoning)) %>% 
  group_by(uniqueID) %>% # for duplicate grids, take the mean of the zone score
  mutate(zone_score = mean(zone_score)) %>%
  ungroup() %>%
  distinct(uniqueID, .keep_all = TRUE)
  
    
```


```{r aggregate power grid, message=FALSE, warning=FALSE}

# buffer the building and aggregate building centroid into the buffer
# if outside the buffer, get the nearest neighbor
elec_grid <- st_join(net_centroid, 
                electricity %>% st_buffer(dist = 1640)) %>% 
  group_by(uniqueID.x) %>% 
  summarize(elect_use = mean(electricity_2021))

unjoin <- elec_grid %>% 
  filter(is.na(elect_use)) %>% 
  mutate(neighbor_index = 1:n()) 
  
nn_indices <- st_nn(unjoin,
                    electricity, k = 1)
nearest_neighbors <- electricity[unlist(nn_indices), ] %>% mutate(neighbor_index = 1:n())

unjoin <- unjoin %>%
  st_drop_geometry() %>% 
  left_join(nearest_neighbors %>% select(electricity_2021, neighbor_index), 
            by = c("neighbor_index" = "neighbor_index"))

elec_grid <- elec_grid %>% 
  left_join(., unjoin %>% 
              dplyr::select(uniqueID.x, electricity_2021), by = c("uniqueID.x" = "uniqueID.x")) %>% 
  mutate(elect_use = case_when(
    is.na(elect_use) ~ electricity_2021,
    TRUE ~ elect_use
  )) %>% 
  dplyr::select(-electricity_2021) %>% 
  st_drop_geometry()

fishnet_nowater <- fishnet_nowater %>% 
  left_join(., elec_grid, by = c("uniqueID" = "uniqueID.x"))
```


```{r save all aggregated criteria}

st_write(fishnet_nowater, here("data", "derived", "all_criteria_raw.geojson"), driver = "GeoJSON")

```



# Multi-Criteria Decision Analysis

In it's simplest term, MCDA can be summarized into three steps.
1. Transform indicator data into standard **geographic unit**. This include: aggregate or dissolve for nested relationship, area-weighted re-aggregation for un-nested relationship, count points within polygons, etc. 
2. Transform indicator data into standard **measuring unit**. This include: rank or percentile, z-score, some other functions. 
3. Combine indicator data into a composite score with certain weights. 

## Self-Assigned Weight + WSM

The Weighted Sum Method is one of the simplest and most commonly used MCDA techniques. In WSM, each criterion is assigned a weight based on its importance, and each alternative is scored based on these criteria. The final score for each alternative is calculated by summing the products of the scores and their respective weights.

1. Scale all the criteria from 0 to 1. 
2. Invert the scale for certain criteria if necessary. 
3. Assign a weight to each criterion reflecting its relative importance. The sum of weights should equal 1.
4. For each alternative, multiply the performance score of each criterion by its weight and sum the results.
5. Rank the alternatives based on their total weighted scores.

```{r weighted sum}

scale_values <- function(x){(x-min(x))/(max(x)-min(x))}
WSM <- fishnet_nowater %>% 
  mutate(scl_density = scale_values(density),
         scl_acc = scale_values(TSFCA),
         scl_drivingpop = scale_values(driving_pop),
         scl_evse = scale_values(evse.nn),
         scl_elec = scale_values(elec_cnt),
         scl_parking = scale_values(parking.nn),
         scl_policefire = scale_values(policefire.nn),
         scl_elect_use = scale_values(elect_use)
         ) %>% # now inverse some scales
  mutate(scl_evse_re = 0 - scl_evse + 1,
         scl_acc_re = 0 - scl_acc + 1, 
         scl_parking_re = 0 - scl_parking + 1,
         scl_policefire_re = 0 - scl_policefire + 1) %>% 
  mutate(score = 0.12*scl_density + 0.18*scl_evse_re + 0.15*scl_elec + 0.1*scl_drivingpop + 0.1*scl_parking_re + 0.05*scl_policefire_re + 0.05*zone_score + 0.12*scl_acc_re +0.13*scl_elect_use) # self assign weight
```


```{r visualize weighted sum}

custom_palette <- c("#C4C4C4", "#B2BF95", "#80A676", "#88BFBF", "#7EA1BF")

quantiles <- classIntervals(WSM$score, n = 5, style = "quantile")
WSM$quantile <- cut(WSM$score, breaks = quantiles$brks, include.lowest = TRUE)

WSM_plot <- ggplot()+
  geom_sf(data=fishnet, color=NA, fill="white") +
  geom_sf(data=WSM, color = NA, aes(fill=quantile)) + 
  scale_fill_manual(values = custom_palette) + 
  labs(title = "Weighted Sum Method Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )
```


```{r save a version}

st_write(WSM, here("data", "derived", "WSM.geojson"), driver = "GeoJSON")

```


## Self Assigned Weight + TOPSIS

TOPSIS is a method that identifies solutions from a finite set of alternatives based on their geometric distance from an ideal solution. The ideal solution has the best performance values for all criteria, and the negative-ideal solution has the worst. This approach is employed widely for four main reasons: 1) the logic is rational and understandable, 2) the computation process is straightforward, 3) the concept permits the pursuit of the best alternatives for each criterion depicted in a simple mathematical form, 4) the importance weights are incorporated into the comparison procedure. 

1. Construct the decision matrix where each element represents the performance score of each alternatives.
2. Normalization
3. Determine the ideal and negative ideal solution. In other words, the best and worst possible value for each criterion. 
4. Calculate the separation measure as in euclidean distance calculation. 
5. Measure relative closeness: how close each alternative is to the ideal solution relative to its distance from the negative-ideal solution. 
6. Ranking 

```{r construct performance table}

performanceTable <- fishnet_nowater %>% 
  dplyr::select(-c(uniqueID)) %>% 
  st_drop_geometry() %>% 
  mutate(density = scale_values(density),
         TSFCA = scale_values(TSFCA),
         driving_pop = scale_values(driving_pop),
         evse.nn = scale_values(evse.nn),
         elec_cnt = scale_values(elec_cnt),
         parking.nn = scale_values(parking.nn),
         policefire.nn = scale_values(policefire.nn),
         elect_use = scale_values(elect_use)
         )
uniqueIDs <- fishnet_nowater$uniqueID
rownames(performanceTable) <- uniqueIDs
```


```{r topsis operations}

# reverse is min
weights <- c(0.18, 0.12, 0.1, 0.12, 0.15, 0.1, 0.05, 0.05, 0.13)
criteriaMinMax <- c("min", "max", "max", "min", "max", "min", "min", "max", "max")
# positiveIdealSolutions <- c(0.17, 0.60, 0.70, 0, 0.50, 0.10, 0.14, 0.3,0.007)
overall1 <- TOPSIS(performanceTable, weights, criteriaMinMax)
TOPSIS <- data.frame(uniqueID = names(overall1), values = overall1) %>%  
  mutate(uniqueID = as.integer(uniqueID))
  
TOPSIS_toplot <- fishnet_nowater %>% 
  left_join(TOPSIS, by = "uniqueID")

```


```{r visualize topsis result}

quantiles <- classIntervals(TOPSIS_toplot$values, n = 5, style = "quantile")
TOPSIS_toplot$quantile <- cut(TOPSIS_toplot$values, breaks = quantiles$brks, include.lowest = TRUE)

TOPSIS_plot <- ggplot()+
  geom_sf(data=fishnet, color=NA, fill="white") +
  geom_sf(data=TOPSIS_toplot, color=NA, aes(fill=quantile))+ 
  scale_fill_manual(values = custom_palette) + 
  labs(title = "TOPSIS Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )

```

However, the TOPSIS method presents certain drawbacks. One of the problems attributable to TOPSIS is that it can cause the phenomenon known as rank reversal. In this phenomenon the alternatives’ order of preference changes when an alternative is added to or removed from the decision problem. In some cases this may lead to what is called total rank reversal, where the order of preferences is totally inverted, that is to say, that the alternative considered the best, with the inclusion or removal of an alternative from the process, then becomes the worst. 


```{r save a version}

st_write(TOPSIS_toplot, here("data", "derived", "TOPSIS.geojson"), driver = "GeoJSON")

```


## Self Assigned Weights + PROMETHEE

PROMETHEE is a ranking method based on pairwise comparisons of alternatives according to each criterion, considering both preference and indifference thresholds. Preference threshold refers to the minimum difference between the performance of two alternatives on a given criterion such that the decision-maker has a strict preference for one alternative over the other. **If the performance difference is greater than or equal to p, the decision-maker clearly prefers one alternative over the other for that criterion.** The indifference threshold q represents the maximum difference between the performance of two alternatives on a given criterion such that the decision-maker is indifferent between them. **If the performance difference is less than or equal to q, the decision-maker considers the two alternatives to be effectively equal with respect to that criterion.** A preference function, which define how much better one site is over another for each criterion (a linear function where the lower the accessibility, the better), is then applied based on the two thresholds. 

1. Calculate the difference in performance scores for each criterion between all pairs of alternatives.
2. Set up an indifference threshold.
3. Set up a preference threshold. 
4. Decide a preference function
5. Apply the preference function using the indifference and preference threshold to convert the differences into preference value, which is typically betweeen 0 and 1. 
6. For each alternative, sum the preference values across all criteria to obtain the aggregated preference indices.
7. Calculate the positive and negative outranking flows for each alternative, which represent how much an alternative outranks others and is outranked by others, respectively.
8. Use the net outranking flow (positive flow minus negative flow) to rank the alternatives.

Now, this indifference and preference thresholds are based solely on quantitative analysis on the distribution of the data. We computed the IQR for all criteria and use that as the preference threshold. We divide the IQR by three and use that as the indifference threshold. 

```{r determine indifference and preference threshold}

IQR(performanceTable$evse.nn) / 3
IQR(performanceTable$density) / 3
IQR(performanceTable$driving_pop) /3 
IQR(performanceTable$TSFCA) / 3 
IQR(performanceTable$elec_cnt) / 3
IQR(performanceTable$parking.nn) / 3
IQR(performanceTable$policefire.nn) / 3
IQR(performanceTable$zone_score) / 3
IQR(performanceTable$elect_use) / 3

```


```{r promethee operations}

criteriaWeights <- c(0.18, 0.12, 0.1, 0.12, 0.15, 0.1, 0.05, 0.05, 0.13)
names(criteriaWeights)<-colnames(performanceTable)

criteriaMinMax <- c("min", "max", "max", "min", "max", "min", "min", "max", "max")
names(criteriaMinMax)<-colnames(performanceTable)

preferenceFunction<-c("Usual","U-shape","V-shape","Level","V-shape-Indiff","Gaussian", "Level","V-shape-Indiff", "Level")
# The gaussParameter vector in the code specifies the standard deviation (s) for the Gaussian preference 
# function for each criterion.
gaussParameter<-c(0.25,1,2,0,0,0,0,0,0)
names(gaussParameter)<-colnames(performanceTable)

#Preference threshold
preferenceThreshold<-c(0.31, 0.21, 0.12, 0, 0.17, 0.19, 0.23, 0.075, 0.005367344)
names(preferenceThreshold)<-colnames(performanceTable)

#Indifference threshold
indifferenceThreshold<-c(0.1, 0.07, 0.041, 0, 0.058, 0.065, 0.078, 0.025, 0.001789115)
names(indifferenceThreshold)<-colnames(performanceTable)

performanceTable <- as.matrix(performanceTable)
promethee <- PROMETHEEOutrankingFlows(performanceTable, preferenceFunction,preferenceThreshold,
indifferenceThreshold,gaussParameter,criteriaWeights,criteriaMinMax)
```


```{r clean up promethee table}

promethee_df <- data.frame(uniqueID = names(promethee$outrankingFlowsPos), 
           Pos = unname(promethee$outrankingFlowsPos),
           Neg = unname(promethee$outrankingFlowsNeg)) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(net = Pos-Neg) %>% 
  dplyr::select(uniqueID, net)

PROMETHEE_toplot <- fishnet_nowater %>% 
  left_join(promethee_df, by = "uniqueID")
  
```


```{r visualize promethee}
quantiles <- classIntervals(PROMETHEE_toplot$net, n = 5, style = "quantile")
PROMETHEE_toplot$quantile <- cut(PROMETHEE_toplot$net, breaks = quantiles$brks, include.lowest = TRUE)

PROMETHEE_plot <- ggplot()+
  geom_sf(data=fishnet, color=NA, fill="white") +
  geom_sf(data=PROMETHEE_toplot, color=NA, aes(fill=quantile)) + 
  scale_fill_manual(values = custom_palette) + 
  labs(title = "PROMETHEE Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )


```


```{r save a version}

st_write(PROMETHEE_toplot, here("data", "derived", "PROMETHEE.geojson"), driver = "GeoJSON")

```

## AHP + TOPSIS

To determine weights using the Analytic Hierarchy Process (AHP), start by defining the decision problem and structuring it into a hierarchy, including the goal at the top, followed by criteria and sub-criteria (if any), and the alternatives at the bottom. Construct pairwise comparison matrices for the criteria by comparing each pair and assigning relative importance values on a scale from 1 to 9. Normalize the matrices by dividing each element by the sum of its column, then calculate the priority vector (weights) by averaging the normalized values across each row. 

```{r set up pairwise matrix}

crit <- c("density","evse","ev", "drivingpop", "parking", "policefire", "zoning", "acc", "electricity")

criteriaWeightsPairwiseComparisons <- matrix(c(1.0, 3.0, 2.0, 1/2, 1/2, 1/4, 1/3, 1.5, 1.2,
                                               1/3, 1.0, 1/1.2, 1/4, 1/4, 1/8, 1/6, 1/2, 1/1.5, 
                                               1/2, 1.2, 1.0, 1/3, 1/3, 1/6, 1/6, 1/2, 1/1.5, 
                                               2, 4, 3, 1.0, 1/1.2, 1/4, 1/5, 1.5, 1.8,
                                               2, 4, 3, 1.2, 1.0, 1/4, 1/4, 2, 1.5, 
                                               4, 8, 6, 4, 4, 1, 1.2, 3, 2,
                                               3, 6, 6, 5, 4, 1/1.2, 1, 3, 2,
                                               1/1,5, 2, 2, 1/1.5, 1/2, 1/3, 1/3, 1, 1.5, 
                                               1/1.2, 1.5, 1.5, 1/1.8, 1/1.5, 1/2, 1/2, 1/1.5, 1
                                               ),
nrow=length(crit),
ncol=length(crit),
dimnames=list(crit,crit))


```


```{r compute weights based on the matrix}

column_sums <- colSums(criteriaWeightsPairwiseComparisons)

# Divide each element by the sum of its column
normalized_matrix <- sweep(criteriaWeightsPairwiseComparisons, 2, column_sums, FUN = "/")
weights <- rowMeans(normalized_matrix)

# Output the weights
print(weights)

```

```{r apply topsis operations}

# reverse is min
weights <- c(0.249, 0.114, 0.105, 0.098, 0.1909, 0.074, 0.033, 0.033, 0.102)
criteriaMinMax <- c("min", "max", "max", "min", "max", "min", "min", "max", "max")
overall2 <- TOPSIS(performanceTable, weights, criteriaMinMax)
AHP_TOPSIS <- data.frame(uniqueID = names(overall2), values = overall2) %>%  
  mutate(uniqueID = as.integer(uniqueID))
  
AHP_TOPSIS_toplot <- fishnet_nowater %>% 
  left_join(AHP_TOPSIS, by = "uniqueID") 


```


```{r visualize ahp results}

quantiles <- classIntervals(AHP_TOPSIS_toplot$values, n = 5, style = "quantile")
AHP_TOPSIS_toplot$quantile <- cut(AHP_TOPSIS_toplot$values, breaks = quantiles$brks, include.lowest = TRUE)

AHP_TOPSIS_plot <- ggplot()+
  geom_sf(data=fishnet, color=NA, fill="white") +
  geom_sf(data=AHP_TOPSIS_toplot, color=NA, aes(fill=quantile)) + 
  scale_fill_manual(values = custom_palette) +
  labs(title = "AHP Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )


```

```{r save a version}

st_write(AHP_TOPSIS_toplot, here("data", "derived", "AHP.geojson"), driver = "GeoJSON")

```


# Analyses

## Mapping Spatial Consistency 

```{r fig.height=12, fig.width=12}

library(ggpubr)
ggarrange(WSM_plot, TOPSIS_plot, PROMETHEE_plot, AHP_TOPSIS_plot, ncol = 2, nrow = 2)

```

## Measure Rank Alignment 

```{r combine all ranks together}

TOPSIS <- TOPSIS %>%
  mutate(topsis_rank = min_rank(desc(values)))

WSM <- WSM %>% 
  dplyr::select(c(uniqueID, score)) %>% 
  mutate(wsm_rank = min_rank(desc(score)))

PROMETHEE <- promethee_df %>% 
  mutate(promethee_rank = min_rank(desc(net)))

AHP_TOPSIS <- AHP_TOPSIS %>% 
  mutate(ahptopsis_rank = min_rank(desc(values)))
 
rank_compare <- WSM %>% 
  dplyr::select(c(uniqueID, wsm_rank)) %>% 
  left_join(., TOPSIS %>% dplyr::select(uniqueID, topsis_rank), by = "uniqueID") %>% 
  left_join(., PROMETHEE %>% dplyr::select(uniqueID, promethee_rank), by = "uniqueID") %>% 
  left_join(., AHP_TOPSIS %>% dplyr::select(uniqueID, ahptopsis_rank), by = "uniqueID") %>% 
  mutate(avg_rank = as.integer((wsm_rank + topsis_rank + promethee_rank + ahptopsis_rank) / 4))

```


```{r spearman rank correlation}

options(scipen=999)
# Define a function to compute correlation and p-value
compute_spearman <- function(x, y) {
  cor_test <- cor.test(x, y, method = "spearman")
  return(c(correlation = cor_test$estimate, p_value = cor_test$p.value))
}

# Compute correlations and p-values for all comparisons
global_cor_results <- data.frame(
  Comparison = c("AHP vs TOPSIS", "AHP vs WSM", "AHP vs PROMETHEE",
                 "TOPSIS vs WSM", "TOPSIS vs PROMETHEE", "WSM vs PROMETHEE"),
  Correlation = c(
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$topsis_rank)[1],
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$wsm_rank)[1],
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$promethee_rank)[1],
    compute_spearman(rank_compare$topsis_rank, rank_compare$wsm_rank)[1],
    compute_spearman(rank_compare$topsis_rank, rank_compare$promethee_rank)[1],
    compute_spearman(rank_compare$wsm_rank, rank_compare$promethee_rank)[1]
  ),
  P_Value = c(
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$topsis_rank)[2],
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$wsm_rank)[2],
    compute_spearman(rank_compare$ahptopsis_rank, rank_compare$promethee_rank)[2],
    compute_spearman(rank_compare$topsis_rank, rank_compare$wsm_rank)[2],
    compute_spearman(rank_compare$topsis_rank, rank_compare$promethee_rank)[2],
    compute_spearman(rank_compare$wsm_rank, rank_compare$promethee_rank)[2]
  )
)
cor.test(rank_compare$ahptopsis_rank, rank_compare$topsis_rank, method = "spearman")

```


Comparing across TOPSIS, PROMETHEE, and weighted sum in ranking the alternatives by looking at the distribution of differences in the rank, we found that:
1. The differences in rank between weighted sum and promethee method is the smallest. 
2. TOPSIS is leading to rank reversal issues for some parts of Philadelphia. In other word, a few grids that were ranked of lower priority in PROMETHEE and WSM are ranked of much high priority in TOPSIS and vice versa. Closer examination of these grids reveal that they are located in the outskirt of Philadelphia, mainly industrial areas that use a lot of electricty power. TOPSIS assumes that criteria are independent of each other. When a new alternative is introduced or an existing one is removed, the distance to the ideal and negative-ideal solutions can change. The presence of extreme values (very high or very low) can significantly influence the ideal and negative-ideal solutions. MCDA is sensitive to the quality of our data. 
3. Assign weights directly has under-ranked several sites than using the AHP. 

```{r neighborhood?}

# rank consistency by neighborhood? 
neighborhood <- st_read(here("data", "raw", "PhillyPlanning_Neighborhoods", "PhillyPlanning_Neighborhoods.shp")) %>% st_transform('ESRI:102728') 
  

final <- fishnet_nowater %>% 
  left_join(., rank_compare %>% st_drop_geometry, by = "uniqueID") %>% 
  left_join(net_centroid %>% 
               st_intersection(neighborhood %>% dplyr::select(NAME)) %>% 
              st_drop_geometry(), by = "uniqueID") %>% 
  mutate(NAME = ifelse(is.na(NAME), "NOT APPLICABLE", NAME))

```

1. Difference between WSM and TOPSIS
2. Difference between WSM and PROMETHEE
3. Difference between WSM and AHP
4. Difference between TOPSIS and PROMETHEE
**5. Difference between TOPSIS and AHP**
6. Difference between PROMETHEE and AHP

```{r}

# if negative, it means that the first one over-rank the second one
# if positive, it means that the first one under-rank the second one
compare<- final %>% 
  mutate(wsm_topsis = wsm_rank - topsis_rank, # compare the rank between wsm and topsis
         wsm_promethee = wsm_rank - promethee_rank, # compare the rank between wsm and promethee
         wsm_ahp = wsm_rank - ahptopsis_rank, # compare the rank between wsm and ahp
         promethee_topsis = promethee_rank - topsis_rank, # compare the rank between topsis and promethee
         topsis_ahp = topsis_rank - ahptopsis_rank, # compare the rank between topsis and ahp
         promethee_ahp = promethee_rank - ahptopsis_rank)  # compare the rank between promethee and topsis


```


```{r compare difference distribution 2}

custom_palette1 <- c("#B2BF95", "#80A676", "#88BFBF", "#7EA1BF", "grey", "pink")
my_sort <- c("promethee_topsis", "topsis_ahp", "wsm_topsis", "wsm_promethee", "wsm_ahp", "promethee_ahp")

temp <- compare %>% 
  st_drop_geometry() %>% 
  dplyr::select(c(uniqueID, promethee_topsis, topsis_ahp, wsm_topsis, wsm_promethee, wsm_ahp, promethee_ahp, NAME)) %>%
  pivot_longer(cols = c("promethee_topsis", "topsis_ahp", "wsm_topsis", "wsm_promethee", "wsm_ahp", "promethee_ahp"),  
               names_to = "variable",        
               values_to = "value") %>% 
  ggplot(aes(x = as.numeric(factor(variable, levels = my_sort)))) +
  geom_half_violin(aes(x = as.numeric(factor(variable, levels = my_sort)) + 0.1,
                       y = value, fill = factor(variable, levels = my_sort)),
                   side = 'r', scale = "width", cex = 0.8, lwd=0.3) +  
  geom_boxplot(aes(x = as.numeric(factor(variable, levels = my_sort)) + 0.1,
                   y = value, fill = factor(variable, levels = my_sort)),
               outlier.colour = "black", width = 0.1, cex = 0.8, lwd=0.3, outlier.size=0.6) +  
  geom_jitter(aes(x = as.numeric(factor(variable, levels = my_sort)) - 0.1,
                  y = value, color = factor(variable, levels = my_sort)),
              width = 0.1, size = 0.2, stroke = 0.8) +
  scale_fill_manual(values = custom_palette1) + 
  scale_color_manual(values = custom_palette1) +  
  labs(title = "Distribution of Differences in Rank Between Different Methods",
       x = "Method",
       y = "Difference") +
  theme_bw()
  theme(axis.text.x=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_blank(),
        )


# my_sort <- c("topsis_promethee", "topsis_ahp", "wsm_topsis", "wsm_promethee", "wsm_ahp", "promethee_ahp")
ggsave("/Users/emzhou/Desktop/Paper/svg/violin.svg", temp, width = 10, height = 7)

```

```{r}
# map distribution of the criteria
final %>% 
  st_drop_geometry() %>%
  dplyr::select(c(density, TSFCA, driving_pop, evse.nn, elec_cnt, parking.nn, policefire.nn, elect_use, zone_score)) %>%
  pivot_longer(cols = c(density, TSFCA, driving_pop, evse.nn, elec_cnt, parking.nn, policefire.nn, elect_use, zone_score),  
               names_to = "variable",        
               values_to = "value") %>%
  ggplot(aes(x = value)) +
  geom_density(aes(fill = variable), alpha = 0.5) +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Distribution of Criteria",
       x = "Value",
       y = "Density")

summary(final$elect_use)
```


```{r }


final_sankey <- final %>%
  mutate(topsis_percentile = percent_rank(topsis_rank),
         wsm_percentile = percent_rank(wsm_rank),
         promethee_percentile = percent_rank(promethee_rank),
         ahp_percentile = percent_rank(ahptopsis_rank)) %>%
  mutate(
    topsis_cat = cut(topsis_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")),
    wsm_cat = cut(wsm_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")),
    promethee_cat = cut(promethee_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")),
    ahp_cat = cut(ahp_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%"))
  )
  
final_sankey %>% 
  dplyr::count(topsis_cat, wsm_cat, promethee_cat, ahp_cat) %>%
  ggplot(., aes(axis1 = topsis_cat, axis2 = ahp_cat, axis3 = promethee_cat, axis4 = wsm_cat, y = n)) +
  scale_x_discrete(limits = c("TOPSIS", "AHP", "PROMETHEE", "WSM"), 
                   expand = c(.05, .5)) +
  geom_alluvium(aes(fill = topsis_cat), alpha = 0.8) +  # Fill based on initial ranking
  geom_stratum(color = "grey") + 
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 3) + 
  scale_fill_manual(values = custom_palette) +
  theme_bw() +
  labs(title = "Rank Transitions Between Global and Local TOPSIS Methods",
       x = "Ranking Method", y = "Grid Count")

```

```{r}

# percent of identical positions in different decision-making methods
final_sankey %>% 
  dplyr::select(topsis_cat, wsm_cat, promethee_cat, ahp_cat) %>%
  #filter(promethee_cat == wsm_cat) %>%
  #filter(ahp_cat == wsm_cat) %>%
  # filter(promethee_cat == ahp_cat) %>%
  #filter(topsis_cat == ahp_cat) %>%
  #filter(topsis_cat == wsm_cat) %>%
  #filter(topsis_cat == promethee_cat) %>%
  filter(wsm_cat == ahp_cat & ahp_cat == promethee_cat) 

```

Pct of identical positions between TOPSIS and WSM: 772 / 1212 = 63.7%
Pct of identical positions between TOPSIS and PROMETHEE: 750 / 1212 = 61.9%
Pct of identical positions between WSM and PROMETHEE: 848 / 1212 = 70.0%
Pct of identical positions between TOPSIS and AHP-Topsis: 805 / 1212 = 66.4%
Pct of identical positions between WSM and AHP-Topsis: 696 / 1212 = 57.4%
Pct of identical positions between PROMETHEE and AHP-Topsis: 654 / 1212 = 54.0%
Pct of identical positions between TOPSIS, AHP-TOPSIS, and WSM: 568 / 1212 = 46.9%
Pct of identical positions between TOPSIS, AHP-TOPSIS, and PROMETHEE: 548 / 121 = 45.2%
Pct of identical positions between WSM, AHP-TOPSIS, and PROMETHEE: 545 / 1212 = 45.0%
Pct of identical positions between all four methods: 472 / 1212 = 38.9%



## Visualize Spatial Rank Variability 

```{r fig.height=7, fig.width=10}

custom_palette2 <- c("#e8623f", "white", "#57748d")


compare %>% 
  dplyr::select(c(uniqueID, promethee_topsis, topsis_ahp, wsm_topsis, wsm_promethee, wsm_ahp, promethee_ahp, NAME)) %>%
  pivot_longer(cols = c("promethee_topsis", "wsm_topsis",  "wsm_ahp", "promethee_ahp"),  
               names_to = "variable",       
               values_to = "value") %>% 
  # filter(variable == "promethee_topsis" | variable == "wsm_topsis") %>% 
  ggplot() +
  geom_sf(aes(fill = value), color = "grey") +
  labs(title = "Ranking Difference between TOPSIS and Other Approaches") +
  facet_wrap(~ variable, labeller= labeller(variable = c(
    `promethee_topsis` = "PROMETHEE vs. TOPSIS",
    `topsis_ahp` = "TOPSIS  vs. AHP",
    `wsm_topsis` = "WSM  vs. TOPSIS",
    `wsm_promethee` = "WSM  vs. PROMETHEE",
    `wsm_ahp` = "WSM  vs. AHP",
    `promethee_ahp` = "PROMETHEE  vs. AHP"))) +
  scale_fill_gradientn(colors = custom_palette2, 
                       name = "Rank Difference", 
                       na.value = "grey50") + 
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank()    )



```




```{r}

st_write(compare, here::here("data", "derived", "compare.geojson"), driver = "GeoJSON")


```


WSM: easy to understand and implement, but assumes that criteria are independent of each other and ranking is highly dependent on the weights

TOPSIS: easy to understand and implement and is more comprehensive than the simple WSM, but it also assumes that criteria are independent of each other, is sensitive to rank reveral issues, and requires additional input froms stakeholders to decide upon the positive and negative ideal scenario. 

PROMETHEE: the most robust statistical mode, but requires careful selection of preference functions and preference function, indifference threshold 

AHP: breaks complex decisions into manageable part, but pairwise comparison can be subjective and biased  and that maintaining consistency in comparison can be challenging. 

# Local MCDA Approach


First, we want to state the difference between global and local MCDA. Both of them are crucial to decision-making in various ways. Global MCDA applies the same criteria weights and decision rules uniformly across the entire study area. It assumes that all areas have the same priorities and conditions for EVSE placement. However, it ignores local variations in socio-economic factors, infrastructure, or mobility behavior.

Local MCDA, on the other hand, recognizes that different parts of the city have different needs, which is typical for decision-making in an urban context. It divides the city into meaningful clusters (based on urban form, population, mobility, etc.). Following this, it conducts MCDA within each cluster using localized weights or decision rules.

Second, why is this important in Philadelphia? Philadelphia is a diverse city with distinct neighborhoods that vary in population density, transit reliance, socioeconomic conditions, and EV adoption. Applying a single set of MCDA weights citywide may not reflect neighborhood-specific needs. For example: Center City & University City: High density, transit-dependent, lower car ownership, meaning population density might be a more important criterion. Northeast Philadelphia: Lower density, higher car ownership, perhaps parking availability might be more important. In suburban areas: EV adoption rates might be more important.With local MCDA, we can ensure context-specific decision-making, which aligns better with urban planning priorities and policies. If we could separate the city into different clusters based on the conditions of the criteria we have, then drawing on our local knowledges, we would be able to assign different weights to different criteria for each cluster. This would allow us to better reflect the local conditions and priorities of each area. 


Third, at what scale should we implement local MCDA? In other words, what do we mean by "local"? This depends on the decision context and the spatial scale of the problem. There are many ways to solve the problem: census tract, neighborhoods, zip codes, etc. In our case, we chose a more systematic way to classify urban areas. We plan to first use PCA to reduce dimensionality by extracting key latent factors from the 9 variables. This also helps to remove redundancy (e.g., population density & driving population may be correlated). Then, we use k-means clustering to group similar grid cells into clusters based on the PCA-transformed data. This ensures that each cluster represents spatial areas with similar urban characteristics. Finally, we perform separate MCDA for each cluster, using different criteria weights. This allows us to capture local variations in EVSE suitability and make more informed decisions. This is also called spectrum clustering, it clusters things better. 

Finally, we would also like to point out that doing so has its particular implications in terms of GIS application, planning, and policymaking. In terms of smart cities planning, having targeted EVSE deployment strategies for high-density vs. suburban areas can help optimize the use of limited resources and ensures that underserved communities receive fair access to EV. Considering global MCDA - it only draws attention to a particular part of the city, which while indeed in need of more EVSE, this might overlook the needs of other neighborhoods. A neighborhood ranked slighly lower in the global MCDA might still need a charging station, or even maybe ranked higher in the local MCDA. As for policymaking, doing so allows us to have cluster-specific incentives for EV adoption or charging subsidies.

## Cluster Determination

```{r pca for 9 criteria}

# Start with fishnet_nowater dataset as the base layer (all criteria aggregated to the fishnet level)
# perform pca analysis
set.seed(123)
pca_fishnet<- prcomp(fishnet_nowater %>% st_drop_geometry %>% select(-uniqueID), 
                     center = TRUE, 
                     scale. = TRUE)

# loadings
pca_loadings <- pca_fishnet$rotation

# pc scores
pc_scores <- as.data.frame(pca_fishnet$x)

# combine dataframe
fishnet_PCA <- cbind(uniqueID = fishnet_nowater$uniqueID, pc_scores)

```

```{r checking}

# unit vector ? 
pc1_norm <- sqrt(sum(pca_loadings[,1]^2))
pc2_norm <- sqrt(sum(pca_loadings[,2]^2))

print(pc1_norm)  # Should be close to 1
print(pc2_norm)  # Should be close to 1

# Compute the dot product
options(digits = 3)
options(scipen = 999)
dot_product <- sum(pca_loadings[,1] * pca_loadings[,2])
print(dot_product) 

# any correlation?
cor(pc_scores[, 1], pc_scores[, 2])

```


```{r pve plot}

pve <- summary(pca_fishnet)$importance[2, ]  # Proportion of variance explained

pve_df <- data.frame(
  PC = seq_along(pve),  # Principal component number
  PVE = pve             # Proportion of variance explained
)

ggplot(pve_df, aes(x = PC, y = PVE)) +
  geom_line(color = "#6EA6C3", size = 1) +   # Blue connecting line
  geom_point(size = 3, color = "#DC6D85") +  # Red points, similar to pch=16 in base R
  theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 15, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8)) +
  labs(title = "PVE Scree Plot of PCA with All 9 Variables",
       x = "Number of PCs",
       y = "Proportion of Variance Explained") +
  scale_x_continuous(breaks = seq(1, length(pve), 1))  

```


```{r}

pve_cumulative <- summary(pca_fishnet)$importance[3, ]  
pve_cumulative

```


```{r biplot}

library(ggbiplot)
ggbiplot(pca_fishnet, 
              choices = c(1, 2),  
              obs.scale = 1, 
              var.scale = 1, 
              groups = rep("All Data", nrow(pca_fishnet$x)),  # Trick to color points
              ellipse = FALSE,   
              circle = FALSE) + 
  scale_color_manual(values = c("#6EA6C3")) +  # Change dot color
  theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 15, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8)) +
  theme(plot.title = element_text(size = 15, face = "bold"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "none") + 
  labs(title = "PCA Biplot", 
       x = "PC1 (First Principal Component)", 
       y = "PC2 (Second Principal Component)")  + 
  geom_hline(yintercept = 0, color = "grey", linewidth = 0.4) + 
  geom_vline(xintercept = 0, color = "grey", linewidth = 0.4)
```


```{r}

# let's use 4 pcs to summarize the overall variability 
# what is an optimal number of clusters? 
k.values <- 1:8
fviz_nbclust(pc_scores[, 1:4], kmeans, method = "wss")

```


```{r}
# let's double check the cluster
# Compute the best number of clusters using 30 indices
nb_result <- NbClust(pc_scores[, 1:4], distance = "euclidean", min.nc = 2, max.nc = 6, method = "kmeans")

# Get the best number of clusters
best_k <- nb_result$Best.nc[1]
print(best_k)
```


```{r run k means}

# since four clusters looks optimal, we will use k = 4
# also we are using 4 PC, which explains 65% of variability 

# run k-means
kmeans_fishnet <- kmeans(pc_scores[, 1:4], centers = 4)

# assign cluster
fishnet_clustered <- data.frame(fishnet_PCA[, 1:4], Cluster = as.factor(kmeans_fishnet$cluster))

# cluster number
table(fishnet_clustered$Cluster)

# combine with original fishnet data
fishnet_clustered <- left_join(fishnet_clustered, fishnet_nowater, by = "uniqueID")
```


```{r}

#  4 clusters with 4 PCs
fishnet_clustered <- st_as_sf(fishnet_clustered)

fishnet_clustered %>% 
  ggplot() +
  geom_sf(aes(fill = Cluster), color = "black") +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "Clustered Fishnet") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )

```

```{r cluster summary}


# figure out common features in a cluster?
# what are the commonalities within the cluster? 

cluster_summary <- fishnet_clustered %>%
  group_by(Cluster) %>%
  dplyr::summarize(
    avg_density = mean(density, na.rm = TRUE),
    avg_driving_pop = mean(driving_pop, na.rm = TRUE),
    avg_evse_nn = mean(evse.nn, na.rm = TRUE),
    avg_TSFCA = mean(TSFCA, na.rm = TRUE),
    avg_elec_cnt = mean(elec_cnt, na.rm = TRUE),
    avg_parking_nn = mean(parking.nn, na.rm = TRUE),
    avg_police_fire = mean(policefire.nn, na.rm = TRUE),
    avg_zoning = mean(zone_score, na.rm = TRUE),
    avg_electric_use = mean(elect_use, na.rm = TRUE),
    count = n()  # Number of grid cells in each cluster
  )

cluster_summary 
```

## Cluster Characteristics and Local Weights

```{r}

ggplot(fishnet_clustered, aes(x = PC1, y = PC2, color = Cluster)) +
  geom_point() +
  labs(title = "Cluster Separation in PCA Space", x = "PC1", y = "PC2") +
  theme_bw() +
  theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 15, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))
```


```{r}

fishnet_clustered %>%
  st_drop_geometry() %>%
  dplyr::select(Cluster, driving_pop, evse.nn, elec_cnt, parking.nn) %>%
  mutate(
    evse.nn = scale_values(evse.nn),
    driving_pop = scale_values(driving_pop),
    elec_cnt = scale_values(elec_cnt),
    parking.nn = scale_values(parking.nn),
  ) %>%
  pivot_longer(cols = -Cluster, names_to = "Variable", values_to = "Scaled_Value") %>%
  ggplot(aes(x = Scaled_Value, fill = as.factor(Cluster))) +
  geom_density(alpha = 0.6) +
  scale_fill_brewer(palette = "Set3") +
  facet_wrap(~Variable, scales = "free") +
  labs(
    title = "Density Distribution of Scaled Criteria by Cluster",
    x = "Scaled Value",
    y = "Density",
    fill = "Cluster"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    legend.position = "bottom",
    strip.text = element_text(size = 11, face = "bold")
  )

```

Now we need to separate the fishnet data into 4 clusters based on the PCA-transformed data. We will then conduct MCDA for each cluster separately. This allows us to capture local variations in EVSE suitability and make more informed decisions. 

```{r separate cluster}

# separate fishnet data into 4 clusters
cluster1 <- fishnet_clustered %>% filter(Cluster == 1)
cluster2 <- fishnet_clustered %>% filter(Cluster == 2)
cluster3 <- fishnet_clustered %>% filter(Cluster == 3)
cluster4 <- fishnet_clustered %>% filter(Cluster == 4)

```

Cluster 1 has the highest average population density, driving population, registered EV, are more accessible to existing EVSE, and have more parking spaces nearby, more commercial space. Grids that fall into this cluster are mostly in Center City and University City.

Cluster 2 has low population density and driving age population, but does have some existing EVSE and electric vehicle registration, most of these places are industrial area that consumes high electricity power. Grids that fall into this cluster are mostly suburbs in Northeast Philadelphia and South Philadelphia, a lot of those are in close proximity to industrial area and the airport. 

Cluster 3 has moderate population density, but more than half of those are in driving age. It is also in the suburb area, but the more affluent area in northwestern Philadelphia like Chestnut Hill where people commute by car. Distance to EVSE is far with low accessibility, but a significant of them have electric vehicle. Also most grids are far from any existing parking lot. Most of those are residential area. Also uses a lot of electricity. 

Cluster 4 also has high population density, and more than half of the population are in driving age. It has a fair number of registered electric vehicle and are mostly residential area. However, accessibility to the existing EVSE is low and mostly are residential area. These are the socio-economically disadvantaged area of Philadelphia, like Kensington, and North Philadelphia. 

That said, from a policymaking aspect, the weight we place on each criterion should be different for each cluster. 

According to the city of Philadelphia and the original conceptualization of the plan, we need to find places to implement new EVSE where there's a large potential demand (densely populated, high driving age population, and high registered EV) but low current accessibility; where there is existing parking location to use, high electricity usage to support the EVSE, and where the neighborhood is safe and suitable for EVSE installation. Among those, prioritize commercial and residential areas. 

Currently, for global MCDA we assigned 0.12 weight to population density, **0.18 to existing EVSE proximity**, **0.15 to registered EVs**, **0.1 to percentage of driving age population**, **0.1 to parking spot proximity**, 0.05 to neighborhood safety score, 0.05 to zoning suitability, 0.12 to EVSE  accessibility, and 0.13 to electricity usage. 

In Cluster 1 (Dense Urban Core)
- high demand for EVSE as justified by its population density, high driving-age population, high registered EVs, and many commercial area
- lower need for new parking spots: More existing infrastructure, so parking spot availability may be less critical.
- Priority: Address accessibility gaps and leverage high demand.


In Cluster 2 (Industrial & High Power Consumption Area)
- low demand for EVSE as justified by its low population density, low driving-age population, and low registered EVs. existing EVSE and EV ownership present but not as high as Cluster 1. 
- high power consumption can support EVSE
- Priority: Focus on leveraging high electricity availability and industrial parking for EVSE deployment.


In Cluster 3 (Affluent Residential Suburb)
- moderate demand for EVSE as justified by its moderate population density, but high driving-age population and many EV owners. 
- low existing EVSE accessibility and public parking spot availability
- Priority: Expand EVSE infrastructure due to high EV ownership but poor accessibility.

In Cluster 4 (Disadvantaged Residential Areas)
- High demand but poor accessibility: High population & driving-age population but few EVSE options.
- Socio-economic considerations: EV adoption is still growing, but accessibility should be prioritized.
- More parking limitations and infrastructure concerns.

We acknowledge that there's subjectivity involved in assigning weights to each criterion. But the goal here is not to be perfect on the weights but to understand how the outcomes can be different when we apply different weights and when we consider local variations. 

We will perform three sets of comparison. 
- Global MCDA (Global Weights) ranks grids based on criteria importance at the citywide scale. This approach ensures consistency across Philadelphia but may overlook local variations.
- Local MCDA (Local Weights) adapts weights to the unique characteristics of each cluster. This method highlights localized priorities that might be underrepresented in a citywide analysis.
- Local MCDA (Global Weights) applies the citywide weight structure within each cluster, allowing us to observe how the same criteria influence different local contexts.


We will use the TOPSIS method to test out the local MCDA model, since TOPSIS was one of the methods used in global MCDA. Running WSM and PROMETHEE locally would introduce more layers of variation, making it harder to isolate the impact of local weight adjustments versus methodological differences. 
Also, TOPSIS avoids WSM's linearity issue. WSM is purely additive, meaning it treats all weighted values as independent and does not consider how alternatives compare to an ideal or worst-case scenario.This can be problematic in localized MCDA where extreme values may exist (e.g., a cluster with very high EV ownership but poor accessibility).
PROMETHEE is a pairwise ranking method, meaning it compares each alternative to every other alternative, which can introduce ranking biases when dealing with smaller, more homogeneous clusters.TOPSIS is computationally less expensive than PROMETHEE, making it easier to run on multiple clusters without excessive processing time.


## Local MCDA Using TOPSIS

```{r construct performance table for local MCDA}

# cluster 1
performanceTable_c1 <- cluster1 %>% 
  dplyr::select(-c(uniqueID, Cluster, PC1, PC2, PC3)) %>% 
  st_drop_geometry() %>% 
  mutate(density = scale_values(density),
         TSFCA = scale_values(TSFCA),
         driving_pop = scale_values(driving_pop),
         evse.nn = scale_values(evse.nn),
         elec_cnt = scale_values(elec_cnt),
         parking.nn = scale_values(parking.nn),
         policefire.nn = scale_values(policefire.nn),
         elect_use = scale_values(elect_use)
         )
uniqueIDs_c1 <- cluster1$uniqueID
rownames(performanceTable_c1) <- uniqueIDs_c1

# cluster 2 
performanceTable_c2 <- cluster2 %>% 
  dplyr::select(-c(uniqueID, Cluster, PC1, PC2, PC3)) %>% 
  st_drop_geometry() %>% 
  mutate(density = scale_values(density),
         TSFCA = scale_values(TSFCA),
         driving_pop = scale_values(driving_pop),
         evse.nn = scale_values(evse.nn),
         elec_cnt = scale_values(elec_cnt),
         parking.nn = scale_values(parking.nn),
         policefire.nn = scale_values(policefire.nn),
         elect_use = scale_values(elect_use)
         )
uniqueIDs_c2 <- cluster2$uniqueID
rownames(performanceTable_c2) <- uniqueIDs_c2


# cluster 3 
performanceTable_c3 <- cluster3 %>% 
  dplyr::select(-c(uniqueID, Cluster, PC1, PC2, PC3)) %>% 
  st_drop_geometry() %>% 
  mutate(density = scale_values(density),
         TSFCA = scale_values(TSFCA),
         driving_pop = scale_values(driving_pop),
         evse.nn = scale_values(evse.nn),
         elec_cnt = scale_values(elec_cnt),
         parking.nn = scale_values(parking.nn),
         policefire.nn = scale_values(policefire.nn),
         elect_use = scale_values(elect_use)
         )
uniqueIDs_c3 <- cluster3$uniqueID
rownames(performanceTable_c3) <- uniqueIDs_c3

# cluster 4
performanceTable_c4 <- cluster4 %>% 
  dplyr::select(-c(uniqueID, Cluster, PC1, PC2, PC3)) %>% 
  st_drop_geometry() %>% 
  mutate(density = scale_values(density),
         TSFCA = scale_values(TSFCA),
         driving_pop = scale_values(driving_pop),
         evse.nn = scale_values(evse.nn),
         elec_cnt = scale_values(elec_cnt),
         parking.nn = scale_values(parking.nn),
         policefire.nn = scale_values(policefire.nn),
         elect_use = scale_values(elect_use)
         )
uniqueIDs_c4 <- cluster4$uniqueID
rownames(performanceTable_c4) <- uniqueIDs_c4


```


```{r local topsis operations}


# we need one set for each cluster
# reverse is min 
# by order, this is evse.nn, density, driving_pop, TSFCA, elec_cnt, parking.nn, policefire.nn, zone score, elect_use
weights_original <- c(0.18, 0.12, 0.10, 0.12, 0.15, 0.10, 0.05, 0.05, 0.13)
weights_c1 <- c(0.12, 0.14, 0.12, 0.16, 0.16, 0.08, 0.06, 0.06, 0.10)
weights_c2 <- c(0.14, 0.08, 0.08, 0.10, 0.12, 0.12, 0.06, 0.10, 0.20)
weights_c3 <- c(0.20, 0.10, 0.12, 0.14, 0.18, 0.06, 0.06, 0.06, 0.08)
weights_c4 <- c(0.16, 0.14, 0.12, 0.18, 0.12, 0.08, 0.04, 0.06, 0.10)

names(weights_c1) <- colnames(performanceTable_c1)
names(criteriaMinMax) <- colnames(performanceTable_c1)
names(weights_c4) <- colnames(performanceTable_c4)
names(criteriaMinMax) <- colnames(performanceTable_c4)

# check total 
sum(weights_c1)  # Should be 1
sum(weights_c2)  # Should be 1
sum(weights_c3)  # Should be 1
sum(weights_c4)  # Should be 1

criteriaMinMax <- c("min", "max", "max", "min", "max", "min", "min", "max", "max")

## cluster 1
overall_c1 <- TOPSIS(performanceTable_c1, weights_c1, criteriaMinMax)
TOPSIS_c1 <- data.frame(uniqueID = names(overall_c1), values = overall_c1) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c1 = min_rank(desc(values)),
         topsis_percentile_c1 = percent_rank(topsis_rank_c1) * 100)

# test for original
overall_c1_org <- TOPSIS(performanceTable_c1, weights_original, criteriaMinMax)
TOPSIS_c1_org <- data.frame(uniqueID = names(overall_c1_org), values = overall_c1_org) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c1_org = min_rank(desc(values)),
         topsis_percentile_c1_org = percent_rank(topsis_rank_c1_org) * 100)

## cluster 2
overall_c2 <- TOPSIS(performanceTable_c2, weights_c2, criteriaMinMax)
TOPSIS_c2 <- data.frame(uniqueID = names(overall_c2), values = overall_c2) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c2 = min_rank(desc(values)),
         topsis_percentile_c2 = percent_rank(topsis_rank_c2) * 100)

# test for original
overall_c2_org <- TOPSIS(performanceTable_c2, weights_original, criteriaMinMax)
TOPSIS_c2_org <- data.frame(uniqueID = names(overall_c2_org), values = overall_c2_org) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c2_org = min_rank(desc(values)),
         topsis_percentile_c2_org = percent_rank(topsis_rank_c2_org) * 100)

## cluster 3
overall_c3 <- TOPSIS(performanceTable_c3, weights_c3, criteriaMinMax)
TOPSIS_c3 <- data.frame(uniqueID = names(overall_c3), values = overall_c3) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c3 = min_rank(desc(values)),
         topsis_percentile_c3 = percent_rank(topsis_rank_c3) * 100)

# test for original
overall_c3_org <- TOPSIS(performanceTable_c3, weights_original, criteriaMinMax)
TOPSIS_c3_org <- data.frame(uniqueID = names(overall_c3_org), values = overall_c3_org) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c3_org = min_rank(desc(values)),
         topsis_percentile_c3_org = percent_rank(topsis_rank_c3_org) * 100)

## cluster 4
overall_c4 <- TOPSIS(performanceTable_c4, weights_c4, criteriaMinMax)
TOPSIS_c4 <- data.frame(uniqueID = names(overall_c4), values = overall_c4) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c4 = min_rank(desc(values)),
         topsis_percentile_c4 = percent_rank(topsis_rank_c4) * 100)

# test for original
overall_c4_org <- TOPSIS(performanceTable_c4, weights_original, criteriaMinMax)
TOPSIS_c4_org <- data.frame(uniqueID = names(overall_c4_org), values = overall_c4_org) %>%  
  mutate(uniqueID = as.integer(uniqueID)) %>% 
  mutate(topsis_rank_c4_org = min_rank(desc(values)),
         topsis_percentile_c4_org = percent_rank(topsis_rank_c4_org) * 100)

```


might not need this 
```{r original percentile}

# Get the original rank and compute percentile rank 
# then join it with the result from local MCDA
# to compare whether the percentile rank is higher or lower than the global MCDA
final_percentile <- final %>% 
  dplyr::select(NAME, uniqueID, avg_rank, topsis_rank, wsm_rank) %>%
  mutate(
    avg_percentile = percent_rank(avg_rank) * 100,
    topsis_percentile = percent_rank(topsis_rank) * 100)
  
```

might not need this as well
```{r compare global and local }

# quick check for cluster 1
TOPSIS_c1_cp  <- TOPSIS_c1 %>% 
  left_join(TOPSIS_c1_org %>% dplyr::select(uniqueID, topsis_percentile_c1_org), by = "uniqueID") %>%
  left_join(final_percentile %>% dplyr::select(uniqueID, 
                                               NAME, 
                                               topsis_percentile, 
                                               avg_percentile), by = "uniqueID") %>%
  dplyr::select(NAME, uniqueID, avg_percentile, topsis_percentile, topsis_percentile_c1, topsis_percentile_c1_org)

# quick check for cluster 2
TOPSIS_c2_cp  <- TOPSIS_c2 %>% 
  left_join(TOPSIS_c2_org %>% dplyr::select(uniqueID, topsis_rank_c1_org, topsis_percentile_c2_org), by = "uniqueID") %>%
  left_join(final_percentile %>% dplyr::select(uniqueID, 
                                               NAME, 
                                               topsis_percentile, 
                                               avg_percentile
                                          ), by = "uniqueID") %>%
  dplyr::select(NAME, uniqueID, avg_percentile, topsis_percentile, topsis_percentile_c2, topsis_percentile_c2_org)

# quick check for cluster 3
TOPSIS_c3_cp  <- TOPSIS_c3 %>% 
  left_join(TOPSIS_c3_org %>% dplyr::select(uniqueID, topsis_percentile_c3_org), by = "uniqueID") %>%
  left_join(final_percentile %>% dplyr::select(uniqueID, 
                                               NAME, 
                                               topsis_percentile, 
                                               avg_percentile
                                          ), by = "uniqueID") %>%
  dplyr::select(NAME, uniqueID, avg_percentile, topsis_percentile, topsis_percentile_c3, topsis_percentile_c3_org)

# quick check for cluster 4
TOPSIS_c4_cp  <- TOPSIS_c4 %>% 
  left_join(TOPSIS_c4_org %>% dplyr::select(uniqueID, topsis_percentile_c4_org), by = "uniqueID") %>%
  left_join(final_percentile %>% dplyr::select(uniqueID, 
                                               NAME, 
                                               topsis_percentile, 
                                               avg_percentile
                                          ), by = "uniqueID") %>%
  dplyr::select(NAME, uniqueID, avg_percentile, topsis_percentile, topsis_percentile_c4, topsis_percentile_c4_org)

```


```{r global and local score combined }

global_local <- TOPSIS %>% 
  mutate(global_val = values) %>% 
  mutate(global_rank = topsis_rank) %>% 
  left_join(TOPSIS_c1 %>% 
              mutate(c1_val = values, cluster_1 = "1") %>% 
              dplyr::select(uniqueID, c1_val, cluster_1), by = "uniqueID") %>%
  left_join(TOPSIS_c1_org %>% 
              mutate(c1_val_org = values) %>% 
              dplyr::select(uniqueID, c1_val_org), by = "uniqueID") %>%
  left_join(TOPSIS_c2 %>% 
              mutate(c2_val = values, cluster_2 = "2") %>% 
              dplyr::select(uniqueID, c2_val, cluster_2), by = "uniqueID") %>%
  left_join(TOPSIS_c2_org %>%
              mutate(c2_val_org = values) %>% 
              dplyr::select(uniqueID, c2_val_org), by = "uniqueID") %>%
  left_join(TOPSIS_c3 %>% 
              mutate(c3_val = values, cluster_3 = "3") %>% 
              dplyr::select(uniqueID, c3_val, cluster_3), by = "uniqueID") %>%
  left_join(TOPSIS_c3_org %>%
              mutate(c3_val_org = values) %>% 
              dplyr::select(uniqueID, c3_val_org), by = "uniqueID") %>%
  left_join(TOPSIS_c4 %>% 
              mutate(c4_val = values, cluster_4 = "4") %>% 
              dplyr::select(uniqueID, c4_val, cluster_4), by = "uniqueID") %>%
  left_join(TOPSIS_c4_org %>%
              mutate(c4_val_org = values) %>% 
              dplyr::select(uniqueID, c4_val_org), by = "uniqueID") %>%
  mutate(local_val = coalesce(c1_val, c2_val, c3_val, c4_val)) %>%
  mutate(cluster = coalesce(cluster_1, cluster_2, cluster_3, cluster_4)) %>%
  mutate(local_rank = min_rank(desc(local_val))) %>%
  mutate(local_val_org = coalesce(c1_val_org, c2_val_org, c3_val_org, c4_val_org)) %>%
  mutate(local_rank_org = min_rank(desc(local_val_org))) %>%
  dplyr::select(uniqueID, global_val, global_rank, local_val, local_rank, local_val_org, local_rank_org, cluster)
  

```


## Mapping Global and Local MCDA Results

```{r fig.height=8, fig.width=12}

global_local_toplot <- fishnet_nowater %>% 
  dplyr::select(uniqueID, geometry) %>%
  left_join(global_local, by = "uniqueID")

quantiles_glo <- classIntervals(global_local_toplot$global_val, n = 5, style = "quantile")
quantiles_loc <- classIntervals(global_local_toplot$local_val, n = 5, style = "quantile")
quantiles_loc_org <- classIntervals(global_local_toplot$local_val_org, n = 5, style = "quantile")

global_local_toplot$quantile_glo <- cut(global_local_toplot$global_val, breaks = quantiles_glo$brks, include.lowest = TRUE)
global_local_toplot$quantile_loc <- cut(global_local_toplot$local_val, breaks = quantiles_loc$brks, include.lowest = TRUE)
global_local_toplot$quantile_loc_org <- cut(global_local_toplot$local_val_org, breaks = quantiles_loc_org$brks, include.lowest = TRUE)

glo <- ggplot()+
 geom_sf(data=fishnet, color=NA, fill="white") +
 geom_sf(data = global_local_toplot, aes(fill = quantile_glo), color = NA) +
  scale_fill_manual(values = custom_palette) +
  labs(title = "Global TOPSIS Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )

loc <- ggplot()+
 geom_sf(data=fishnet, color=NA, fill="white") +
 geom_sf(data = global_local_toplot, aes(fill = quantile_loc), color = NA) +
  scale_fill_manual(values = custom_palette) +
  labs(title = "Local TOPSIS Score", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )


loc_org <- ggplot()+
 geom_sf(data=fishnet, color=NA, fill="white") +
 geom_sf(data = global_local_toplot, aes(fill = quantile_loc_org), color = NA) +
  scale_fill_manual(values = custom_palette) +
  labs(title = "Local TOPSIS Score Original", fill = "Bins") +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, size=0.8)
        )


ggarrange(glo, loc, loc_org, ncol = 2, nrow = 2)

```

## Compute Spatial Statistics

```{r}

library(spdep)

# global moran's I for global MCDA
coords <-  st_coordinates(global_local_toplot %>% st_centroid()) 
neighborList <- knn2nb(knearneigh(coords, 8))
spatialWeights <- nb2listw(neighborList, style="W")

moranTest_glo <- moran.mc(global_local_toplot$global_val, 
                      spatialWeights, nsim = 999)

moranTest_loc <- moran.mc(global_local_toplot$local_val, 
                      spatialWeights, nsim = 999)

expected_index <- mean(moranTest_glo$res)
variance <- sum((moranTest_glo$res - expected_index)^2) / (length(moranTest_glo$res) - 1)
z_score <- (moranTest_glo$statistic - expected_index) / sqrt(variance)

expected_index_loc <- mean(moranTest_loc$res)
variance_loc <- sum((moranTest_loc$res - expected_index_loc)^2) / (length(moranTest_loc$res) - 1)
z_score_loc <- (moranTest_loc$statistic - expected_index_loc) / sqrt(variance_loc)

glo_moran_summary <- data.frame(Statistic = c("Moran's Index", "Expected Index", "Variance", "z-score", "p-value"),
                                Value = c(moranTest_glo$statistic, expected_index, variance, z_score, moranTest_glo$p.value))

loc_moran_summary <- data.frame(
  Statistic = c("Moran's Index", "Expected Index", "Variance", "z-score", "p-value"),
  Value = c(moranTest_loc$statistic, expected_index_loc, variance_loc, z_score_loc, moranTest_loc$p.value)
)


ggplot(as.data.frame(moranTest_loc$res[c(1:999)]), aes(moranTest_loc$res[c(1:999)])) +
  geom_histogram(binwidth = 0.01, fill = "black") +
  geom_vline(aes(xintercept = moranTest_glo$statistic), colour = "#6E0E0A" ,lwd=1) +
  labs(title = "Observed and Permuted Moran's I",
       subtitle = "Observed Moran's I in Red",
       x = "Moran's I",
       y = "Count") +
  theme_light() +   
theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=10))

glo_moran_summary
loc_moran_summary
```


```{r}

library(sfdep)
library(spdep)
# Local Getis Ord 
neigh_nbs <- global_local_toplot %>% 
  mutate(
    nb = st_knn(geometry, k = 8),  
    wt = st_weights(nb),              
    neigh_lag = st_lag(global_val, nb, wt)  
  )


gi_hot_spots <- neigh_nbs %>% 
  mutate(Gi = local_g_perm(global_val, nb, wt, nsim = 999)) %>% 
  unnest(Gi)


gi_hot_spots <- gi_hot_spots %>%  
  dplyr::select(gi, p_folded_sim, uniqueID) |> 
  mutate(
    classification = case_when(
      # Classify based on the following criteria:
      gi > 0 & p_folded_sim <= 0.01 ~ "Very hot",
      gi > 0 & p_folded_sim <= 0.05 ~ "Hot",
      gi > 0 & p_folded_sim <= 0.1 ~ "Somewhat hot",
      gi < 0 & p_folded_sim <= 0.01 ~ "Very cold",
      gi < 0 & p_folded_sim <= 0.05 ~ "Cold",
      gi < 0 & p_folded_sim <= 0.1 ~ "Somewhat cold",
      TRUE ~ "Insignificant"
    ),    # Convert 'classification' into a factor for easier plotting
    classification = factor(
      classification,
      levels = c("Very hot", "Hot", "Somewhat hot",
                 "Insignificant",
                 "Somewhat cold", "Cold", "Very cold")
    )
  )

gi_hot_spots %>% 
  ggplot() + 
  geom_sf(aes(fill = classification), color = "black", lwd = 0.1) +
  scale_fill_brewer(type = "div", palette = 6, name = "Classification") +
  labs(title = "Spatial Hotspots and Coldspots") +
  theme(legend.position="bottom",
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks =element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth=0.8)
        )

st_write(gi_hot_spots, "/Users/emzhou/Desktop/Paper/svg/globalhotspot.geojson", driver = "GeoJSON")
```


## Compute Distribution Statistics 

```{r}

# visualize the distribution of global and local MCDA scores
ggplot(global_local, aes(x = global_val)) +
  geom_histogram(binwidth = 0.01, fill = "black") +
  geom_vline(aes(xintercept = mean(global_val)), colour = "#6E0E0A" ,lwd=1) +
  labs(title = "Global MCDA Score Distribution",
       subtitle = "Mean Score in Red",
       x = "MCDA Score",
       y = "Count") +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=10))

# visualize the distribution of global and local MCDA scores
ggplot(global_local, aes(x = local_val)) +
  geom_histogram(binwidth = 0.01, aes(fill = cluster)) +
  geom_vline(aes(xintercept = mean(global_val)), colour = "#6E0E0A" ,lwd=1) +
  labs(title = "Global MCDA Score Distribution",
       subtitle = "Mean Score in Red",
       x = "MCDA Score",
       y = "Count") +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=10))
```

```{r}
summary(TOPSIS_c1$values)
summary(TOPSIS_c2$values)
summary(TOPSIS_c3$values)
summary(TOPSIS_c4$values)
```


This highlights some fundamental difference between how scores are distributed when moving from a single unified ranking (global) to localized evaluations.

- Local MCDA scores are relative within each cluster.If Cluster 1 consists of grids that are all highly ranked globally, but they are similar to each other, then the range of local scores will be compressed.This could lead to Cluster 1 having the lowest local MCDA scores on average because the best grids in that cluster do not stand out as much relative to their peers. Cluster 1 has a tight range of values (0.160 to 0.592), which suggests that grids in Cluster 1 are more homogeneous and do not have strong local contrasts.

- In Cluster 3 and Cluster 4, the range of scores is much wider. This suggests that there’s greater differentiation among grids in these clusters, which can lead to higher local scores when some grids outperform others significantly.

- This result suggests that while Cluster 1 is ranked highest in the global MCDA, it may be internally more uniform, meaning that within its local context, no grid significantly outperforms others. In a global setting, Cluster 1's grids might dominate the rankings. In a local setting, they are competing only within their cluster, where the differences between the best and worst-performing grids are small.Conversely, clusters with more variation internally (e.g., Cluster 3 and 4) will see bigger score differences in local MCDA. This supports the idea that local MCDA captures within-cluster variation better than global MCDA, while global MCDA favors clusters that perform well across all criteria.

Homogeneity Compresses Scores. All grids in cluster 1 have good access to EVSE. Whereas in cluster 3 and 4, this suggests more variations and diverse cluster. If Cluster 1 had the best infrastructure or services in the global MCDA, it dominated the global rankings.

Local MCDA is sensitive to intra-cluster variance.If Cluster 3 and Cluster 4 have a mix of very strong and weak grids, then some grids will rise to the top.

## Furhter Comparisons Within Cluster

```{r}

# Cluster 1
cor(global_local %>% filter(cluster == "1") %>% dplyr::select(global_val, local_val, local_val_org), method = "spearman")

# Cluster 2
cor(global_local %>% filter(cluster == "2") %>% dplyr::select(global_val, local_val, local_val_org), method = "spearman")

# Cluster 3
cor(global_local %>% filter(cluster == "3") %>% dplyr::select(global_val, local_val, local_val_org), method = "spearman")

# Cluster 4
cor(global_local %>% filter(cluster == "4") %>% dplyr::select(global_val, local_val, local_val_org), method = "spearman")

```
```{r}
print(results)

```


This analysis explores the relationship between local and global TOPSIS approaches across four clusters, evaluating how different weighting strategies impact decision-making outcomes. Correlation coefficients (r) and R-squared (R²) values provide insights into the alignment between methods and the extent to which local preferences (either local or global weights) deviate from global decision-making frameworks.

In Cluster 1: 
The strong correlation between Local TOPSIS (Global Weights) and Local TOPSIS (Local Weights) (r = 0.901, R² = 0.811) suggests that local alternatives are relatively stable when weights are shifted between local and global perspectives. This means local rankings do not change dramatically when weights are aggregated at a higher level. 
However, the weak relationship between Local TOPSIS (Global Weights) and Global TOPSIS (r = 0.207, R² = 0.106) suggests that local preferences differ significantly from the global ranking. This could indicate that the criteria importance differs significantly between the local and global perspectives.
The moderate agreement between Local TOPSIS (Local Weights) and Global TOPSIS (r = 0.870, R² = 0.305) suggests that while some commonalities exist, the local weighting approach captures unique regional priorities not reflected in the global ranking.

> A global model might not fully capture the localized nuances that influence the ranking of alternatives, and using global weights at a local level could introduce bias or misalignment.

In Cluster 2:
High correlation across all comparisons (r > 0.8) suggests that rankings remain largely stable across local and global weighting strategies.
The strong agreement between Local TOPSIS (Local Weights) and Global TOPSIS (r = 0.979, R² = 0.868) indicates that local priorities are well-aligned with the global ranking, suggesting that the criteria importance at the global level is well-reflected in this region.
The high agreement between Local TOPSIS (Global Weights) and Global TOPSIS (r = 0.802, R² = 0.711) indicates that using global weights in local analyses still produces relatively consistent results, meaning that global prioritization does not significantly distort local rankings.

In Cluster 3:
Strong agreement between all comparisons suggests that rankings are relatively consistent across weighting strategies, though not as perfectly aligned as Cluster 2.
The high R² values (> 0.7) indicate that local variations in weighting still have some impact on rankings, but not to a degree that would create large discrepancies between global and local perspectives.

In Cluster 4: 
Local TOPSIS (Global Weights) and Local TOPSIS (Local Weights) are highly correlated (r = 0.962, R² = 0.925), meaning that local rankings remain stable regardless of weighting choice.
The weak correlation between Local TOPSIS (Global Weights) and Global TOPSIS (r = 0.491, R² = 0.406) suggests that local rankings are not well-aligned with the global approach, meaning that local conditions heavily influence decision-making outcomes.
The moderate agreement between Local TOPSIS (Local Weights) and Global TOPSIS (r = 0.837, R² = 0.416) indicates that while some overlap exists, global rankings alone may not be reliable for decision-making in this cluster.

> In this case, a strictly global approach could misrepresent local priorities. Local weightings should be prioritized in decision models to ensure that rankings reflect actual on-the-ground conditions.


See an example here: 
842 Frankford. in Cluster 4
Ranked 529 in Global
Ranked 24 in Local
This is located at the intersection of Tacony Creek Park and Junita Park, close to Adams Ave. There is a shopping mall with big parking lots. 

Average global rank in cluster 1 is 228
Average global rank in cluster 2 is 666
Average global rank in cluster 3 is 744
Average global rank in cluster 4 is 694

In Cluster 1, 165 out of 246 grids were in the top 20% in global model. 
In Cluster 2, 25 out of 190 grids were in the top 20% in global model.
In Cluster 3, 28 out of 269 grids were in the top 20% in global model. 
In Cluster 4, 25 out of 507 grids were in the top 20% in global model. 

In Cluster 4, there were 219 out of 507 grids ranked not in the top 20% in the global model but were in the top 20% in the local model. 

Using local MCDA, we will mostly capture the original top20% in cluster 1. But more importantly, a lot of areas in cluster 4 that were originally neglected for their low rank will now be given more attention. 

```{r}

global_local%>%
  mutate(global_percentile = percent_rank(global_rank),
         local_percentile = percent_rank(local_rank),
         local_org_percentile = percent_rank(local_rank_org)) %>%
  mutate(
    global_cat = cut(global_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")),
    local_cat = cut(local_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%")),
    local_org_cat = cut(local_org_percentile, breaks = c(-1, 0.2, 0.4, 0.6, 0.8, 1.1), labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%"))
  ) %>% 
  filter(cluster == "1") %>%
  filter(global_cat == "0-20%")
  # filter(global_cat == "0-20%" & local_cat != "0-20%")
```

When running global MCDA, every grid competes against all grids. When running Cluster specific MCDA, grids only compete within their cluster, causing a redistribution of relative rankings A grid that was middle-rank globally may appear much worse locally if it no longer competes with weaker grids outside its cluster. In other words, in the global model, these ideal/worst solutions are calculated from all grids, meaning that certain grids may appear closer to the ideal simply because they are better than the worst performers globally. Also in the global model, percentile rankings are spread across all grids, whereas in Cluster 1, only a subset of grids exist.If a grid was ranked 60% globally but now competes against only the densest, highest-priority locations, it might drop into a much lower percentile


Rank reversal highlights the influence of localized criteria adjustment, where grids that were ranked high globally are now low-ranked within Cluster 1 demonstrates how local prioritization can shift the relative importance of criteria. This is an expected outcome of context-sensitive decision-making as the relative distances to ideal and anti-ideal solutions shift when the evaluation is confined to specific clusters. A good example could be that globally, a grid may rank well due to high EV registration and population density.Locally, in a cluster where most grids already have high EV registration and we put more priority on other criteria, this same grid no longer stands out, leading to a lower ranking.

That said, local MCDA helps to avoid Over-Investment in already well-served areas. If only global MCDA were used, policymakers might favor areas that already rank high globally, potentially reinforcing existing disparities.Local MCDA allows a more targeted, equity-focused approach, ensuring that EVSE expansion reaches areas with high demand but lower current accessibility.For example, in cluster 4. global MCDA might rank some grids lower due to low current EV ownership, even though they urgently need better accessibility. Local MCDA increases the weight of EVSE accessibility and population density, ensuring these areas are prioritized for investment.

Therefore, we should take a multi-scale approach (global + local) . This ensures that investments in EVSE infrastructure align with both system-wide efficiency and local equity goals. Global MCDA determines which neighborhoods or districts should receive EVSE investment. Local MCDA ensures that within those neighborhoods, resources are allocated to grids that need it the most, rather than just reinforcing existing infrastructure.


# Public Private Partnership

We decide to find places in Philadelphia that is consistently ranked as more suitable for new EVSE among all MCDA models. Through consulting the output of our  site visits, various stakeholder meetings, and budget considerations, we selected one of those sites in South Philadelphia and proposed a public-private partnership model between Philadelphia’s OIT and local grocery stores to install and maintain the EVSEs. We subsequently conducted financial analyses for the cost and revenue of breakdowns and designed a phased implementation of EVSE infrastructure given the current site conditions. 

```{r load acme locations}

acme <- read.csv(here::here("data", "raw", "acme.csv")) %>% 
  st_as_sf(., coords = c("Longitude", "Latitude"), crs = 4326) %>% 
  mutate(longitude= unlist(map(geometry, 1)),
         latitude = unlist(map(geometry, 2)))

```


```{r leaflet plots}

pal <- colorBin(palette = custom_palette, # change the color palette here
                domain = final$avg_rank, # change the parameter name after $
                bins = c(0, 200, 400, 600, 800, 1000, 1200)) # change the bins here 


labs <- paste0("<strong>", final$NAME,
               "</strong><br/>ID: ", final$uniqueID, 
               "</strong><br/>Avg Rank: ", final$avg_rank, 
               "</strong><br/>WSM Rank: ", final$wsm_rank, 
               "</strong><br/>TOPSIS Rank: ", final$topsis_rank, 
               "</strong><br/>PROMETHEE Rank: ", final$promethee_rank, 
               "</strong><br/>AHP-TOPSIS Rank: ", final$ahptopsis_rank
               ) 

labs_acme <- paste0("</strong><br/>ACME Address ", acme$Address)
acme_marker_color <- "black"
```

```{r}
leaflet <- final %>%
  st_transform("EPSG:4326") %>% 
  leaflet() %>%
  setView(lng = -75.1652, lat = 40.0126, zoom = 10.5) %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addPolygons(weight = 0.2,
              color = "black", 
              opacity = 0.5,
              fillColor = ~pal(avg_rank), # change the name of the parameter in the bracket
              popup = labs,
              fillOpacity = 2) %>%
  addCircleMarkers(
    data = acme,  # Replace with your actual point data
    lng = ~longitude,  # Replace with the appropriate variable in your point dataset for longitude
    lat = ~latitude,  # Replace with the appropriate variable in your point dataset for latitude
    weight = 0.5,
    radius = 4,
    fillColor = "black",  # Set the fill color for the markers
    fillOpacity = 1,
    popup = ~labs_acme  # Replace with the appropriate variable in your point dataset for popup text
  ) %>% 
  addLegend(position = "bottomright",
            pal = pal,
            values = ~avg_rank, # change the name of the parameter after ~
            title = "EVSE Suitability Rank") %>% 
  addControl(
    html = paste0('<div style="background: white; padding: 2px; border-radius: 2px;">
                    <i style="background:', acme_marker_color, '; width: 15px; height: 15px; display: inline-block; border-radius: 50%;"></i> ACME<br>
                    </div>'),
    position = "bottomright"
  )

leaflet
```



# Takeaways

GIS-based MCDA is a robust criteria based methodology that support multiple criteria and statistical models at once, which allows for more in-depth decision making in the planning field. 

Our study demonstrate a spatially informed starting point to identify potential areas for new EVSE that contribute to Philadelphia's sustainable transportation goals. 

However, our study also highlighted several challenges in using MCDA. This includes agreeing on the input criteria, the weighting schemes, and various other inputs required for MCDA models. Considering the number of goals for this particular decision, equity,sustainable transportation, affordability, etc.,  what should be included and how important they are become important questions to consider and would significantly alter the result. 

In addition, we would like to point out that there's no perfect decision making model. Methods that are more comprehensive and robust mathematically 1) requires more decision inputs, which is more time consuming, introduces more inconsistencies from stakeholders, and introduces more subjectivity 2) could be more computationally intensive, 3) less intuitive to non-experts. The best practice is to not rely on a single method in decision making. Beyond using MCDA models, it is essential to check the actual site condition before making final decisions. Note that there are things that spatial analysis cannot capture.

Moreover, our study have also highlighted several recurrent challenges in geospatial model for decision making. Specifically, the modifiable areal unit problem (MAUP) is a source of uncertainty, considering that out raw data all comes with different spatial unit and needs to be re-aggregated before proceed. In addition, ecological fallacy could be an issues as we are making inference about a small portion of the population using tract level statistics for the whole population.Moreover, errors and uncertainties may also arise when dealing with missing information, when we make assumptions/spatial interpolations about a neighboring grids' situation.