Packages R utilisés

Premièrement durant mon stage, j’ai dû utiliser un nombre assez importants de packages et de fonctions, les voici :

spatstat : Spatstat est un package permettant d’analyser des semis de points. Ses fonctionnalités comportent notamment de l’analyse exploratoire de données, des méthodes de moédélisation et d’ajustement, de la simulation de processus ponctuels spatiaux ou encore des tests statistiques.

maptools : package permettant l’implémentation de fonctions pour la manipulation de données géographiques.

spdep : Spdep (Spatial dependence: weighting schemes, statistics and models) pour dépendance spatiale est un package apportant un ensemble de fonctions permettant de créer des objets matriciels de pondération spatiale.

factorextra : le package factorextra permet d’extraire et de visualiser les résultats des analyses de données multivariées (analyse en composantes principales, analyse factorielle multiple hiérarchique…).

fpc : Fpc (Flexible Procedures for Clustering) est un package permettant diverses méthodes d’agrégation et de validation d’agrégation (regroupement par régression linéaire,…)

Présentation des données

Dans cette analyse, nous étudierons la répartition spatiale de feux de forêts dans le bassin méditerranéen français en raison de son climat avec des étés chauds et secs et des hivers frais et humides et plus particulièrement ici dans le département du Var pour les années 2016 et 2017.

plot(Feux83_1617,cols=c("blue","red"),legend=TRUE)

 

FIGURE 1 : départs de feux dans le département du Var pour l’année 2016(X1) et l’année 2017(X2).

I. Analyse de la répartition spatiale d’un semis de points

L’objectif de cette partie est de présenter des méthodes graphiques, des tests statistiques de la concentration spatiale et des méthodes de détection d’agrégats sur un semis de point. Les tests statistiques sont souvent construits avec pour hypothèse H0 une répartition spatiale totalement aléatoire appelée CSR pour Complete Spatial Randomness. Le modèle se rapportant à cette hypothèse est le modèle de processus de Poisson homogène.

Un processus de Poisson d’intensité (\(\lambda\)) est un processus ponctuel homogène isotrope pour lequel la disposition des points est complètement aléatoire à chaque réalisation.

PARLEZ PLUS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Visualisations graphiques

Un semis de points est un ensemble de données donnant les emplacements spatiaux observés d’entités très diverses : personnes, animaux, événements, objets…

Nous nous intéresserons dans cette rubrique à différentes méthodes de visualisation graphiques des données nous permettant d’évaluer la structure spatiale de la répartition.

Considérons la représentation des occurences de feux de forêts suivante :

Feux83_1617$marks <- as.factor(Feux83_1617$marks)
Feux.2016 <- split(Feux83_1617)$X1
plot(Feux.2016,pch=1,legend=FALSE)

FIGURE 2 : départs de feux dans le département du Var seulement pour l’année 2016

Comptage du nombre des points dans des quadrats

Un moyen simple de vérifier l’hétérogénéité consiste à vérifier si des régions de même superficie contiennent un nombre à peu près égal de points.[1, page 168,8]
Nous verrons qu’un test statistique se rapporte à cette méthode statistique.

Q <- quadratcount(Feux83_1617, nx= 6, ny=4) 
plot(Q)
plot(Feux83_1617,pch=1,legend=FALSE,add=TRUE)


FIGURE 3 : Comptage du nombre points dans des quadrats pour les départs de feux du Var

Intensité

On peut définir l’intensité de la façon suivante : Soit Kh un noyau de bande passante h et x un point de R²,l’intensité lissée estimée en x définie par :

D’un point de vue pratique, le lissage spatial est une modélisation locale qui repose sur le choix de paramètres.

 

Le noyau décrit la façon dont le voisinage est appréhendé.
La fenêtre de lissage est le paramètre fondamental de l’analyse. Elle quantifie la «taille» du voisinage.

FIGURE 4 – Schématisation du lissage spatial avec un noyau

De plus si l’on soupçonne que l’intensité peut être hétérogène, elle peut être estimée de manière non paramétrique par des techniques d’estimation des noyaux ou encore par le comptage dans des quadrats. [8,12]

Intensité non paramétrique à noyau

On dénombre plusieurs techniques pour calculer les noyaux, cependant deux d’entre elles sont plus fréquemment utilisés :

Noyau gaussien :

Noyau quadratique :

Le noyau gaussien prend en compte l’ensemble des points de la zone d’étude tandis que le noyau quadratique donne un poids plus élevé aux points les plus proches qu’aux points éloignés.
De plus une des difficultés majeures est le choix de la fenêtre de lissage h, qui conditionne le niveau de lissage de l’intensité. Plusieurs méthodes existent bw.diggle, bw.ppl, bw.cvl, etc.Dans le reste de cette analyse nous privilégierons celle basée sur bw.CVL établie par Cronie et Van Lieshout. [8,12]

par(mfrow=c(1,2))
plot(density.ppp(Feux83_1617, sigma = bw.diggle(Feux83_1617),edge=T),
     main=paste("h =",round(bw.diggle(Feux83_1617),2)))  
plot(density.ppp(Feux83_1617, sigma = bw.CvL(Feux83_1617),edge=T),
     main=paste("h =",round(bw.CvL(Feux83_1617),1)))

FIGURE 5 – Différentes méthodes de lissage pour la visualisation

On peut ainsi voir la présence de zones plus denses que d’autres. A gauche la méthode de lissage basée sur une méthode de Diggle (Utilisation d’une technique de validation croisée pour sélectionner une fenêtre de lissage) et à droite une méthode de lissage fondé sur la formule de Campbell. On remarque ainsi bien la présence de zones plus denses que d’autre. Le choix le plus significatif est la méthode de droite où l’on peut voir de façon plus claire l’intensité de notre jeu de données. [1,8,12]

On peut de plus rajouter un contour afin de mieux visualiser les différentes intensités.

FV <- density(Feux83_1617)
plot(FV)
contour(FV, add=TRUE)

FIGURE 6 – Contour des intensités des départs de feux pour le département du Var

Espacement

Des informations supplémentaires sur un semis de point peuvent souvent être révélées en mesurant les distances ou les plus courtes distances dans le semis.

Diagramme de Stienen

 

Le diagramme de Stienen d’un semis de point est obtenu en traçant un cercle autour de chaque point de notre jeu de données, d’un diamètre égal à la distance de son plus proche voisin. [1, page 255 section Spacing]

stienen(Feux83_1617)

FIGURE 7 – Diagramme de Stienen pour les départs de feux pour le département du Var

Dirichlet et liste de voisinage

 

De la même façon on peut représenter un graphique de dirichlet en ajoutant une liste de voisinage basée sur des mesures entre points.[4]

FeuxV.sp <- dirichlet(Feux83_1617) 
FeuxV.sp <- as(FeuxV.sp, "SpatialPolygons")
FEUXVNeigh <- as.SpatialPoints.ppp(Feux83_1617)
FeuxVcoord <- FEUXVNeigh@coords
neighbours <- dnearneigh(FEUXVNeigh,d1=0,d2=4000)
plot(FeuxV.sp, border = 'lightgrey')
plot(neighbours, coordinates(FEUXVNeigh), add=TRUE, col='red')

FIGURE 8 – Graphique de Dirichlet et liste de voisinage basée sur le plus proche voisin pour les départs de feux dans le département du Var

Test pour la détection globale de répartition spatiale non totalement aléatoire

On s’intéresse ici à l’existence d’une hétérogénéité globale de la distribution spatiale. De nombreux tests dont les statistiques sont basées sur des distances existent : entre couples de points, au plus proche voisin, entre un point quelconque et le point au semis le plus proche. L’objectif de ces tests est de détecter une répartition spatiale non poissonienne sans localiser les zones où l’on s’éloigne le plus de l’hypothèse CSR. (Complete Spatial Randomness) Nous ne ferons pas une description exhaustive de ces tests, cependant en voici quelques-uns :

Adéquation à un modèle poissonien

Afin de savoir si un semis de point a une tendance aléatoire, répulsive ou agrégée on peut utiliser des test d’adéquation à un modèle poissonien.

Test de quadrat

Comme nous avons pu le voir précédemment, on peut diviser l’espace en quadrats. Par ce biais, on peut alors tester l’hypothèse d’adéquation à un modèle poissonien en utilisant un test de Chi².

La statistique du test de Chi² est la suivante :

Avec nj représentant le nombre de points observé dans un quadrat et mu la valeur moyenne de points sous l’hypothèse d’un modèle poissonien.

quadrat.test(Feux83_1617)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  Feux83_1617
## X2 = 61.888, df = 23, p-value = 4.058e-05
## alternative hypothesis: two.sided
## 
## Quadrats: 24 tiles (irregular windows)

On rejette ainsi l’hypothèse H0 qui est “la distribution de notre semis de point est aléatoire”, avec moins de 5% de chance de se tromper. Cependant on sait seulement que notre semis de point n’est pas distribué de façon aléatoire mais on ne sait actuellement pas si la distribution de notre semis de point est agrégée ou répulsive.

Test basé sur les distances

1. Loi sous H0 est connue (ou approchée)

 

Les tests qui suivent sont fondés dans le cas où la loi sous H0 est approchée.

Test de Clark-Evans

 

Le test de Clark-Evans se base sur la distance moyenne au plus proche voisin divisé par la distance moyenne attendue (sous l’hypothèse H0) au plus proche voisin. [1, page 259]

R suit une loi normale centrée réduite, et sa valeur peut être comparée à la valeur critique de la loi de Student, au seuil considéré (1,96 au seuil de 5%)

On note :

La distance attendue par rapport au voisin le plus proche dans le cas d’un processus de poisson d’intensité lambda.

clarkevans.test(Feux83_1617)
## 
##  Clark-Evans test
##  No edge correction
##  Z-test
## 
## data:  Feux83_1617
## R = 0.81329, p-value = 2.014e-06
## alternative hypothesis: two-sided

La présence du terme Z-test désigne le fait que la statistique de test suive une loi normale sous l’hypothèse nulle. De même en réalisant ce test et en obtenant cette p-valeur on ne sait pas encore si notre semis de point est répulsif ou s’il y a la présence d’un agrégat. Pour cela on peut préciser a notre test “alternative =”clustered"".

Test d’Hopkins et Skellman

 

Ce test est relativement proche du test de Clark-Evans, il se base cependant sur l’espace vide (i.e la distance d’un point de référence fixe de location u dans la fenêtre a son plus proche voisin.) et non plus sur la distance moyenne attendue.[1, page 259]

hopskel.test(Feux83_1617)
## 
##  Hopkins-Skellam test of CSR
##  using F distribution
## 
## data:  Feux83_1617
## A = 0.48203, p-value = 1.149e-11
## alternative hypothesis: two-sided

Le test d’Hopkins et Skellman compare la valeur A à la distribution de Fisher de paramètres (2m,2m).

1. Loi sous H0 est connue (ou approchée)

 

Dans ce cas on utilise des tests de Monte-Carlo, c’est-à-dire qu’on calcule la statistique de tests sur des échantillons simulés sous H0. Les simulations de Monte-Carlo testent l’hypothèse suivante: la valeur d’un indicateur est le résultat d’un processus aléatoire. Si l’hypothèse est rejetée, la valeur de l’indicateur est alors considérée comme étant significative.

Fonction G

 

Dans le cas de la fonction G on se place dans le cas d’une maille prédéfinie sur chaque point du semis. On construit la courbe G(r) de la fraction des points dont le plus proche voisin est à une distance inférieure ou égale à r. La définition de G est assez intuitive : en présence d’agrégats, on s’attend à voir les courbes augmenter rapidement ; au contraire, si les points se repoussent, peu de points auront leur premier voisin à petite distance. La fonction G décrit le voisinage des points. On réalise ici une envelope c’est-à-dire des simulations de la fonction, le nsim correspond ici au nombre de simulations soit 39 (39 simulations équivalent à un seuil de 5%.)

Genv <- envelope(Feux83_1617,Gest,nsim=39,fix.n=TRUE)
plot(Genv,main="Fonction du plus proche voisin")

FIGURE 9 – Utilisation de la fonction G pour caractériser un semis de point agrégé.

On rejette ici l’hypothèse H0 de CSR au seuil de 5%.

Fonction K de Ripley

 

La fonction K de Ripley fournit plus d’informations que la fonction G puisqu’elle prend en compte toutes les distances, donc tous les voisins et non seulement le plus proche. On définit K(r) comme une fonction étant égale au rapport du nombre de voisins à une distance <= r sur l’intensité dans le cas homogène. [2,3]

FIGURE 10 – Fonction K pour l’étude d’un semis de point

On peut conclure en rejetant l’hypothèse H0 de CSR au seuil de 5%.

Il existe encore de nombreuses fonctions basées sur des méthodes de distances différentes (fonction L, fonction F,etc…) [1, page 203 page 277, 2,3]

Autre alternative

On dénombre plusieurs autres tests basés notamment sur une covariable, la présence d’une covariable signifie que nous sommes plus dans le cas d’un test de CSR ou le processus de Poisson est homogène mais dans le cas d’un processus de Poisson inhomogène. [1 page 382] :

  • Kolmogorov-Smirnov

  • Anderson-Darling

  • Cramer-Von Mises

Généralisation à un autre type de modèle

On peut de même généraliser ces tests à un autre type de modèle, ainsi la qualité de l’ajustement d’un modèle statistique désigne le degré d’ajustement du modèle aux données observées. Pour tester l’hypothèse d’adéquation à un autre type de modèle on peut utiliser les tests suivants :

Test de Dao Genton

Supposons que l’hypothèse nulle soit un semis de point avec un paramètre thêta dont la valeur est inconnue et doit être estimée. Le problème essentiel du test de Monte Carlo est que la p-valeur p ne soit pas uniformément distribuée sous H0. Ainsi dans le test de Dao Genton la procédure habituelle du test de Monte Carlo appliquée aux données est également appliquée à chacun des modèles de points simulés [1 page 406]

Voici la p-valeur associé à ce test :

dg.test(Feux83_1617) #Trés long a lancer car 19 simulations Monte-Carlo (30-40 min)
## 
##  Dao-Genton adjusted goodness-of-fit test
##  based on Diggle-Cressie-Loosmore-Ford test of CSR
##  First stage: Monte Carlo test based on 19 simulations
##  Summary function: K(r)
##  Reference function: theoretical
##  Alternative: two.sided
##  Interval of distance values: [0, 22174.3180264945]
##  Test statistic: Integral of squared absolute deviation
##  Deviation = observed minus theoretical
##  Second stage: nested, 18 simulations for each first-stage simulation
## 
## data:  X
## p0 = 0.05, p-value = 0.05

Test de Diggle Cressie Loosmore Ford

On retrouve aussi le test de Diggle Cressie Loosmore Ford ou DCLF. Celui-ci se base sur la distance L. La statistique du test DLCF est ainsi la zone située entre les fonctions L empiriques et théoriques.

dclf.test(Feux83_1617,Lest,nsim=39,fix.n=TRUE,verbose=FALSE) #Long a lancer
## 
##  Diggle-Cressie-Loosmore-Ford test of CSR
##  Monte Carlo test based on 39 simulations with fixed number of points
##  Summary function: L(r)
##  Reference function: theoretical
##  Alternative: two.sided
##  Interval of distance values: [0, 22174.3180264945]
##  Test statistic: Integral of squared absolute deviation
##  Deviation = observed minus theoretical
## 
## data:  Feux83_1617
## u = 4.3642e+10, rank = 1, p-value = 0.025

