Commit 6278ce19 authored by Hadrien Commenges's avatar Hadrien Commenges

autocor

parent 2391f397
---
title: "Autocorrélation spatiale"
author: "H. Commenges"
date: "2020"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, results = TRUE, warning = FALSE, message = FALSE, error = FALSE)
```
## Objectifs
- Créer des matrices de contiguité
- Calculer et interpréter des mesures d'autocorrélation spatiale
- Cartographier l'autocorrélation spatiale
## Données
- carroyage de 200 mètres produit par l'Insee sur l'espace Paris et petite couronne (départements 75, 92, 93, 94). La variable d'intérêt sera `REVUC` : le revenu moyen par unité de consommation dans chacun des carreaux
## Instructions
**0. Charger les packages nécessaires**, *a minima* le package `sf` et le package `spdep`
```{r}
library(sf)
library(spdep)
library(cartography)
library(leaflet)
```
**1. Charger les données**
```{r}
revCarr <- readRDS(file = "DATA/carroyage.Rds")
```
**2. Produire la liste des voisins à partir de l'objet spatial**
```{r}
listNei <- poly2nb(pl = revCarr, row.names = revCarr$idk, snap = 1, queen = TRUE)
head(listNei)
```
**3. Explorer la liste des voisins (dénombrer les voisins)**
```{r}
numberNei <- sapply(listNei, length)
table(numberNei)
barplot(table(numberNei),
border = "grey", bg = "grey",
main = "Distribution du nombre de voisins",
xlab = "Nombre de voisins",
ylab = "Nombre de communes")
```
**4. Transformer la liste de voisins en liste de poids**
```{r}
listWgt <- nb2listw(neighbours = listNei, style = "W", zero.policy = TRUE)
head(listWgt$weights)
```
**5. Afficher le graphique de Moran et calculer la statistique de Moran**
```{r, eval=FALSE}
moran.plot(revCarr$REVUC, listw = listWgt, quiet = TRUE, zero.policy = FALSE)
```
```{r}
moran.test(revCarr$REVUC, listw = listWgt, zero.policy = TRUE)
```
**6. Récupérer les valeurs moyennes chez les voisins (*lagged values*) et reproduire le graphique de Moran**
```{r}
revCarr$LAGREV <- lag.listw(x = listWgt, var = revCarr$REVUC, zero.policy = TRUE)
avgValue <- mean(revCarr$REVUC)
plot(revCarr$REVUC, revCarr$LAGREV, pch = 20, xlab = "Revenu par UC de la commune", ylab = "Revenu moyen dans le voisinage")
abline(h = avgValue, col = "firebrick")
abline(v = avgValue, col = "firebrick")
```
**7. Cartographier les cadrans en statique : haut-haut, bas-bas, haut-bas, bas-haut**
*Discrétiser sur la base du revenu et du revenu moyen chez les voisins et afficher une carte qualitative.*
```{r}
revCarr$TYPEGO <- ifelse(revCarr$REVUC > avgValue, 1, 0)
revCarr$TYPNEI <- ifelse(revCarr$LAGREV > avgValue, 1, 0)
revCarr$TYPE <- paste(revCarr$TYPEGO, revCarr$TYPNEI, sep = "_")
typoLayer(x = revCarr,
var = "TYPE",
col = c("black", "springgreen", "black", "firebrick"),
border = NA)
```
**8. Cartographier les cadrans avec leaflet**
```{r}
colPal <- c("firebrick", "black", "black", "springgreen")
leaflet() %>%
addProviderTiles("Stamen.TonerBackground") %>%
addPolygons(data = st_transform(revCarr, crs = 4326),
stroke = FALSE,
fillColor = colPal[as.numeric(as.factor(revCarr$TYPE))],
fillOpacity = 0.6)
```
------------
This diff is collapsed.
"CHAR1","CHAR2"
"lexi","sloan"
"lexi","karev"
"owen","yang"
"owen","altman"
"sloan","torres"
"sloan","altman"
"torres","arizona"
"torres","karev"
"derek","grey"
"karev","izzie"
"o'malley","izzie"
"torres","o'malley"
"yang","colin"
"yang","preston"
"karev","kepner"
"sloan","addison"
"karev","addison"
"derek","addison"
"sloan","nancy"
"karev","olivia"
"o'malley","olivia"
"grey","o'malley"
"karev","mrs. seabury"
"chief","adele"
"chief","ellis grey"
"ellis grey","thatch grey"
"susan grey","thatch grey"
"bailey","tucker"
"izzie","hank"
"izzie","denny"
"grey","finn"
"grey","steve"
"bailey","ben"
"lexi","avery"
---
title: "Analyse de graphe (introduction)"
author: "H. Commenges"
date: "2020"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, results = TRUE, warning = FALSE, message = FALSE, error = FALSE)
```
## Objectifs
- Apprendre les principales mesures globales de description d'un graphe
- Apprendre les principales mesures locales caractérisant les noeuds et arcs d'un graphe
- Se familiariser avec la structure de données de type `igraph`
## Données
- reśeau de relations sexuelles de la série [Grey's Anatomy](http://fr.wikipedia.org/wiki/Grey%27s_Anatomy). Exercice tiré du blog [BabelGraph](http://www.babelgraph.org/wp).
Plusieurs questions peuvent être posées à partir de ce réseau pour introduire les principales mesures de description du réseau, en particulier en termes de centralité. Si une maladie sexuellement transmissible (MST) se diffusait dans le réseau :
- à qui (**noeud**) devrait-on faire en priorité un test pour détecter la MST ?
- à qui (**noeud**) devrait-on en priorité donner un paquet de préservatifs (masculins ou féminins) ?
- qui (**noeud**) serait le/la plus à même d'échapper à la maladie ?
- quelle relation (**arc**) faudrait-il rompre pour limiter l'épidémie ?
## Instructions
**0. Charger les packages nécessaires**, *a minima* le package `igraph`
```{r}
library(igraph)
library(tidyverse)
```
**1. Charger les données**
```{r}
greysTab <- read.csv("DATA/ga_edgelist.csv", stringsAsFactors = FALSE)
greysNet <- graph.data.frame(greysTab, directed = FALSE)
summary(greysNet)
```
**2. Représenter le graphe en faisant varier les paramètre graphiques et l'algorithme de visualisation**
```{r}
plot(greysNet,
vertex.size = 10,
vertex.color = "grey",
vertex.frame.color = NA,
vertex.label.color = "black",
vertex.label.cex = 0.7,
edge.color = "firebrick",
edge.width = 2,
layout = layout_in_circle(greysNet))
plot(greysNet,
vertex.size = 10,
vertex.color = "grey",
vertex.frame.color = NA,
vertex.label.color = "black",
vertex.label.cex = 0.7,
edge.color = "firebrick",
edge.width = 2,
layout = layout_with_fr(greysNet))
```
**3. Calculer les mesures globales de description du graphe**
- nombre de noeuds
- nombre d'arcs
- nombre de composantes connexes
- rapport entre le nombre d'arcs et le nombre de noeuds (indice de complexité, indice $\beta$)
- diamètre du graphe : longueur du plus long des plus courts chemins
- rapport entre le nombre d'arcs et le nombre maximum d'arcs (indice de connectivité, indice $\gamma$)
- le nombre maximum d'arcs est égal **pour un graphe non planaire et non dirigé** à $\frac{V(V-1)}{2}$
```{r}
vcount(greysNet)
ecount(greysNet)
ecount(greysNet) / vcount(greysNet)
components(greysNet)
diameter(greysNet)
ecount(greysNet) / (vcount(greysNet) * (vcount(greysNet) - 1) / 2)
```
**3. Calculer les mesures locales de description du graphe**
- degré : nombre de voisins, nombre d'arcs incidents
- degré pondéré : nombre d'arcs incidents pondérés par un attribut (quel attribut ?)
```{r}
nodesDegree <- degree(greysNet)
summary(nodesDegree)
table(nodesDegree)
```
De la longueur des plus courts et du trajet des plus courts chemins sont dérivés deux indicateurs de centralité :
- la centralité de proximité (*closeness*)
- la centralité d'intermédiarité (*betweenness*)
```{r}
closeNodes <- closeness(graph = greysNet)
betweenNodes <- betweenness(graph = greysNet)
print(closeNodes)
print(betweenNodes)
```
La centralité d'intermédiarité peut aussi être calculée sur les arcs du réseau :
```{r}
betweenEdges <- edge.betweenness(greysNet)
print(betweenEdges)
```
**4. Représenter ces mesures sur le graphe**
```{r}
# proximité
plot(greysNet,
vertex.size = 5000 * closeNodes,
vertex.color = "grey",
vertex.frame.color = NA,
vertex.label = V(greysNet)$name,
vertex.label.color = "black",
vertex.label.cex = 0.7,
edge.color = "firebrick",
edge.width = 2)
# intermédiarité
plot(greysNet,
vertex.size = 0.1 * betweenNodes,
vertex.color = "grey",
vertex.frame.color = NA,
vertex.label = V(greysNet)$name,
vertex.label.color = "black",
vertex.label.cex = 0.7,
edge.color = "firebrick",
edge.width = 2)
# intermédiarité des arcs
plot(greysNet,
vertex.size = 3,
vertex.color = "grey",
vertex.frame.color = NA,
vertex.label = V(greysNet)$name,
vertex.label.color = "black",
vertex.label.cex = 0.7,
edge.color = "firebrick",
edge.width = 0.1 * betweenEdges)
```
**5. Répondre aux questions posées en introduction : quelle mesure répond à quelle question ?**
- à qui (**noeud**) devrait-on faire en priorité un test pour détecter la MST ?
- à qui (**noeud**) devrait-on en priorité donner un paquet de préservatifs (masculins ou féminins) ?
- qui (**noeud**) serait le/la plus à même d'échapper à la maladie ?
- quelle relation (**arc**) faudrait-il rompre pour limiter l'épidémie ?
------------
This diff is collapsed.
---
title: "Ségrégation"
author: "H. Commenges"
date: "2020"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, results = TRUE, warning = FALSE, message = FALSE, error = FALSE)
```
## Objectifs
- Calculer les principaux indicateurs de ségrégation
- Prendre conscience des apports et des limites de ces indicateurs
- Cartographier les indicateurs de ségrégation
Dans la littérature scientifique il existe un très grand nombre de mesures de ségrégation, et dans l'écosystème R il existe un grand nombre de fonctions déjà implémentées qui permettent de ces mesures. Parmi les packages spécifiques : `OasisR`, `seg`, `segregation`, `ineq`. Pour faciliter l'apprentissage, on calcule ici un faible nombre de mesures et on les calcule sans mobiliser de fonctions tout-en-un.
## Données
- données du recensement américain de 2010 au niveau des *census tracts* (données compilées par [Arthur Gailes](https://belonging.berkeley.edu/arthur-gailes)).
- fond de carte des *census tracts* (TIGER/Line, en ligne sur le [site du Census Bureau](https://www.census.gov)).
## Instructions
**0. Charger les packages nécessaires**, *a minima* les packages `sf`, `ineq`.
```{r}
library(sf)
library(leaflet)
library(ineq)
library(cartography)
library(reshape2)
library(tidyverse)
```
**1. Charger et représenter les données**
```{r}
sfba <- readRDS(file = "DATA/bayarea.Rds")
plot(sfba$geometry)
leaflet() %>%
addTiles() %>%
addPolygons(data = st_transform(sfba, crs = 4326),
color = "grey",
weight = 1,
fillColor = "grey",
fillOpacity = 0.3)
```
**2. Sélectionner le comté de San Francisco et explorer les distributions statistiques et spatiales des populations**
```{r}
sfCity <- sfba %>% filter(COUNTYFP10 == "075")
plot(sfCity$geometry)
```
```{r}
# calcul des pourcentages
sfCity <- sfCity %>%
mutate(PCTHIS = hispanic / total_pop * 100,
PCTWHI = white / total_pop * 100,
PCTBLA = black / total_pop * 100,
PCTASI = asian / total_pop * 100)
# distributions statistiques
hist(sfCity$PCTHIS, main = "Distribution du pourcentage d'hispanos", breaks = 15)
hist(sfCity$PCTWHI, main = "Distribution du pourcentage de blancs", breaks = 15)
hist(sfCity$PCTBLA, main = "Distribution du pourcentage de noirs", breaks = 15)
hist(sfCity$PCTASI, main = "Distribution du pourcentage d'asiatiques", breaks = 15)
# planche d'histogrammes avec ggplot2
pctLong <- melt(sfCity,
id.vars = "GEOID10",
measure.vars = c("hispanic", "white", "black", "asian"),
variable.name = "group",
value.name = "pct")
ggplot(pctLong) +
geom_histogram(aes(x = pct), color = "grey80", fill = "grey40") +
facet_wrap(~group) +
theme_bw()
# distributions spatiales
choroLayer(x = sfCity, var = "PCTHIS",
method = "quantile", nclass = 6,
col = carto.pal(pal1 = "wine.pal", n1 = 6),
border = "grey", lwd = 0.5,
legend.title.txt = "Pourcentage d'hispanos",
legend.pos = "topleft")
choroLayer(x = sfCity, var = "PCTASI",
method = "quantile", nclass = 6,
col = carto.pal(pal1 = "turquoise.pal", n1 = 6),
border = "grey", lwd = 0.5,
legend.title.txt = "Pourcentage d'asiatiques",
legend.pos = "topleft")
```
**3. Calculer et représenter la distribution spatiale différentielle des hispanos par rapport à la distribution globale**
Calculer le pourcentage de la population hispano qui réside dans chaque zone (pourcentage en colonne : nombre d'hispanos de la zone i rapporté au nombre total d'hispanos), puis calculer le pourcentage de la population non hispano qui réside dans chaque zone (pourcentage en colonne), puis calculer pour chaque zone la différence entre les deux pourcentages. Enfin, représenter cette différence.
```{r}
# calcul des proportions et différence
sfCity$total_nonhis <- sfCity$total_pop - sfCity$hispanic
sfCity$PROPNONHIS <- 100 * sfCity$total_nonhis / sum(sfCity$total_nonhis)
sfCity$PROPHIS <- 100 * sfCity$hispanic / sum(sfCity$hispanic)
sfCity$PROPDIFF <- sfCity$PROPHIS - sfCity$PROPNONHIS
hist(sfCity$PROPDIFF)
# représentation cartographique de la différence
choroLayer(x = sfCity, var = "PROPDIFF",
breaks = c(-1, -0.5, 0, 0.5, 1, 2, 3),
col = carto.pal(pal1 = "blue.pal", n1 = 2, pal2 = "red.pal", n2 = 4, middle = FALSE),
border = "grey", lwd = 0.5,
legend.title.txt = "Différentiel hispano",
legend.values.rnd = 2,
legend.pos = "topleft")
```
**4. Calculer et représenter la courbe de ségrégation (*segregation curve*)**
```{r}
# calcul des pourcentages cumulés
segCurve <- sfCity %>%
st_drop_geometry() %>%
select(GEOID10, PROPHIS, PROPNONHIS) %>%
arrange(PROPHIS) %>%
mutate(CUMNONHIS = cumsum(PROPNONHIS),
CUMHIS = cumsum(PROPHIS))
# représentation avec plot()
plot(segCurve$CUMNONHIS, segCurve$CUMHIS, type = "l",
xlab = "Proportion cumulée de non hispanos",
ylab = "Proportion cumulée d'hispanos")
abline(a = 0, b = 1)
# représentation avec ggplot()
ggplot(segCurve) +
geom_line(aes(CUMNONHIS, CUMHIS), color = "firebrick") +
geom_abline(slope = 1, intercept = 0, color = "grey30") +
xlab("Proportion cumulée de non hispanos") +
ylab("Proportion cumulée d'hispanos") +
theme_bw()
```
**5. Calculer l'indice de Hoover (ou D de dissimilarité)**
Demi-somme de la valeur absolue des différences de proportion (valeur entre 0 et 1 ou entre 0 et 100 si les pourcentages ont été multipliés par 100).
```{r}
# calcul de Hoover (ou Duncan)
sum(abs(sfCity$PROPDIFF / 100)) / 2
```
**6. Calculer et représenter l'indice d'entropie relative pour quantifier le mélange des groupes**
```{r}
# implémentation de la fonction
compute_entropy <- function(x){
freqTab <- x / sum(x)
obsEntropy <- -sum(freqTab * log2(freqTab))
maxEntropy <- log2(length(x))
relEntropy <- obsEntropy / maxEntropy
return(relEntropy)
}
# transformation du tableau en matrice d'effectifs
matGroups <- sfCity %>%
select(4:8) %>%
st_drop_geometry() %>%
as.matrix()
# application de la fonction ligne-à-ligne (row-wise)
sfCity$RELH <- apply(matGroups, 1, compute_entropy)
# représentation cartographique
choroLayer(x = sfCity, var = "RELH",
method = "quantile", nclass = 6,
col = carto.pal(pal1 = "sand.pal", n1 = 6),
border = "grey",
lwd = 0.5,
legend.title.txt = "Entropie relative",
legend.values.rnd = 2,
legend.pos = "topleft")
```
**7. Représenter l'indice d'entropie relative dans Leaflet**
```{r}
# création d'une palette de couleurs
colPal <- carto.pal(pal1 = "sand.pal", n1 = 6)
brks <- quantile(sfCity$RELH, probs = seq(from = 0, to = 1, by = 0.2))
sfCity$RELHDISCRET <- cut(sfCity$RELH, breaks = brks, labels = 1:5, include.lowest = TRUE, right = FALSE)
leaflet() %>%
addProviderTiles("Stamen.Toner") %>%
addPolygons(data = st_transform(sfCity, crs = 4326),
color = "grey",
weight = 0.8,
fillColor = colPal[sfCity$RELHDISCRET],
fillOpacity = 0.7,
label = paste0("Entropie relative : ", round(sfCity$RELH, 2)))
```
**9. Calculer l'indice de Hoover standardisé**
```{r}
allGeoids <- rep(x = sfCity$GEOID10, times = sfCity$total_pop)
allHis <- rep("hispanic", sum(sfCity$hispanic))
allNonhis <- rep("total_nonhis", sum(sfCity$total_nonhis))
length(allGeoids) == length(allHis) + length(allNonhis)
tabShuffle <- data.frame(GEOIDS10 = allGeoids, GROUP = c(allHis, allNonhis))
# implémentation de la fonction pivot et calcul de dissimilarité
cast_and_compute_dissimilarity <- function(x, idvar, groupvar){
tabWide = dcast(data = x,
formula = formula(eval(paste(idvar, "~", groupvar))),
fun.aggregate = length, value.var = groupvar)
absDif <- abs((tabWide[, 2] / sum(tabWide[, 2])) - (tabWide[, 3] / sum(tabWide[, 3])))
hoover <- sum(absDif) / 2
return(hoover)
}
# simulation stochastique
listHoover <- numeric()
for(i in 1:100){
tabShuffle$GROUP <- sample(x = tabShuffle$GROUP,
size = length(tabShuffle$GROUP),
replace = FALSE)
tempHoover <- cast_and_compute_dissimilarity(x = tabShuffle, idvar = "GEOIDS10", groupvar = "GROUP")
listHoover <- append(listHoover, tempHoover)
}
hist(listHoover, breaks = 15, main = "Valeurs de Hoover simulées (n = 100)")
obsHoover <- sum(abs(sfCity$PROPDIFF / 100)) / 2
avgHoover <- mean(listHoover)
sdHoover <- sd(listHoover)
standHoover <- (obsHoover - avgHoover) / sdHoover
paste("Indice de Hoover standardisé (Cortese el al.) :", round(standHoover, digits = 2))
```
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment