QGIS: des Ellipses de Déviation Standard (SDE), un plugin, « Standard Deviational Ellipse », des scripts R (processing) et Python et une approche critique...

	Niveau	Intermédiaire
	Logiciels utilisés	QGIS R Python QGIS Plugin: Standard Deviational Ellipse
	Plateforme	Windows | Mac | Linux | FreeBSD

Parmi les commandes d’ArcGIS qui manquaient dans QGIS, figurait l’ Ellipse de Déviation Standard (Standard Deviational Ellipse ou SDE). Ce n’est plus le cas actuellement puisque c’était déjà possible de le faire avec des scripts R dans la Boîte à outils de Traitements (Processing) avec des packages dédiés et le nouveau plugin Standard Deviational Ellipse. Ses premiers résultats ont suscité des réactions de ma part (problem with the rsd values compared with the aspace or phonTools R packages ) qui ont amené Håvard Tveite à améliorer et compléter son plugin. 

Je vais essayer ici d’expliquer le pourquoi des choses, de manière visuelle et sans trop de statistiques, car ce concept a sérieusement été remis en cause récemment par Gong, J (2002), Clarifying the standard deviational ellipse. Geographical Analysis 34, 155–167.(pdf) et Wang B, Shi W, Miao Z (2015) Confidence Analysis of Standard Deviational Ellipse and Its Extension into Higher Dimensional Euclidean Space. PLoS ONE 10(3): e0118537. doi:10.1371/journal.pone.0118537,

C’est quoi ?

Si la dispersion d’une variable en une dimension est facile à calculer et à visualiser (histogrammes, etc.), c’est moins le cas pour des distributions bivariées (ou multivariées). La solution classique est de calculer séparément les écarts types de chaque variable et de représenter leur distribution par des ellipses (ou des ellipsoïdes) dont les axes sont ces écarts types. La méthode permet aussi de mesurer une tendance directionnelle (angle de l’ellipse).

Dans la réalité ce sont en fait des isocontours d’une surface, similaires à des courbes de niveau.

Et le graphique réalisé avec le package R psych  qui illustre le fait que cette démarche n'est pas limitée au domaine géospatial, loin de là.

Cette ellipse est donc théoriquement déterminée par trois mesures :

1. calcul du point « moyen » du nuage des points (cela peut être le centroïde, le barycentre ...)
2. calcul de l’angle de rotation de l’ellipse (orientation)
3. calcul de la dispersion par les écarts-types des coordonnées x et y à partir du  point « moyen ». Ceux-ci vont définir les axes de l’ellipse (c’est pourquoi ESRI nomme aussi l’ellipse  « ellipse de l’écart type ».

Toutes les combinaisons sont possibles: affecter une masse (valeur) à chaque point, utiliser des moyennes pondérées, des médianes ou la covariance des données (pondérée ou non). Le calcul peut aussi varier en fonction de la répartition des points et il est possible d’utiliser soit la distance euclidienne soit d'autres distances ( comme la Distance de Mahalanobis ou la Distance de Manhattan).        

Je n’illustrerai ici que le cas de valeurs non pondérées (le reste n'est qu'un facteur d'échelle et de rotation).

Statistisquement parlant, l'objectif recherché est d'obtenir les valeurs de  la classique Règle 68-95-99.7 (3 sigma) empirique ( graphique réalisé avec Python et le module matplotlib) dans le cas d'une distribution univariée, et pas bivariée (Wang et al. (2015) contestent ces valeurs et en proposent des nouvelles pour les distributions multivariées, voir plus bas)

Un peu d’histoire ou comment une ellipse n’est pas une ellipse, mais...

Proposée dès 1926 par Lefever (Measuring Geographic Concentration by Means of the Standard Deviational Ellipse. American Journal of Sociology. 1926; 32(1):88–94), il a vite été démontré que le résultat n’était pas une ellipse (Furfey PH. A Note on Lefever’s "Standard Deviational Ellipse". American Journal of Sociology. 1927; 33(1):94–8), mais à ce qui est nommé actuellement une « Courbe de Déviation Standard ».  Depuis lors, bien évidemment, de nombreuses corrections ont été apportées pour tenter de résoudre le problème avec des résultats variables, mais la controverse continue... 

Il faut noter que le concept est né sous l’hypothèse que les variables observées suivaient la loi normale, mais que cette condition n’est plus indispensable avec les algorithmes utilisés actuellement. Elle est théoriquement dérivée de la distribution bivariée (ou bidimensionnelle)

A l’heure actuelle, dans le monde geospatial, ne subsistent que deux grands algorithmes et des solutions intermédiaires:

• celui proposé par Yuill, R. S.(1971) dans The Standard Deviational Ellipse: An Updated Tool for Spatial Description, Geografiska Annaler 53B(1),28-39) qui est une énième adaptation réussie de la solution de Lefever, avec de nombreuses variantes postérieures;
• celui implanté  dans le  logiciel CrimeStat III (Ned Levine, 2010, A Spatial Statistics Program for the Analysis of Crime Incident Locations (version 3.3).  Ned Levine & Associates, Houston, TX.; National Institute of Justice, Washington, DC.), spécialement dans la note 4 de Part II: Spatial Description (pdf).

Wang et al. (2015) ont récemment proposé un nouvel algorithme basé sur la décomposition spectrale de la covariance, plus exact selon eux.

Pourtant les géomaticiens continuent à utiliser cette méthode dans de nombreux domaines, sans se préoccuper ni de l’algorithme utilisé, ni de la signification statistique réelle du résultat.

Préambules statistiques: l’éllipse de covariance

Au départ, pour toute distribution bivariée, il est possible de calculer l’ellipse de covariance (ou ellipse de dispersion, ou ellipse des erreurs, ou...). Une matrice de covariance peut en effet être décomposée en (voir A geometric interpretation of the covariance matrix, par example):

• vecteurs propres qui permettent de calculer les 2 angles de l’ellipse. Le premier vecteur (violet) est la direction dans laquelle les données varient le plus, le second (bleu) est orthogonal au premier.
• valeurs propres qui sont les 2 diamètres de l’ellipse. Les rayons sont les racines carrées des valeurs propres et correspondent aux écarts types modifiés (translation + rotation) de x et y, σx' et σy’

La courbe limite de l’ellipse (Δχ2) correspond à un contour de probabilité de 1 écart type avec une probabilité de ± 39,35 % (c'est le genre de phrase qui nécessiterait de très nombreuses explications qu’il n’est pas possible de développer ici et je vous renvoie à vos manuels de statistiques).

Un script Python qui permet de le faire est disponible à par_covar.py. Le même résultat peut être obtenu avec les packages R ellipse:ellipse, siar: standard.ellipse , PhonTools: sdellipse, mvdalab: ellipse.mvdalab et autres.

Résultats:

En multipliant les valeurs obtenues par des facteurs d’échelle appropriés (loi normale, voir How to draw a covariance error ellipse?), il est donc possible de calculer des ellipses d’équiprobabilités. Ceci sera fait avec le package R car: dataEllipse() qui permet de superposer les contours d'une distribution bivariée normale de deux variables sur un semis de points. Un script R pour tracer ces ellipses dans la Boîte à Outils est disponible à car_ellipses.rsx.

Et il est possible de constater que l’ellipse de covariance correspond bien à une probabilité de ± 39,35 %.

Ces ellipses vont nous servir de repère pour superposer les différentes ellipses calculées dans la suite. Elles ne sont en effet que des adaptions et c'est pourquoi les coordonnées du point « moyen »,  les angles et l’excentricité restent les mêmes pour toutes les solutions.

L'interface du plugin

L’ellipse de Yuill

Les calculs d’estimation des rayons/écarts types (standard deviation) sont:

Elle est calculée avec l'option Yull du plugin.

Résultats (pour rappel, les valeurs des angles et de l’excentricité sont les mêmes que celles obtenues pour l’ellipse de covariance): 

 C’est aussi l’ellipse calculée par: 

• l’ancien script en Avenue sellipse.ave du livre Statistical Analysis with ArcView GIS;
• les versions d’ArcGIS inférieures à la 10.x (ArcMap 9.3: How Directional Distribution: Standard Deviational Ellipse (Spatial Statistics) works).

L’ellipse de CrimeStats III

Ned Levine, dans le cas de CrimeStats III, démontre dans la note 4 (voir référence plus haut) qu’avec les traitements précédents, l’ellipse résultante est trop petite (d’un facteur racine carrée de 2) pour deux raisons statistiques, le calcul des axes et le fait que les autres solutions ne sont pas correctes au niveau des degrés de liberté (n-2 au lieu de n). Les résultats sont beaucoup plus grands qu’avec l’algorithme de Yull.

Elle peut être calculée avec le package R aspace: calc_sde()  (un script R pour la Boîte à Outils de Traitement (Processing) est disponible à aspace_SDE.rsx ) ou avec l'option CrimeStat du plugin.

Résultats: 

La solution d’ArcGIS  10.x

Entre la version 9.x et la version 10.x, ArcGIS a modifié son algorithme en utilisant une des corrections de CrimeStats (racine carrée de 2) sans autre explication ni justification (ArcMap 10.x: Fonctionnement de Directional Distribution (Standard Deviational Ellipse))

Elle est calculée avec les options Yuill + sqrt(2) correction  du plugin et avec divers packages R (ellipse pour 1 écart type ici)

Résultats:

Le plugin offre aussi la possibilité d’utiliser la deuxième correction de CrimeStats seule avec l’option Yull + DF correction.
 

L'ellispse de Wang et al. (2015)

Le  nouvel algorithme proposé par Wang et al (2015) se base aussi sur la matrice de covariance, mais en modifiant le cacul  (1/n au lieu de 1/n-1 pour les connaisseurs). Ensuite le traitement est le même, voir le script ellipse_wang.py.

Résultats (c'est le même que celui de Yuill): 

L'important est qu'ils calculent les valeurs de  la classique règle empirique des 3 sigmas pour des distributions multivariées. Ce sont les valeurs approximatives de 39,35%, 84,47% et 98,89% et non 68,27%, 95,45% et 99.73% comme le stipule ESRI.

Conclusions

Que d’ellipses, me direz-vous, mais elles ne se différencient que par les résultats du calcul des axes de l’ellipse et donc de l’aire. Du  fait de la signification statistique des résultats, c’est là où naissent les problèmes et les controverses.

• dans le cas d’ESRI, il est stipulé que « un polygone d’ellipse d’un écart type recouvre approximativement 68 pour cent des entités» (que ce soit avec les résultats de la version 9.x ou ceux la version 10.x) , en supposant que variables concernées suivent une distribution spatiale normale;
• Cette conclusion est remise en question par Wang et all (2015) qui opinent que les conclusions d’ESRI correspondent à la classique Règle 68-95-99.7 (3 sigma) des distributions univariées normales, mais pas pour des  distributions bivariées (ou multivariées). Ils proposent donc des nouvelles valeurs: "Obviously, such confusing interpretation may mislead the GIS users to believe the univariate 3-sigma rule remains valid in two-dimensional Euclidean space, or even higher dimensions". Ils calculent des valeurs correctes selon eux;
• la solution de CrimeStat s'éloigne des précédentes;
• les ellipses basées sur la seule matrice de covariance sont théoriquement les plus correctes.

Alors quelle est la bonne ellipse ? Le nouveau plugin permet de toutes les figurer (hormis les eclipses de covariance et les axes) et à vous choisir la bonne...

Pour ceux qui ne s’intéressent qu’aux orientations (angles et excentricité), comme moi, elles sont toutes bonnes, mais pour ceux qui traitent les surfaces et leurs interprétations statistiques, c’est autre chose... Il serait nécessaire en tout cas de toujours spécifier la solution retenue comme dans la majorités des articles scientifiques.

Et oui, il n'y a donc pas une seule Elipse de Deviation Standard comme pourraient le penser les utisateurs d'ArcGIS avec leur commande standardisée. Je remercie aussi Håvard Tveite d'avoir tenu compte de mes remarques (j'utilise depuis longtemps ce genre d'ellipse dans un contexte non spatial).

En route pour les ellipses de covariance 3D et les ellipses de Deviation Standard bayesiennes.

Annexe:

Sommaire des résultats obtenus avec les points utilisés

Pour rappel, les coordonnées du point « moyen », les angles et l'excentricité sont les mêmes pour toutes les solutions (77,12° par rapport au N et 0,9144).

Tous les traitements présentés ont été effectués sur Mac OS X ou Linux avec QGIS 2.14.x, R version 3.x et Python 2.7.x. ArcGIS 10.1 sur Windows a été utilisé pour obtenir des ellipses pour comparer les résultats. LaTeX a été utilisé pour les équations.

licence Creative Commons Paternité-Pas d'Utilisation Commerciale-Pas de Modification 2.0 France