Le type de test effectué dépend de l’argument X (le jeu de données) :

  • Si X est une sorte de modèle de points, alors un test CSR sera effectué. Autrement dit, l’hypothèse nulle est : “le modèle de points est complètement aléatoire”

  • Si X est un modèle de processus ponctuel ajusté, alors un test d’adéquation du modèle ajusté sera effectué. L’objet du modèle contient le modèle de points de données auquel il a été ajusté à l’origine. l’hypothèse nulle est “le modèle de points de données est une réalisation du modèle”.

Test pour la détection locale d’agrégation spatiale

Il existe quelques méthodes permettant d’estimer l’existence d’agrégats.

Test du rapport de vraisemblance

Le test du rapport de vraisemblance est optimal pour détecter une différence d’intensité à l’intérieur et à l’extérieur du cercle. Ces cercles sont définis de telle sorte qu’à chaque localisation spatiale u?, on dessine un cercle centré en u d’un rayon étant égal à r.

Dans notre cas voici ce qu’on obtient : [1, page 190]

On note nin le nombre de points dans chaque cercle de rayon b(u,r), nout le nombre de points à l’extérieur du cercle. Ain correspond ainsi à l’aire de la partie des cercles dans la région d’observation et Aout l’aire de la partie des cercles en dehors de la région d’observation.

par(mfrow=c(1,2))
LRTSCVL=scanLRTS(Feux83_1617,r=bw.CvL(Feux83_1617))
plot(pvals <- eval.im(pchisq(LRP,df=1,lower.tail = FALSE)))
X=pvals<0.05
plot(X,ncolours = 2)

FIGURE 11 – A gauche p-valeur du test du rapport de vraisemblance, à droite localisation de ses p-valeurs quand p<0.05

Détection d’agrégats

Dans cette section, nous nous intéressons à l’analyse de processus ponctuels, et de manière précise à la détection de zones de plus forte d’intensité des processus, appelées agrégats.

clusterset=clusterset(Feux83_1617,what = "domain")
plot(clusterset)

FIGURE 12 – Estimation non paramétrique (Allard-Fraley) de la région de haute intensité [1, page 192]

nnclean=nnclean(Feux.2016,k=7,plothist=TRUE)
plot(nnclean)

Le modèle de points qui en résulte est identique au jeu de données originel, à l’exception de deux colonnes de marques : un facteur qui classe le point en caractéristiques et en bruits. Le fait de rajouter plothist=TRUE voir première image donne la probabilité ajustée que le point appartienne à la composante d’agrégats du modèle de mélange. Le choix de k=7 étant arbitraire, le choix optimal de k dépend des caractéristiques du modèle,donc k est généralement choisi par tâtonnement.

FIGURE 13 – Histogramme des distances du septième plus proche voisin et modèle ajusté de densité du mélange.

FIGURE 14 – Nettoyage du voisin le plus proche pour les départs de feux de forêts dans le département du Var avec présence de région à haute intensité.

Partitionnement spatial

Le partitionnement spatial peut nous permettre d’identifier des zones avec des répartitions spatiales proches.

Algorithme de partitionnement

Le Clustering par partitionnement est une méthode de regroupement utilisée pour classer les observations, au sein d’un ensemble de données, en plusieurs groupes en fonction de leurs similarités.[5] Il y a plusieurs types de méthodes :

On verra ici l’algorithme des k-médoïdes. En statistiques, un médoïde est le représentant le plus central d’une classe.

L’algorithme PAM(Partitioning Around Medoids) est basé sur la recherche de k objets représentatifs ou de médoïdes parmi les observations de l’ensemble des données.

FeuxVar <- as.data.frame.ppp(Feux83_1617)
FeuxV <- FeuxVar
FeuxV <- na.omit(FeuxV)
FeuxV=FeuxV[,c(-3)]
FeuxV=scale(FeuxV)

fviz_nbclust(FeuxV,kmeans,
             method = "silhouette")
pam.res <- pam(FeuxV, 6)
pm <- cbind(FeuxVar, cluster = pam.res$cluster)
fviz_cluster(pam.res, data = FeuxV,
             repel=FALSE,
             geom="point",
             ellipse.type = "convex",
             pch=19,
             alpha=0.5,
             ggtheme = theme_minimal())

