---
title: "MTA Use Case"
author: "Ronan Ysebaert"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    css: styles.css
    toc: yes
    toc_depth: 2
css: styles.css
bibliography: mta.bib
link-citations: yes
vignette: >
  %\VignetteIndexEntry{MTA Use Case} 
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown} 
editor_options:   
  chunk_output_type: console
---

<br>

This vignette describes a typical use case of the `MTA` package with a concrete example: the analysis of income inequalities in the Greater Paris Metropolitan Area. This example plays all the functions of the MTA package with appropriate visualizations (plots, maps).

<br>

# 1. Introducing case-study and MTA

## 1.1 Multiscalar Territorial Analysis for Policy Study

Beginning in the MTA conceptual framework, several issues must be considered: the study area, the territorial hierarchy and the selected indicator.

* **Metropolis of Greater Paris (_MGP_) : a study area making sense in a political context**.    
The _MGP_ has been created the 1st January 2016 after a preparatory mission and aims at finding solutions to improve and develop the potential of the Paris metropolitan area to be more present at international level (following the initiatives that took place in London, New-York or Tokyo). The [official website](https://www.metropolegrandparis.fr/fr "Métropole du Grand Paris - Mission de préfiguration") (last visit: 2022-06) of this infrastructure provides all the details of this new level of governance. Namely, this spatial planning project proposes to _"develop a stronger solidarity between territories, reduce territorial inequalities and propose a rebalancing in term of access to household, employment, training, services and amenities"_. For this emerging zoning of spatial planning, the need for indicators and analysis of territorial inequalities is high. The MTA methodology seems to be particularly adapted to propose territorial evidences in this context.

* **From the IRIS to the _MGP_, a lot of hierarchical territorial divisions are available**.   
The lowest territorial level existing in France associated with public statistical information is _IRIS_ _(Ilots Regroupés pour l'Information Statistique)_. The reference size targeted for each IRIS is [2000 inhabitants](https://www.insee.fr/fr/metadonnees/definition/c1523 "IRIS definition, INSEE") (last visit: 2022-06). Each municipality over 10 000 inhabitants and a high proportion of municipalities between 5 000 and 10 000 inhabitants are split in IRIS. It is possible to merge these IRIS in _municipalities_ (36 000 territorial units in France). These municipalities can also be aggregated in _Établissements Publics Territoriaux_(equivalent to _Établissement de Coopération intercommunale, EPCI)_ or in _départements_. Each territorial context has specific competences in spatial planning term. Thus, the _MGP_ is constituted by a large number of territorial divisions. It allows the calculation of a high number of deviations according to these territorial contexts. **This analysis is mainly focused on *municipalities* and *Établissements Publics Territoriaux (_EPT_)* ** but could be extended to other territorial levels. 

* **Indicators relevant at this scale of analysis**.  
The chosen indicator (average income tax reference per households) is one of the key drivers of spatial planning policies at local level in France. It is one of the indicators considered relevant for city policy ([_Politique de la Ville_](https://www.insee.fr/fr/statistiques/2500477 "Key indicators of the French urban policy"), last visit: 2022-06).  

This vignette proposes some concrete outputs for improving the knowledge on income inequalities and proposing solutions for a better geographical distribution of wealth in the _MGP_ area. Relevant statistics will be computed using MTA functionalities and plotted on graphics and maps.  


## 1.2 Dataset and complementary Packages

The example dataset, `GrandParisMetropole` includes all the required information to create 

* com is a sf object. It includes the geometries, identifiers (DEPCOM) and names (LIBCOM) of the 150 _MGP_ municipalities
 identifiers and names of the _EPT_ of belonging of each municipality (EPT and LIB_EPT), identifier of the *département* of belonging of each municipality (DEP) and the statistical indicators for the multiscalar inequality analysis, the numerator (amount of income tax references in Euros, INC) and the denominator (number of tax households, TH) (see ?com).

* EPT is a sf object. It includes the geometries of the 12 _EPT_ of the _MGP_ (see ?EPT). 

* cardist is a matrix including the a time distance matrix between municipalities computed using `osrm` package (see ?cardist). 

Two additional packages are used: [`mapsf`](https://riatelab.github.io/mapsf/) for thematic mapping purposes and `ineq` for computing inequality indexes and Lorenz Curve Plot.

````{r input, warning = FALSE, cache = FALSE}
# load packages
library(sf)
library(MTA)
library(mapsf)
library(ineq)

# load dataset
com <- st_read(system.file("metroparis.gpkg", package = "MTA"),
               layer = "com", quiet = TRUE)
ept <- st_read(system.file("metroparis.gpkg", package = "MTA"), 
               layer = "ept", quiet = TRUE)

# set row names to municipalities names
row.names(com) <- com$LIBCOM
````

## 1.3 Context Maps

Context maps highlight the territorial organization of the study area and provide some first insights regarding the spatial patterns introduced by the indicator used for the analysis. 

### 1.3.1 Study Area

The 150 municipalities of the _MGP_ belongs to 12 intermediate territorial divisions: the 
_Établissements Publics Territoriaux_ (EPT). This territorial zoning follows approximately the delineation of the _départements_ : _Seine-Saint-Denis_ (in blue on the map below); Paris (in purple); _Hauts-de-Seine_ (in green) and _Val-de-Marine (in orange). 

````{r plot_zonings, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.align = 'center'}
# label / colors management
LIBEPT <- c("Paris", "Est Ensemble", "Grand-Paris Est", 
            "Territoire des aeroports", "Plaine Commune", "Boucle Nord 92", 
            "La Defense","Grand Paris Sud Ouest", "Sud Hauts-de-Seine", 
            "Val de Bievres - Seine Amond - Grand Orly",
            "Plaine Centrale - Haut Val-de-Marne - Plateau Briard",
            "Association des Communes de l'Est Parisien")

# colors
cols <- c("#cfcfcf", # Grey(Paris)
          "#92C8E0", "#7BB6D3", "#64A4C5", "#458DB3", # Blues (dept 93)
          "#A6CC99", "#8CBB80", "#71A966", "#4E9345", # Greens (dept 92)
          "#F38F84", "#EF6860", "#EA3531") # Reds (dept 94)

colEpt <- data.frame(LIBEPT, cols)

# zoning
mf_map(x = com, var="LIBEPT", type = "typo",
       val_order = LIBEPT, pal = cols, lwd = 0.2, border = "white",
       leg_pos = "left", leg_title = "EPT of belonging")

mf_map(ept, col = NA, border = "black", add = TRUE)

# layout 
mf_layout(title = "Territorial Zoning of the MGP",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"))
````




### 1.3.2 Map numerator and denominator

The numerator (number of tax households) and the denominator (total amount of income tax references) building the reference ratio (average income per tax households) are firstly plotted.

````{r INCDeINC_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.show="hold"}
# layout
opar <- par(mfrow = c(1, 2))

# numerator map
com$INCM <- com$INC / 1000000
mf_map(com,  col = "peachpuff",  border = NA)
mf_map(ept,  col = NA,  border = "black", add = TRUE)

mf_map(x = com, var = "INCM", type = "prop",
       col =  "#F6533A", inches = 0.15, border = "white",
       leg_pos = "topleft", leg_val_rnd =  0,
       leg_title = "Amount of income taxe reference\n(millions of euros)")

# layout
mf_layout(title = "Numerator - Amount of income tax reference",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)


# denominator map
mf_map(com,  col = "peachpuff",  border = NA)
mf_map(ept,  col = NA,  border = "black", add = TRUE)

mf_map(x = com, var = "TH", type = "prop",
       col =  "#515FAA", inches = 0.15, border = "white",
       leg_pos = "topleft", leg_val_rnd =  -2,
       leg_title = "Number of tax households", add = TRUE)

# layout
mf_layout(title = "Denominator - Tax households",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
par(opar)
````

Without surprise, the highest amounts of tax households and income are located 
in the central area of the _MGP_ (Parisian _Arrondissements_). That being said, these 
two maps suggest an unequal distribution of income when looking at the distribution of 
population in Paris suburbs. 


### 1.3.3  Ratio 

````{r ratio_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# ratio
com$ratio <- com$INC / com$TH

# ratio map
mf_map(x = com, var = "ratio", type = "choro",
       breaks = c(min(com$ratio, na.rm = TRUE), 
                  20000, 30000, 40000, 50000, 60000, 
                  max(com$ratio, na.rm = TRUE)),
       pal = c("#FCDACA", "#F6A599", "#F07168", "#E92B28", "#C70003", "#7C000C"),
       border = "white", lwd = 0.2, leg_pos = "topleft", leg_val_rnd = 0,
       leg_title = paste0("Average amount of income tax",
                          "\nreference per households\n(in euros)")) 

# EPT borders
mf_map(ept, col = NA, border = "black", add = TRUE)

# Label min and max
mf_label(x = com[which.min(com$ratio),], var = "LIBCOM", cex = 0.6, font = 2,
         halo = TRUE)
mf_label(x = com[which.max(com$ratio),], var = "LIBCOM", cex = 0.6, font = 2,
         halo = TRUE)

# layout 
mf_layout(title = "Ratio - Income per tax households, 2013",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
````

The _MGP_ area is characterized by high income inequalities. For the 150 municipalities 
of this area, the values extend from 14 730 _(La Courneuve)_ to 96 310 euros 
(*Paris_7th Arrondissement*). 53 municipalities of the _MGP_ area (35 % of the municipalities) 
are below the French average, i.e. 25 660 euros. The lagging households are mainly 
concentrated into the north part of _MGP_ area.  Highest values are concentrated 
in the Western part of Paris and its suburbs.  



## 1.4 Introducing the MTA Functions

### 1.4.1 General, territorial and spatial deviations
`MTA` package introduces three contexts to monitor territorial inequalities: the 
general, the territorial and the spatial deviations. 

The **global deviation** is dedicated to the analysis of inequalities using a 
value of reference. In this example the global deviation refers to the inequalities 
existing between each municipality in regard to the whole _MGP_ average.

The **territorial deviation** consists in measuring the inequalities existing for 
each basic territorial unit as regards to an intermediate territorial level of 
reference. In this case-study, it will be the _EPT_ of belonging of each municipality. It implies to include beforehand in the input dataset a territorial level of belong for each territorial unit. 

The **spatial deviation** is a measure of inequalities taking into account the neighborhood as a reference. It allows to measure the deviation existing between basic territorial units using three possible parameters : the territorial contiguity (order 1, 2 or n), the spatial neighborhood (territorial units located at less than X kilometers) or functional distances (territorial units located at less than Y minutes by road for instance). In `MTA`, territorial contiguity and spatial neighborhood measures are directly calculated. Functional distances must be uploaded separately through a dataframe structured with the following fields : id1, id2, distance measure. 
For this case-study, the contiguity criteria (order 1) will be used for the calculation of the spatial deviation. It gives a good proxy of proximity relationships according to the size and the homogeneity of the municipalities of the _MGP_ area. Other measures could be also adapted, such as time-distances by road (municipalities located at less than 15 minutes by car) or by public transportation.  


### 1.4.2 Relative and absolute deviations
In MTA, two methods are implemented to measure statistical differences to a given context of reference: **a relative deviation** and an **absolute deviation**. In `MTA`, each indicator is considered as a ratio defined by a numerator (GDP for instance) divided by a denominator (population for instance). 

The **relative deviation** states the position of each region as regards to a reference context. It is based on the following calculation:
Relative deviation (Region i) = 100 * ((Numerator(Region i)/Denominator(Region i)) /(reference ratio))
Territorial units characterized by a context of reference below index 100 are under the average of a given reference context, and reciprocally. 

The **absolute deviation** specifies which process of redistribution should be realized in absolute terms to achieve a perfect equilibrium for the ratio of reference in the global, the territorial or the spatial context. It is calculated as below:
Absolute deviation (Region i) = Numerator (Region i) - (reference ratio * denominator (Region i))
In concrete terms, it examines how much amount of the numerator should be moved in order to reach equidistribution, for each territorial unit, taking into account as a reference the selected deviation context value.
More generally, absolute deviation must be considered as a statistical tool to discuss on the amplitude of existing territorial inequalities. It is obvious that reaching a perfect equilibrium between territorial units is highly sensitive and is scarcely a policy objective itself. It is nevertheless interesting to consider the amount of money or population affected by the inequality to have in hand a concrete material for leading discussions on this delicate spatial planning issue. 


### 1.4.3 Cartography of relative and absolute deviations
In this vignette, color palette proposed for displaying on map the relative deviations are the ones suggested by the HyperAtlas tool (blue palette = under the average; red palette = above the average). But other diverging palettes (green/red, etc.) could be also be used for displaying MTA results on maps.  

Absolute deviations are highlighted using proportional circles on maps. It examines which amount of the numerator should be moved to the poorest municipalities to reach equidistribution. Circles displayed in a red means that the territorial unit have to contribute a given amount of numerator to achieve the equilibrium in a given context; and reciprocally blue circles means that the territorial unit have to receive frmo the others. 

It is quite easy to reproduce or adapt these maps using the functions proposed by the `mapsf` package.


### 1.4.4 Synthesis 
The analysis of the three deviations provides generally a lot of information. The synthesis functions `bidev` and `mst` helps to summarize values of the global, territorial and spatial deviations and allow to answer to some basic and interesting questions:

* Which territorial units are positioned in a favorable/lagging situation for 2 or 3 deviations?
* Which territorial units are positioned in a contradictory situations ? This is especially interesting from a social and policy points of view (NIMBY phenomena, characterizing the opposition by residents to a proposal for a new development because it is close to them). 

Moreover, some additional functions are implemented to ease the mapping (`map_bidev` and `map_mst`) or proposing appropriate plots (`plot_bidev` and `plot_mst`) for visualizing these results. 

<br>

# 2. Relative deviations

This section aims at describing how computing MTA relative deviations and creating appropriate maps and plots. 

## 2.1 Global deviation

First of all, each territorial unit is compared to the overall study area average (*Métropole du Grand Paris*). 

The code below takes in entry the numerator (INC) and the denominator (TH) of the `com` object and returns the global deviation indicators using `gdev()` function. The relative deviation (`type = "rel"`) is computed. This indicator is afterwards mapped. 

````{r gdevrel_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# general relative deviation
com$gdevrel <- gdev(x = com, 
                    var1 = "INC", 
                    var2 = "TH", 
                    type = "rel")

# Colors for deviations
devpal <-  c("#4575B4", "#91BFDB", "#E0F3F8", "#FEE090", "#FC8D59", "#D73027")

# Global deviation mapping
mf_map(x = com, var = "gdevrel", type = "choro",
       breaks = c(min(com$gdevrel, na.rm = TRUE), 67, 91, 100, 125, 150,
                  max(com$gdevrel, na.rm = TRUE)),
       pal = devpal, border = "white", lwd = 0.2, 
       leg_pos = "topleft", leg_val_rnd = 0,
       leg_title = paste0("Deviation to the global context",
                          "\n(100 = Metropole du Grand Paris average)")) 

# Plot EPT layer
mf_map(ept, col = NA, border = "#1A1A19", lwd = 1, add = TRUE)

# layout 
mf_layout(title = "Global deviation - Tax income per households",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
````

The resulting map highlights strong inequalities: First, it is intersting to know that all the municipalities of the _EPT_ of _Plaine Commune_, _Territoire des Aéroports_ and _Est Ensemble_ are below the average of the _MGP_ area (33 501 euros per households). For the territories of _Val de Bièvres_ and _Grand-Paris Est_, only three municipalities are above the average of the study area. Reversely, all the municipalities of the _Grand-Paris-Sud-Ouest EPT_ are largely above the average of the _MGP_ area. For the other _EPT_ the situation is mixed, depending of the municipalities. 

The `ineq` package proposes some functions to depict global inequalities existing in a study area, such as the Lorenz-curve. 
The `Lc` function takes in entry the numerator and the denominator and returns a Lorenz Curve plot; inequality indexes take in entry the ratio (numerator / denominator) and returns econometric indexes of inequality. 

The code below add some additional graphical parameters to ease the interpretation of this plot. 

````{r lorenz_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.align = 'center', fig.show='hold'}
library(ineq)
#  Concentration of X as regards to concentration of Y
Lc.p <- Lc(com$INC, com$TH)
Lp <- data.frame(cumX = 100 * Lc.p$L, cumY = 100 * Lc.p$p)

# Plot concentrations
opar <- par(mar = c(4,4,1.2,4), xaxs = "i", yaxs = "i", pty = "s")
plot(Lp$cumY, Lp$cumX,  type = "l", col = "red", lwd = 2,
     panel.first = grid(10,10), xlab = "Households (cumulative percentage)",
     ylab = "Income (cumulative percentage)", 
     cex.axis = 0.8, cex.lab = 0.9, ylim = c(0,100), xlim = c(0,100))  

lines(c(0,100), c(0,100), lwd = 2)

# Ease plot reading
xy1 <- Lp[which.min(abs(50 - Lp$cumX)),]
xy2 <- Lp[which.min(abs(50 - Lp$cumY)),]

xy <- rbind(xy1, xy2)

points(y = xy[,"cumX"], 
       x = xy[,"cumY"], 
       pch = 21,  
       cex = 1.5,
       bg = "red")

text(y = xy[,"cumX"], 
     x = xy[,"cumY"],
     label = paste(round(xy[,"cumX"],0), round(xy[,"cumY"],0), sep = " , "),  
     pos = 2,
     cex = 0.6)
par(opar)
````

The curve depicts on its horizontal axis a defined population – e.g., all households – broken down into deciles and ordered from left to right on the horizontal axis (from the lower tax income per household to the higher).
On the vertical axis of the Lorenz curve is shown the cumulative percentage of tax income. 

The reading of the plot reveals these following configurations :

* 50 % of the households earns less than 16 % of the total income.
* 50 % of the total income is held by less than 20 % of households. 


## 2.2 Territorial deviation

The territorial deviation is now calculated to analyze the position of each municipality as regards to a territorial level of reference: the _(EPT)_ of belonging in this case-study.  

The `tdev` function takes in entry the numerator (INC) and the denominator (TH) of the `com` object. The territorial level of belonging in specified with the `key` argument. The territorial relative deviation is returned by the function.     

````{r mdevrel_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# Territorial relative deviation calculation
com$tdevrel <- tdev(x = com, 
                    var1 = "INC", 
                    var2 = "TH", 
                    type = "rel",
                    key = "LIBEPT")

# Cartography
# Territorial deviation mapping
mf_map(x = com, var = "tdevrel", type = "choro",
       breaks = c(min(com$tdevrel, na.rm = TRUE), 67, 91, 100, 125, 150,
                  max(com$tdevrel, na.rm = TRUE)),
       pal = devpal, border = "white", lwd = 0.2, 
       leg_pos = "topleft", leg_val_rnd = 0,
       leg_title = paste0("Deviation to the territorial context",
                          "\n(100 = EPT average)")) 

# Plot EPT layer
mf_map(ept, col = NA, border = "#1A1A19", lwd = 1, add = TRUE)

# Labels to ease comment location
mf_label(x = com[com$LIBCOM %in% c("Le Raincy", "Rungis", "Sceaux",
                                   "Marnes-la-Coquette") ,],
         var = "LIBCOM", cex = 0.6, font = 2, halo = TRUE)

# layout 
mf_layout(title = "Territorial deviation - Tax income per households, 2013",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
````

This map highlights existing inequalities in each _EPT_: The strongest differences in relative terms are located in Paris (opposition between the eastern part and the western part of this _EPT_) and in the _Plaine centrale_ - _Haut Val de Marne EPT_ (opposition between the poorest municipalities located near Paris and the ones located in the periphery). Globally, the richest and the poorest _EPT_ ( _Grand Paris Sud Ouest_ / _Plaine Commune_ and _Territoires des Aéroports du Nord Ouest_) appear relatively homogeneous statistically. In other _EPT_, one municipality appears largely above the average of their _EPT_ of belonging. It is namely the case in _Grand-Paris-Sud-Ouest (Marnes-la-Coquette)_, _Sud Hauts-de-Seine (Sceaux)_, _Grand Orly Seine Biève (Rungis)_ or _Grand Paris Grand Est (Le Raincy_)

Another way to explore characteristics of the territorial deviation consists in analyzing the statistical dispersion (general deviation) by intermediate level _(EPT_ in this case). The best suited graphical representation for this is certainly a boxplot. 

The code below takes in entry the general deviation calculated above and the intermediate levels included in the input dataset. It returns a boxplot displaying the statistical parameters (median, mean, 1st and 3rd quartiles, range, minimum and maximum, extraordinary values) allowing to observe the statistical dispersion existing for each intermediate zoning. To ease the interpretation and the synthesis of the plot, boxplots are ordered by the average of each territorial level. Moreover, the width of the bars are proportional to the number of territorial units included in each intermediate zoning. 

````{r mdev_boxplot, fig.width=7, fig.height=6, warning = FALSE, cache = FALSE, fig.align = 'center', fig.show='hold'}
opar <- par(mar = c(4, 4, 1.2, 2))

# Drop geometries
df <- st_set_geometry(com, NULL)

# Reorder EPT according to gdev value
df$EPT <- with(df, reorder(EPT, gdevrel, mean, na.rm = TRUE))

# Colors management
col <- aggregate(x = df[,"gdevrel"], by = list(LIBEPT = df$LIBEPT),
                 FUN = mean)
col <- merge(col, colEpt, by = "LIBEPT")
col <- col[order(col$x),]
cols <- as.vector(col$cols)

# Drop inexisting levels
df <- droplevels(df)

# Boxplot
bp <- boxplot(df$gdevrel ~ df$EPT, col = cols, ylab = "Global deviation",
              xlab = "Territorial deviation", cex.lab = 0.9,
              varwidth = TRUE, range = 1.5,  outline = TRUE,  las = 1) 

# Horizontal Ablines
abline (h = seq(40, 300, 10), col = "#00000060", lwd = 0.5, lty = 3)

# Plot mean values
xi<- tapply(df$gdevrel, df$EPT, mean, na.rm = TRUE)
points(xi, col = "#7C0000", pch = 19)

# Legend for the boxplot
legend("topleft", legend = rev(as.vector(col$LIBEPT)), pch = 15,
       col = rev(as.vector(col$cols)), cex = 0.8, pt.cex = 1.5)
par(opar)
````

This plot highlights the statistical dispersion existing in each EPT. It confirms globally that wealthier the _EPT_ is, larger the statistical differences between the poorest and the wealthiest territorial units are. In that way, _Plaine Commune_ and _Territoire des Aéroports_ (T6 and T7) appear quite homogeneous in a lagging situation. On the opposite, important differences a revealed in  Paris (minimum = 71 and maximum = 290) or in _La Defense_ and _Grand Paris Sud Ouest_. 

The boxplot highlights also outliers (dots out of the box) in each territorial context. Most _EPT_ are concerned: wealthiest _EPT_ (Paris, _La Défense_, _Grand Paris Sud Ouest_ and _ACEP_, generally characterized by high outliers) and _EPT_ in less favorable situation (_Est Ensemble_, _Grand Paris Est_ and _Sud-Hauts-de-Seine_, which include municipalities with include both low and high outliers. 


## 2.3 Spatial deviation

The spatial deviation is calculated to position territorial units in a neighborhood context. Several criteria may be considered (geographical distance, time-distance). In this example, each municipality will be compared to the average of contiguous municipalities (contiguity order 1).  

The `sdev` function takes in entry the numerator (INC) and the denominator (TH) of the `com` object. The contiguity order 1 is set in the `order` argument. 

````{r localdevrel_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# Spatial relative deviation calculation
com$sdevrel <- sdev(x = com, 
                    xid = "DEPCOM", 
                    var1 = "INC", 
                    var2 = "TH",
                    order = 1,
                    type = "rel")

# Cartography
# Territorial deviation (relative and absolute) cartography
mf_map(x = com, var = "sdevrel", type = "choro",
       breaks = c(min(com$sdevrel, na.rm = TRUE), 67, 91, 100, 125, 150,
                  max(com$sdevrel, na.rm = TRUE)),
       pal = devpal, border = "white", lwd = 0.2, 
       leg_pos = "topleft", leg_val_rnd = 0,
       leg_title = paste0("Deviation to the spatial context",
                          "\n(100 = average of the contiguous",
                          " territorial units - order 1)")) 

# Plot EPT
mf_map(ept, col = NA, border = "#1A1A19",lwd = 1, add = T)

# Labels to ease comment location
mf_label(x = com[com$LIBCOM %in% c("Le Raincy", "Vaucresson", "Sceaux", "Bagneux",
                                   "Marnes-la-Coquette", "Saint-Maur-des-Fosses",
                                   "Puteaux", "Saint-Ouen", "Clichy-sous-Bois", 
                                   "Clichy"),],
         var = "LIBCOM", cex = 0.6, font = 2, halo = TRUE, overlap = FALSE)

# layout 
mf_layout(title = "Spatial deviation - Tax income per households, 2013",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
````

This map highlights local discontinuities existing in the study area. Important local statistical gaps appear in several areas: _Saint-Mandé_ and _Neuilly-sur-Seine_ are characterized by the highest score as regards to their respective neighbors (indexes 173 and 156).Some local _"bastions"_ are revealed in Paris suburbs, such as _Le Raincy_ (index 153), _Sceaux_ (150), _Vaucresson_ (143), _Marnes-la-Coquette_ (142) or _Saint-Maur-des-Fossés_ (140).

On the reverse situation, lower index are observed at _Clichy-sous-Bois_. _Clichy_, _Puteaux_, _Saint-Ouen_ and _Bagneux_. Their average income per household is around 30-40 % below their respective neighborhood.

<br> 
Spatialized indicators are often subject to spatial dependencies (or interactions), which are even stronger than spatial locations are closer. Autocorrelation measures (Moran, Geary, Lisa indexes) allows estimating the spatial dependence between the values of a same indicator at several locations of a given study area. 

The figure displayed below consists in evaluating this spatial autocorrelation by a plot crossing the spatial deviation (Y axis) and global deviation (X axis) values. 

This plot provides interesting inputs for answering to basic questions, such as:
* Are territorial units in favorable situation in a global context also in the same situation at local level?
* Is it possible to detect specific situations as regard to the global trend, meaning the existence of "local bastions" in favorable or lagging situations? 

````{r spat_autocorr_plot, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.align = 'center',  fig.show='hold'}
opar <- par(cex.lab = 1, cex.axis = 0.75, mar = c(4, 4, 1.2, 2))

# Drop geometries
df <- st_set_geometry(com, NULL)

# Spatial autocorrelation
lm <- summary.lm(lm(sdevrel ~ gdevrel, df))

# Equation 
eq <- paste("Spatial Deviation =", 
            round(lm$coefficients["gdevrel","Estimate"], digits = 3),
            "* (Global Deviation) +", 
            round(lm$coefficients["(Intercept)","Estimate"], digits = 3))
rsq <-paste("R-Squared =", 
            round(summary(lm(sdevrel ~ gdevrel, com ))$r.squared, digits = 2))

# Color management
df <- merge(df, colEpt, by = "LIBEPT")

# Plot spatial autocorrelation
plot(df$gdevrel, df$sdevrel,
     ylab = "Local deviation",
     ylim = c(50,260),
     xlab = "Global deviation",
     xlim = c(50,260),
     pch = 20,
     col = as.vector(df$col),
     asp = 1)
abline((lm(df$sdevrel ~ df$gdevrel)), col = "red", lwd =1)

# Specify linear regression formula and R-Squared of the spatial autocorrelation
text(110,60, pos = 4, cex = 0.7, labels = eq)
text(110,55, pos = 4, cex = 0.7, labels = rsq)

abline (h = seq(40,290,10), col = "gray70", lwd = 0.25, lty = 3)
abline (h = seq(50,250,50), col = "gray0", lwd = 1, lty = 1)
abline (v = seq(40,290,10), col = "gray70", lwd = 0.25, lty = 3)
abline (v = seq(50,250,50), col = "gray0", lwd = 1, lty = 1)

# Legend for territorial level
legend("topleft", legend = rev(as.vector(colEpt$LIBEPT)), pch = 15,
       col = rev(as.vector(colEpt$cols)), cex = 0.6, pt.cex = 1.5)
par(opar)
````

The output of the linear model of spatial autocorrelation reveals that the hypothesis of independence is rejected at a probability below than 0.0001. It means that, "everything is related to everything else, but near things are related than distant things" (Tobler, 1970). However, the R-squared of the relation (0.42) suggests that this statistical relation includes outliers very far from the linear regression.

This chart can be modified for analyzing outliers specifically. The code below computes the statistical residuals of the spatial autocorrelation calculated above. Then the residuals are standardized. Finally, a plot is built to display and label significant residuals. 

````{r spat_autocor_res_plt, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.align = 'center',  fig.show='hold'}
opar <- par(cex.lab = 1, cex.axis = 0.75, mar = c(4, 4, 2, 2))

# Standardized residual calculation
lm <- lm(sdevrel ~ gdevrel, df)
df$res <- rstandard(lm)

#risk alpha (0.1 usually)
alpha <- 0.055

# Calculation of the threshold using T-Student at (n-p-1) degrees of freedom
thr <- qt(1 - alpha / 2, nrow(com) - 1)

# Plot residuals
plot(df$sdevrel, df$res,
     xlab = "Local deviation", cex.lab = 0.8,
     ylim = c(-3.5, 3.5),
     xlim = c(40, 200),
     ylab = "Standardized residuals of spatial autocorrelation", 
     cex.lab = 0.8,
     cex.axis = 0.8,
     pch = 20,
     col = as.vector(df$cols))

# Adding thresholds
abline(h = - thr, col = "red")
abline(h = + thr, col = "red")
abline(h = 0, col = "red")

# Detecting exceptional values and labeling them on the plot
ab <- df[df$res < -thr | df$res > thr,]

# Plot residual labels
text(x = ab[,"sdevrel"], y = ab[,"res"], ab[,"LIBCOM"], cex = 0.5, pos = 4)

abline (v = seq(50, 200, 10), col = "gray70", lwd = 0.25, lty = 3)
abline (v = seq(50, 200, 50), col = "gray0", lwd = 1, lty = 1)

# Plot the legend (territorial zoning)
legend("topleft", legend = rev(as.vector(colEpt$LIBEPT)), pch = 15,
       col = rev(as.vector(colEpt$cols)), cex = 0.6, pt.cex = 1.5)
par(opar)
````

Further analysis may be considered in the domain, using for instance Moran indexes et LISA indexes or other distance criteria. This, it is important to remind that the choice of the appropriate threshold or criteria must make sense from a thematic point of view (does a neighborhood of 5 km mean something? Is a distance of 20 minutes by car more appropriated?)

<br>

# 3. Redistributions

Redistributions highlight another way to analyze territorial disparities. The analysis is no more focused on inequalities themselves, but rather than to the statistical process (amount of numerator or denominator) required to reach the equilibrium in each territorial context. 

The meaning of this mathematical calculation is strong in policy term, since it point out the municipalities which may have to contribute / receive the highest to ensure a perfect equidistribution in several territorial contexts. It gives also the amount of the reallocation process.  

The calculation is made for the general and the territorial contexts, using the argument `type = abs`. Two maps are created, displaying the amount of income that should be transfered from the wealthiest to the poorest municipalities. 

````{r redistributions, fig.width=7, fig.height=5.5,  fig.show='hold'}
# general absolute deviation 
com$gdevabs <- gdev(x = com, 
                    var1 = "INC", 
                    var2 = "TH", 
                    type = "abs")

# Territorial absolute deviation calculation
com$tdevabs <- tdev(x = com, 
                    var1 = "INC", 
                    var2 = "TH", 
                    type = "abs",
                    key = "LIBEPT")

# Transform the values in million Euros
com$gdevabsmil <- com$gdevabs / 1000000
com$tdevabsmil <- com$tdevabs / 1000000

# Deviation orientation
com$gdevsign <- ifelse(com$gdevabsmil> 0, "Income surplus", "Income deficit")
com$tdevsign <- ifelse(com$tdevabsmil > 0, "Income surplus", "Income deficit")

# Deviation maps 
opar <- par(mfrow = c(1,2))

# General deviation
# Plot territories
mf_map(com, col = "peachpuff", border = "white", lwd = 0.25)
mf_map(ept, col = NA, border = "#1A1A19", lwd = 1, add = TRUE)

mf_map(x = com, var = c("gdevabsmil", "gdevsign"), type = "prop_typo",
       leg_title = c("Absolute Deviation\n(Income redistribution, euros)",
                    "Redistribution direction"),
       leg_pos = c("topleft", NA), leg_val_rnd = -2,
       val_order = c("Income deficit", "Income surplus"),
       pal =  c("#ff0000","#0000ff"), border = "grey70", 
       add = TRUE)
       
# Labels to ease comment location
mf_label(x = com[com$LIBCOM %in% c("Paris 7e Arrondissement",
                                   "Neuilly-sur-Seine", "Aubervilliers") ,],
         var = "LIBCOM", cex = 0.6, font = 2,
         halo = TRUE, overlap = FALSE)

# Layout map 1
mf_layout(title = "General deviation (Metrople du Grand Paris)",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)

# Territorial deviation
mf_map(com, col = "peachpuff", border = "white", lwd = 0.25)
mf_map(ept, col = NA, border = "#1A1A19", lwd = 1, add = TRUE)

mf_map(x = com, var = c("tdevabsmil", "tdevsign"), type = "prop_typo",
       leg_title = c("Absolute Deviation\n(Income redistribution, euros)",
                    "Redistribution direction"),
       leg_pos = c(NA, "topleft"), leg_val_rnd = -2,
       val_order = c("Income deficit", "Income surplus"),
       val_max = max(abs(com$gdevabsmil)),
       pal =  c("#ff0000","#0000ff"), border = "grey70", 
       add = TRUE)
       
# Labels to ease comment location
mf_label(x = com[com$LIBCOM %in% c("Marnes-la-Coquette",
                                   "Nanterre", "Clichy-sous-Bois") ,],
         var = "LIBCOM", cex = 0.6, font = 2,
         halo = TRUE, overlap = FALSE)

# Layout map 2
mf_layout(title = "Territorial deviation (EPT of belonging)",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
par(opar)
````

In the context of the _MGP_, the *7th Arrondissement of Paris* is the municipality which should contribute the most, all things being equal to its income tax level (1,987 billion Euros of income transfer, 65 % of the total amount of tax income in this municipality). _Neuilly-sur-Seine_ (third position in absolute terms) should transfer 1,9 billion Euros to the lagging municipalities of the _MGP_ area. It corresponds to 62,87% of the total amount of tax income declared in this municipality.

In the other side of the redistribution, La Courneuve should receive 402 million Euros from the wealthiest municipalities (127 % of the current tax income declared). _Aubervilliers_ should receive 793 million Euros, which represents 124 % of the current tax income in this municipality.

The chunks below compute the top 10 theoretical contributors and receivers municipalities (total amount of income and share of current available income) to reach a perfect equilibrium. 

````{r gdev_listing, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# general deviation - Top 10 of the potential contributors in regard 
# Drop geometries
df <- st_set_geometry(com, NULL)
row.names(df) <- df$LIBCOM

# to their total amount of income
df$gdevabsPerc <- df$gdevabs / df$INC * 100
df <- df[order(df$gdevabsPerc, decreasing = TRUE), ]
df[1:10, c("gdevabsmil","gdevabsPerc")]

# general deviation - Top 10 of the potential receivers in regard to 
# their total amount of income
df <- df[order(df$gdevabsPerc, decreasing = FALSE), ]
df[1:10, c("gdevabsmil", "gdevabsPerc")]
````

Looking at the EPT context, the *7th Arrondissement of Paris* is still the municipality which should contribute the most to the poorest municipalities of Paris as regards to the amount of income available in this municipality (1,779 billion Euros of income transfer, 58 % of the amount of income in this municipality). _Marnes-La-Coquette_ (second position) should transfer 38 million Euros to the poorest municipalities of its _EPT_ of belonging _(La Défense)_. Despite the low income mass, it corresponds to 52 % of the total amount of income available in this municipality.

From the other side of the redistribution, _Nanterre_, _Clichy-sous-Bois_ and the *19th Arrondissement of Paris* should receive respectively 1088, 143 and 1926 million euros from the wealthiest municipalities of their _EPT_ of belonging (respectively 88 %, 68 % and 65 % of the total amount of earned income of their households). The highest redistribution for this study area stands for the *20th Arrondissement of Paris* (1,926 billion Euros, 59 % of its total amount of income).

````{r mdev_listing, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE}
# general deviation - Top 10 of the potential contributors in regard 
# Drop geometries
df <- st_set_geometry(com, NULL)
row.names(df) <- df$LIBCOM

# Territorial deviation - Top 10 of the potential contributors
# as regards to their total amount of income
df$tdevabsPerc <- df$tdevabs / df$INC * 100
df <- df[order(df$tdevabsPerc, decreasing = TRUE), ]
df[1:10, c("tdevabsmil", "tdevabsPerc")]

# Territorial deviation - Top 10 of the potential receivers 
# as regards to their total amount of income
df <- df[order(df$tdevabsPerc, decreasing = FALSE), ]
df[1:10, c("tdevabsmil", "tdevabsPerc")]
````

<br>


# 4. Synthesis

This section shows the synthesis functions of the `MTA` package: `bidev`, `plot_bidev` and `map_bidev` are useful to synthetize, plot and map the situation for 2 selected relative deviations. The `mst`, `plot_mst` and `map_mst` are suitable for 3 deviations.   


## 4.1 Two-relative deviations synthesis

The `bidev` function allows firstly to position each territorial unit as regards to 2 deviations (`dev1` and `dev2`) already calculated; and secondly to evaluate the statistical distance to the average (index 100). This function returns a vector which is the combination of two modalities :

* **A**: deviations above index 100 for the two selected contexts. 
* **B**: above 100 for `dev1` and below 100 for `dev2`.
* **C**: below 100 for `dev1` and `dev2`.
* **D**: below 100 for `dev1` and above 100 for `dev2`.

Are combined to modalities A, B, C and D the distance to the average. By default, it corresponds to 25 %, 50 % and 100 % under-above the average. Adapted to normalized index where 100 is average, it creates 4 classes : 

* **ZZ** : Values below 25 % of the average, corresponding to indexes 80 and 125).
* **1** : Values between 25 % and 50 % of the average, corresponding to indexes 67-80 and 125-150.
* **2** : Values between 50 % and 100 % of the average, corresponding to indexes 50-67 and 150-200.
* **3** : Values above 50 % of the average, corresponding to indexes below 50 and above 200. 

At this end, it produces a 13 classes-typology, which is the combination of the two modalities described above (classes A1, B3, C2, etc.). The user is free to modify the breaks proposed by default by the function (`breaks = c(25, 50, 100)`)

To understand the resulting typology, the `plot_mst` is especially adapted. The X-Y scale is defined in logarithm. Index 200 corresponds to territorial units defined to twice the average ; and index 50 twice below the average. 

````{r bidev_plot, fig.width=7, fig.height=8, fig.show='hold'}
opar <- par(mar = c(0, 0, 0, 0))

# Prerequisite  - Compute 2 deviations
com$gdev <- gdev(x = com, var1 = "INC", var2 = "TH")
com$tdev <- tdev(x = com, var1 = "INC", var2 = "TH", key = "EPT")

# EX1 standard breaks with four labels
plot_bidev(x = com, 
           dev1 = "gdev", 
           dev2 = "tdev",
           dev1.lab = "General deviation (MGP Area)",
           dev2.lab = "Territorial deviation (EPT of belonging)",
           lib.var = "LIBCOM",
           lib.val =  c("Marolles-en-Brie", "Suresnes", 
                        "Clichy-sous-Bois", "Les Lilas"))
par(opar)
````

This plot highlights specific areas of interest. In the red quarter, territorial units above the average for the general and territorial deviations, such as _Marolles-en-Brie_(class A3). Reversely, municipalities appearing in the blue quarter are below the average for the 2 deviations, such as _Clichy-sous-Bois_(class C3). 

It shows also territorial units in contradictory situations, such as _Les Lilas_(class D2, green quarter), which is in lagging situation for the general deviation and in favorable situation for the territorial context ((EPT of belonging). Reversely, _Suresnes_ is in favorable situation in the context of the MGP and in lagging situation in its EPT of belonging.   

The `map_bidev` function delivers a list of 3 objects useful to map the results: a `sf` object (`geom`) including the result of the `bidev` function (new column). This object is ordered following the 13 classes. It delivers also a color vector for mapping purpose. It includes only resulting categories: if the dataset does not include "D3" class, the resulting `cols` object will not include dark green color.

For mapping the results, it is recommended to share the plot in two columns : one side for the map, and one side displaying the bidev plot to understand correctly the colors displayed on the map.  

````{r bidev_map, fig.width=7, fig.height=5,  fig.show='hold'}
opar <- par(mfrow = c(1,2))

bidev <- map_bidev(x = com, dev1 = "gdev", dev2 = "tdev",
                   breaks = c(50, 100, 200))

# Unlist resulting function
com <- bidev$geom
cols <- bidev$cols

# Cartography
mf_map(x = com, var = "bidev", type = "typo", 
       pal = cols, lwd = 0.2, border = "white", leg_pos = NA,
       val_order = unique(com$bidev))

mf_map(ept, col = NA, border = "#1A1A19", lwd = 1, add = TRUE)

# Label territories in the C3 category
mf_label(com[com$bidev == "C3",],
         var = "LIBCOM", halo = TRUE)

mf_title( "2-Deviations synthesis:", tab=FALSE)
mf_layout(title = "",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)
mf_title("general and territorial contexts", tab = FALSE, inner = T)

#Associated plot
plot_bidev(x = com,  dev1 = "gdev",  dev2 = "tdev", 
           dev1.lab = "General deviation (MGP Area)",
           dev2.lab = "Territorial deviation (EPT of belonging)",
           lib.var = "LIBCOM", lib.val = "Clichy-sous-Bois", cex.lab = 0.8)
par(opar)
````

This map highlights areas in favorable situation for the MGP area and EPT contexts (red), in lagging situation for these two contexts (blue) or in contradictory situations (yellow and green). Darker is the color, stronger are the distance to the average. 


## 4.2 Three-relative deviations synthesis

The three relative deviations (general, territorial, spatial) can be summarized in a synthetic typology using the `mst` function. It allows, according to a pre-defined threshold (100 being the average), to higlight the territorial units which are under or above this threshold, and for which of the 3 deviations proposed by the MTA package.  

Assuming (for instance) that the treshold is equal to 100 (value corresponding to the average for each deviations) and the superior argument is false (consider the value below the average), the `mst()` function returns a 8-classes typology (mst column) which can be interpreted as follows:

* Class 7: A territorial unit described by a mst value equal to 7 is under (`superior = FALSE`) the average (`threshold = 100`) for all the deviations.
* Classes 1, 2 and 4 : The territorial unit is under the average for one deviation only (1: general deviation, 2: territorial deviation and 4: spatial deviation).
* Classes 3, 5 and 7: The territorial unit is under the average for 2 of the 3 possible deviations (3: general and territorial deviations, 5: general and spatial deviations and 6: territorial and spatial deviations). 
* Class 0: The territorial unit is under the average for none of the deviations, meaning that it is above the average for all the deviations. 

As demonstrated below, this typology is especially useful to highlight territories in lagging or favorable situations, but also territories in contradictory situations. 


### 4.2.1 Wealthiest municipalities

We focus firtstly the analysis on municipalities above 125 % for the three deviations (gdev = *Métropole du Grand Paris* ; tdev = EPT average and sdev = contiguous municipalities, threshold = 125, superior = TRUE). The `mst()` function returns their position according to the three deviations. The subset (mst == 7) returns all the municipalities in this situation. 

````{r synth7_1, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, fig.align = 'center'}
# Prerequisite  - Compute the 3 deviations
com$gdev <- gdev(x = com, var1 = "INC", var2 = "TH")
com$tdev <- tdev(x = com, var1 = "INC", var2 = "TH", key = "EPT")
com$sdev <- sdev(x = com, var1 = "INC", var2 = "TH", order = 1)

# Multiscalar typology - wealthiest territorial units
# Row names = municipality labels
row.names(com) <- com$LIBCOM

# Compute mst
com$mstW <- mst(x = com, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev", 
                threshold = 125, superior = TRUE)

subset(com, mstW == 7, select = c(ratio, gdev, tdev, sdev, mstW), drop = T)
````

The results are afterwards mapped. As `map_bidev` function do, the `map_mst` function takes in entry the 3 pre-calculated deviations and `mst` parameters and returns a list including a `sf` object with the result of the `mst` calculation, and two vectors proposing colors and legend for mapping results.   

````{r map_mst1, fig.width=7, fig.height=5}
# Compute mapmst
mst <- map_mst(x = com, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev", 
               threshold = 125, superior = TRUE)

# Unlist resulting function
com <- mst$geom
cols <- mst$cols
leg_val <- mst$leg_val

# Cartography
mf_map(x = com, var = "mst", type = "typo", 
       border = "white", lwd = 0.2,
       pal = cols, val_order = unique(com$mst), leg_pos = NA)
mf_map(ept, col = NA, border = "black", lwd = 1, add = TRUE)
mf_legend(type = "typo", pos = "topleft", val = leg_val, pal = cols, 
          title = paste0("Situation on General (G)\n",
                         "Terrorial (T) and \n",
                         "Spatial (S) contexts"))
mf_layout(title = "3-Deviations synthesis: Territorial units above index 125",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)

# Add labels for mst = 7
mf_label(x = com[com$mst == 7,], var = "LIBCOM",  halo = TRUE, overlap = FALSE, 
         cex = 0.7)
````

This map highlights the municipalities 25 % above the average for the three deviations (colored in red). They are mainly located in the West part of Paris and of the _MGP_. Some isolated municipalities appear also locally in a favorable situation in the East Part of the _MGP_.

More precisely, 11 municipalities are characterized by an index of average income per household above 125% on the three contexts. It corresponds to 4 *Arrondissements* in Paris _(6e_, _7e_, _8e_ and _16e_) and the municipalities of _Marnes-la-Coquette_, _Neuilly-sur-Seine_, _Sceaux_, _Vaucresson_, _Marolles-en-Brie_, _Saint-Mandé_ and _Santeny_. 

13 municipalities can be considered as "globally" advantaged (above 125% only in a global and/or in a global and in territorial contexts, colored in light orange and in pink). It corresponds to the central *Arrondissements* in Paris _(1st_, _2nd_, _4th_, _9th_) and other municipalities in favorable situation close to wealthier municipalities _(Boulogne-Billancourt_, _Saint-Cloud_, etc).

8 municipalities can be considered as "locally" advantaged (index above 125% only in a spatial context and/or in a spatial and a territorial context, colored in green and yellow). These local poles of wealth are mainly located in the periphery of the _MGP_ area _(Le Raincy_, _Rungis_, _Coubron_, etc.), closed to poorer municipalities. 


### 4.2.2 Municipalities in lagging situation

It is also possible to analyze the reverse situation: the municipalities in lagging situation for this indicator. Thus, the `map_mst()` function takes '80' for the threshold argument and 'FALSE' for the superior argument. This typology will specify the municipalities situated below 25 % of the average for the global and/or the territorial and/or the spatial contexts.  

````{r map_mst2, fig.width=7, fig.height=5}
# Compute mapmst
mst <- map_mst(x = com, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev", 
              threshold = 80, superior = FALSE)

# Unlist resulting function
com <- mst$geom
cols <- mst$cols
leg_val <- mst$leg_val

# Cartography
mf_map(x = com, var = "mst", type = "typo", 
       border = "white", lwd = 0.2,
       pal = cols, val_order = unique(com$mst), leg_pos = NA)
mf_map(ept, col = NA, border = "black", lwd = 1, add = TRUE)
mf_legend(type = "typo", pos = "topleft", val = leg_val, pal = cols, 
          title = paste0("Situation on General (G)\n",
                         "Terrorial (T) and\n",
                         "Spatial (S) contexts"))
mf_layout(title = "3-Deviations synthesis: Territorial units under index 80",
          credits = paste0("Sources : GEOFLA® 2015 v2.1, Apur, impots.gouv.fr",
                           "\nRonan Ysebaert, RIATE, 2021"), 
          arrow = FALSE)

# Add labels for mst = 7
mf_label(x = com[com$mst == 7,], var = "LIBCOM",  halo = TRUE, overlap = FALSE, 
         cex = 0.7)
````

This map shows the municipalities in difficulty regarding the income distribution in the _MGP_. Interesting is to note that all the territories of the _Plaine Commune_, _Territoire des Aéroports_ _EPT_ are below the threshold of 80 % for at least the global context. Most of the municipalities of _EPT_ of _Est Ensemble_, _Grand-Paris-Est_, _Val de Bièvres_ are in the same situation.

The municipalities in lagging situation for the three contexts are located in several areas of the _MGP_ (North-West, East). All are not located not in the immediate periphery of Paris (none of these municipalities are contiguous to Paris). It concerns precisely 6 municipalities: _Bagneux_, _Genevilliers_, _Nanterre_, _Clichy-sous-Bois_, _Bonneuil-sur-Marne_ and _Champigny-sur-Marne_. _Clichy-sous-Bois_ appears especially in a dramatic situation (indexes 47, 60 and 61 for the global, territorial and spatial deviations). 

Most of the municipalities of _Seine-Saint-Denis_ and _Val-de-Marne_ are characterized by a significant lagging situation globally but not locally (classes 1 and 3, in light orange and pink). 

````{r synthesis80_class7}
# Multiscalar typology - Lagging territorial units
# Row names = municipality labels
row.names(com) <- com$LIBCOM

# Compute mst
com$mstP <- mst(x = com, gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev", 
                threshold = 80, superior = FALSE)

# municipalities in lagging situation for the three contexts
subset(com, mstP == 7, select = c(ratio, gdev, tdev, sdev, mstP), drop = T)
````


### 4.2.3 Contradictory situations

As displayed on the map above, some territories appear also in a contradictory situation depending on the deviation used. 

Municipalities characterized by class 3 in the mst synthesis are lagging as regards to the global and the territorial contexts, but not the spatial one. It is especially the arrondissements located in the North-West part of Paris, contiguous to poorest municipalities. 

````{r synthesis80_class3}
# municipalities in lagging situation in the global and territorial contexts
subset(com, mstP == 3, 
       select = c(ratio, gdev, tdev, sdev, mstP), drop = T)
# municipalities in favorable situation in a spatial context or in a spatial and a territorial context 
subset(com, mstP == 4 | mstP == 6, 
       select = c(ratio, gdev, tdev, sdev, mstP), drop = T)
````

4 municipalities are characterized by low indexes in a local context, but not in a global one. It corresponds to the municipalities located in the immediate neighborhood of the wealthiest municipalities of the MGP (classes 4 and 6, in green and yellow): _Paris 15e Arrondissement_, _Suresnes_, _Puteaux_ and in a lesser extent _Malakoff_.




### 4.2.4 Comparing specific territorial units

The `plot_mst` allows to compare selected territorial units on a barplot. This function takes in entry 3 pre-calculated deviations and returns a barplot which shows the position on the three deviations for one (or more) territorial unit(s).

````{r plot_mst, fig.width = 7, fig.height=6, fig.show='hold'}
opar <- par(mar = c(4, 6, 4, 4))
# Synthesis barplot
plot_mst(x = com, 
         gdevrel = "gdev", tdevrel = "tdev", sdevrel = "sdev", 
         lib.var = "LIBCOM", 
         lib.val = c("Neuilly-sur-Seine", "Clichy-sous-Bois", 
                     "Suresnes", "Les Lilas"))
par(opar)
````

The barplot above displays the situation for 4 municipalities of the study area ( _Les Lilas_, _Neuilly-sur-Seine_, _Suresnes_, _Clichy-sous-Bois_) is a good way to propose a visualization for comparing specific territorial units. 


## 4.3 Absolute deviations synthesis

The `mas()` function takes in entry all requested parameters to compute the three deviations, as specified above. It returns a dataframe summarizing the values of absolute deviations, e.g. how much should be redistributed from the poorest to the wealthiest territorial units to ensure a perfect equilibrium of the ratio for the three contexts. Results are expressed both in absolute values (mass of numerator, amount of tax reference in Euros in this case) and as a share of the numerator (x % of the numerator that should be redistributed).  

````{r synthesisabs, fig.width=7, fig.height=5, warning = FALSE, cache = FALSE, eval = TRUE}
# Local redistribution (not yet calculated)
com$sdevabs <- sdev(x = com, xid = "DEPCOM", var1 = "INC", var2 = "TH",
                    order = 1, type = "abs")
com$sdevabsmil <- com$sdevabs / 1000000

# Compute the synthesis DataFrame (absolute deviations)
mas(x = com, 
    gdevabs = "gdevabsmil", 
    tdevabs = "tdevabsmil",
    sdevabs = "sdevabsmil",
    num = "INCM") 
````

For the _MGP_ area, it is 22 billion Euros that should be redistributed from the municipalities in favorable situation to the municipalities in lagging situation. It corresponds to 16.3 % of the total mass of income declared to the taxes. If a policy option consists in ensuring a equilibrium at an intermediate territorial level, such as the _Établissements Publics Territoriaux_, it is 15 billion Euros that should be redistributed (11.6 % of the income mass). If a solution chosen consists in limiting territorial discontinuities in a local context (avoid local poles of wealth or of poverty), it is 10 billion Euros that should be redistributed (7.9 % of the income mass). Ensuring an equilibrium of income in these three territorial contexts are obviously not credible policy options, but it gives some references to monitor the magnitude of territorial inequalities existing in a given study area.   

<br>

# 5. Conclusion and perspectives

The functionalities provided by the `MTA` package are especially useful to enhance territorial inequalities in several perspectives. This example, applied to the *Métropole du Grand Paris* provides some interesting insights to understand challenges raised by income inequalities in this governance structure at the scale of the Parisian metropolitan area. Inputs provided by this vignette are quite simple and aims at using all the functionalities of the `MTA` package, by suggesting some adapted graphical representations. We consider also that it provides useful inputs for territorial monitoring and observation - a necessary step before policy decision making (which solutions for reducing income inequalities in a given area?)

It is also important to remind that the aim of this vignette is not to provide a global picture of territorial inequalities existing in the _MGP_. To have this ambition, further analysis should be led for interrogating other dimensions of the inequalities, such as:

* Comparison dimension: Comparing existing territorial inequalities in other similar metropolitan area (London for instance).
* Time dimension: Evaluating the evolution of territorial inequalities towards the time (are inequalities growing or not?)
* Thematic dimension: Integrating other indicators useful to monitor well-being at local level to this analysis (share of household owners, mean area of the households, structure of households)
* Geographical dimension: Extending the analysis at lower territorial level, such as the IRIS one. It could be useful to discuss on internal inequalities existing in the municipalities of the *Métropole du Grand Paris*. 

That being said, MTA analysis provide a useful methodological base to explore several dimensions of territorial inequalities measures and applications within a policy context:

* Simulation of policy options: it is quite frequent that indicators constitute a basis for implementing policy measures to reduce territorial disparities. The case of the regional policy of the EU is a good example taking into account that most of the funding goes to regions below the statistical threshold of 75 % of the average of the European Union for the GDP per capita criteria at NUTS2 level. The functionalities of the `MTA` package allow in that perspective to simulate quickly the consequences of the use of several alternatives for guiding the funds allocation.

* Highlight contradictions: the fact that MTA functionalities are based on 3 possible measures of territorial inequalities (general, territorial, spatial deviations) leads the analyst or the policy maker to think about several theories regarding the governance of territorial inequalities. A territorial unit situated in a lagging situation at general level and in a favorable situation at local level must not be considered in the same way than a territorial unit characterized by a lagging situation both at general and local levels.