FIGURE 15 – Estimation du nombre optimal d’agrégats basée sur la méthode de la silhouette moyenne

FIGURE 16 – Visualisation des agrégats PAM

Il existe d’autres algorithmes de partitionnement comme la méthode basée les k-means très similaire aux k-médoïdes. Cependant ces méthodes possèdent quelques désavantages :

  • K-Means ne forme que des amas sphériques. Cet algorithme échoue lorsque les données ne sont pas sphériques (c’est-à-dire que la variance est la même dans toutes les directions).

  • L’algorithme de K-Means est sensible aux valeurs extrêmes. Les valeurs extrêmes peuvent fausser les agrégats de K-Means dans une très large mesure.

  • L’algorithme de K-Means exige que l’on spécifie le nombre de clusters en amont.

FIGURE 17 – Visualisation d’un algorithme basé sur les K-means

Une alternative à cette méthode est notamment une méthode basée sur la densité, que l’on verra juste après dans cette section.

Agrégation hiérarchique

La méthode de classification hiérarchique est une approche alternative pour regrouper les objets en fonction de leur proximité. La classification hiérarchique n’exige pas de prédéfinir le nombre de clusters à produire, il faut alors les choisir.

res.hc <- hclust(dist(coords(Feux83_1617)), method = "ward.D2")
fviz_dend(res.hc,k=4,color_labels_by_k = TRUE,rect = TRUE,show_labels = FALSE) #k correspond au choix du nombre de cluster

FIGURE 18 - Dendrogramme basé sur une méthode de Ward

Méthode avancée d’agrégation

Les méthodes d’algorithme de partitionnement (“K-means”, “PAM clustering”) et d’agrégation hiérarchique fonctionnent pour trouver des agrégats de formes sphériques ou convexes. En d’autres termes, elles ne conviennent qu’aux agrégats compacts et bien séparés, elles sont cependant très affectées par la présence de bruits et de valeurs aberrantes dans les données.

DBSCAN (Density-Based Spatial Clustering and Application with Noise), est un algorithme de clustering basé sur la densité qui peut être utilisé pour identifier des agrégats de toutes formes dans un ensemble de données contenant du bruit et des valeurs extrêmes. [6,7]

FIGURE 19 - Résultats d’agrégation de trois bases de données différentes avec DBSCAN

Les deux paramètres essentiels pour un algorithme DBSCAN sont les suivants :

  • eps : Il définit le voisinage autour d’un point de données, c’est-à-dire que si la distance entre deux points est inférieure ou égale à “eps”, ils sont considérés comme voisins. Si la valeur eps est choisie trop petite, une grande partie des données sera considérée comme aberrante. Si elle est choisie très grande, les agrégats fusionneront et la majorité des points de données se trouveront dans les mêmes agrégats. (Une façon de trouver la valeur eps est basée sur le graphique de la distance k).

  • MinPts : Nombre minimum de voisins (points de données) dans un rayon de eps. Plus l’ensemble de données est grand, plus la valeur de MinPts doit être importante. La valeur minimale de MinPts choisit doit être d’au moins 3.

db = dbscan(FeuxVarPoints, eps = 6000,minPts =  5) 
hullplot(FeuxVarPoints,db$cluster)
plot(window, add=T)

FIGURE 20 - Agrégation basée sur la méthode DBSCAN pour les départs de feux du département du Var

L’algorithme DBSCAN produit un résultat assez raisonnable puisqu’on distingue la présence de 6 agrégats. Ces zones distinctes représentent des endroits où la répartition des départs de feux est plus importante.

II. Comparaison de la répartition spatiale dans deux semis de point

Jusqu’à présent, notre analyse se basait sur un seul semis. Mais dans certains cas, les analyses de semis de point impliquent des relations entre plus d’un semis de point. On s’intéresse ici à la comparaison de la répartition spatiale dans 2 semis de point. [11]

Visualisations graphiques

Pour la comparaison de deux semis de point il n’existe pas de méthodes de visualisation bien distinctes d’un semis de point.

plot(Feux83_1617,cols=c("blue","red"),legend=TRUE)

Voir Figure 1. On peut néanmoins représenter les départs de feux dans le département du Var pour l’année 2016 et l’année 2017 en précisant legend=TRUE.

Test pour la détection globale d’agrégation spatiale

De la même façon on peut détecter de façon globale une agrégation spatiale. Pour cela il existe plusieurs techniques.

Statistique de balayage

Une statistique de balayage est la valeur maximale d’un indice de clustering mesurée sur un ensemble de clusters candidats. Comme de nombreuses méthodes de balayage, elle s’appuie sur un rapport de vraisemblance mais elle prend également en compte les corrélations entre les différentes variables. L’idée que nous étudions ici est l’adaptation des statistiques de balayage des tests locaux à des tests globaux en considérant tous les indices de clustering mesurés sur les clusters potentiels en ne retenant que le maximum.

Feux83_1617$marks <- as.factor(Feux83_1617$marks)
Feux.2016 <- split(Feux83_1617)$X1
Feux.2017 <- split(Feux83_1617)$X2
resscan=signscan(Feux.2016,Feux.2017,19)
resscan$pvalvar
## [1] 0.1

Test de ségrégation

L’hypothèse nulle est l’hypothèse plus concrète de l’étiquetage aléatoire : étant donné l’emplacement x, les types d’étiquettes mi sont des variables aléatoires, indépendantes les unes des autres, avec une distribution de probabilité commune P{mi=m}=pm. Cette hypothèse nulle implique que p(m|u) = pm soit constant en fonction de u, pour chaque m, de sorte que les types ne soient pas séparés. Voici la formule liée à ce test statistique : [1 page 580]

Avec phat(m | x[i]) l’estimateur de la probabilité que le i-ème point de données soit de type m, et pbar[m] est la fraction moyenne des points de données qui sont de type m. La statistique T est évaluée pour les données et pour les versions randomisées nsim de X, générées par permutation ou ré-échantillonnage aléatoire des marques.

Feux83_1617$marks <- as.factor(Feux83_1617$marks)
seg=segregation.test(Feux83_1617,nsim = 39) #Nombre de simulation pour un test de MonteCarlo
## 
##  Monte Carlo test of spatial segregation of types
## 
## data:  Feux83_1617
## T = 0.082162, p-value = 0.225

On rejette ici l’hypothèse nulle “le modèle originel et sa version randomisée sont interchangeables, ce qui justifie le test de Monte Carlo” au seuil de 5%.

Une alternative possible est disponible dans le package spatialsegregation écrit par T.Rajala [1 page 581]

Étiquetage aléatoire

Un test d’étiquetage aléatoire est mieux adapté lorsque l’étiquette type de chaque point représente son statut, comme par exemple le statut de la maladie (sain/malade). L’hypothèse nulle d’un étiquetage totalement aléatoire est que le statut de chaque point est déterminé au hasard, indépendamment des autres points.

On peut utiliser un test basé sur les enveloppes, la fonction “dot” est celle qui convient le plus pour évaluer l’étiquetage aléatoire. Dans l’hypothèse d’un processus ponctuel de plusieurs types stationnaires, l’étiquetage aléatoire complet implique que Ji(r)=J(r) où J est la fonction pour le processus ponctuel sans marque. On pourrait d’ailleurs utiliser la fonction K ou la fonction G. Bien que ces choix ne soient pas strictement nécessaires pour la validité du test, ils sont utiles pour interpréter les écarts par rapport à l’hypothèse nulle.

Jdif <- function (X,...,i){
  Jidot <- Jdot(X,...,i=i)
  J <- Jest(X,...)
  dif <- eval.fv(Jidot-J)
  return(dif)
}
EPJ <- envelope(Feux83_1617,Jdif,nsim=39,i="X1",simulate = expression(rlabel(Feux83_1617)))
plot(EPJ)

FIGURE 21 - Test d’échantillon aléatoire pour le jeu de données des départs de feux du département du Var

Ce graphique représente ici une enveloppe basé sur 39 simulations d’étiquetage aléatoires avec i=X1 (Départs de feux pour l’année 2016). Les résultats suggèrent que les marques pour ce jeu de données sont aléatoires.

Test pour la détection locale d’agrégation spatiale

Pour la détection locale d’agrégation spatiale dans le cas de deux semis de point, le nombre de tests statistiques est assez faible, on peut néanmoins citer les suivants :

Risque relatif local

Le risque relatif (RR) est une mesure statistique souvent utilisée en épidémiologie, mesurant le risque de survenue d’un événement dans un groupe par rapport à l’autre.

La fonction de risque relatif est définie ici comme le rapport entre la densité de “cas” et le “contrôle”. En utilisant l’estimation de la densité du noyau pour modéliser ces densités, nous obtenons une estimation exploitable de celles-ci. Cette fonction définit la fonction de risque r de la manière suivante :

Où Fd et Gd désignent respectivement les estimations de la densité de cas et de contrôle.

rr = relrisk(Feux83_1617, 6000)
plot(rr)
plot(eval.im(rr < 0.15),col=c("white","lightsteelblue"))
plot(window,add=TRUE)
plot(Feux83_1617,cols=c("blue","red"),legend=TRUE,add=TRUE)

FIGURE 22 - Localisation des p-valeurs du test du risque relatif quand p<0.05

Alternative

Le package sppt (Spatial Point Pattern Test)

Les tests de ce package mesurent le degré de similarités au niveau local entre deux modèles de points spatiaux et sont basés sur une zone. Ils ne sont pas utilisés pour tester des modèles de points avec les hypothèses nulles de distributions aléatoires, uniformes ou agrégées, mais peuvent être utilisés pour comparer un modèle de point particulier avec ces distributions. L’un des avantages de ces tests est qu’ils peuvent être effectués pour un certain nombre de limites de zones différentes en utilisant les mêmes ensembles de données ponctuelles d’origine.

Détection d’agrégats

Pour la détection d’agrégats on peut retrouver la fonction statscan de balayage définie plus haut. [10]

plot(Feux83_1617,main="",cols=c("blue","red"),legend=FALSE)
symbols(resscan$center[1],resscan$center[2],circles=resscan$radius,inches=F,fg='blue',add=TRUE)

FIGURE 23 - Localisation des occurrences de départs de feux dans le Var

Références

  1. Adrian Baddely-Ege Rubak - Rolf Turner (2016), “Spatial Point Patterns Methodology and Applications with R”, Chapman and Hall/crc.

  2. Eric Marcon. Statistiques spatiales avec applications à l’écologie et à l’économie. Biodiversité et Ecologie. AgroParisTech, 2010. Français.

  3. Analyse de la structure spatiale de semis de points hétérogènes : exemples d’application à des peuplements forestiers

  4. https://r.developpez.com/tutoriels/programmation-graphe/livre-R-et-espace/?page=chapitre-11-initiation-aux-statistiques-spatiales

  5. https://www.datanovia.com/en/blog/cluster-analysis-in-r-practical-guide/

  6. https://rstudio-pubs-static.s3.amazonaws.com/126356_ef7961b3ac164cd080982bc743b9777e.html

  7. https://www.geeksforgeeks.org/dbscan-clustering-in-ml-density-based-clustering/

  8. https://mgimond.github.io/Spatial/point-pattern-analysis.html

  9. http://www.50northspatial.org/clark-evans-test-r/

  10. Global scan methods for comparing two spatial point processes

  11. https://www.seas.upenn.edu/~ese502/NOTEBOOK/Part_I/5_Comparative_Analyses.pdf

  12. Lissage spatial - LAURE GENEBES, AURIANE RENAUD ET FRANÇOIS SÉMÉCURBE,Insee

  13. http://www.sthda.com/english/wiki/wiki.php?id_contents=7940