Tag Archives: 3wen

Des modeles prédictifs en assurance

Cet post est aussi en ligne sur https://hal.archives-ouvertes.fr/hal-02350006, il a été coécrit avec Laurence Barry et Ewen Gallic.

Les compagnies d’assurance émettent des contrats qui prévoient des paiements d’indemnités en cas de survenance d’évènements aléatoires (accident, maladie, décès, etc.). En contrepartie, l’assuré doit s’acquitter d’une prime, dont le montant est déterminé ex-ante, avant le début de la période de couverture. Cette prime se décompose en deux termes : une prime pure (destinée à couvrir les pertes anticipées) et un chargement (incluant des commissions à des agents, divers frais, mais aussi couvrant contre le risque de variabilité des pertes). La prime pure est souvent calculée par classe de risque, et une classification est alors nécessaire.

Assurer une population hétérogène, ou l’importance de la classification

Le regroupement des risques selon diverses informations telles l’âge de l’assuré, son état de santé ou encore sa profession constitue ce que l’on appelle la classification des risques. Cette pratique de segmentation se justifie (à des fins d’admissibilité mais aussi de tarification) par la supposition que les risques sont placés dans des groupes relativement homogènes, au sein desquels les probabilités de survenance sont similaires. Pour Schauer (2006), cette « généralisation », qui vise à voir l’individu sous le prisme de sa classe de risque, de généraliser son comportement à partir de quelques variables explicatives, est probablement la raison d’être de l’actuaire : « to be an actuary is to be a specialist  in  generalization,  and  actuaries  engage  in  a form of decision-making that is sometimes called actuarial ». Statistiquement on cherche une méthode de classification aussi « discriminatoire » que possible[1], en gardant en mémoire que la discrimination est interdite, ce qui rend l’exercice périlleux, et souvent critiqué (nous y reviendrons plus loin).

Les assureurs évoquent souvent deux arguments pour justifier une segmentation. Le premier est qu’elle serait rendue économiquement nécessaire par la concurrence ; ne pas classifier conduit à une anti-sélection, les risques importants restant seuls chez les assureurs qui ne segmentent pas. Dans une telle situation, l’équilibre de marché ne serait pas possible puisque les risques faibles seraient chez un concurrent ayant segmenté. Si le facteur de risque était observable, tant par les assurés que les assureurs, il y aurait un phénomène d’auto-sélection, les assurés à risque faible ayant les polices les moins chères. Cette situation constitue un équilibre séparant de Nash. Mais si le facteur de risque n’est pas observable, un équilibre sous-optimal peut être atteint, résultant d’une externalité négative de cette information non-accessible, à la manière de Wilson (1977), tel que décrit dans Cummins et al. (1982) dans le cas des contrats d’assurance-vie. Cela dit, Kleindorfer & Kunreuther (1980) montrent qu’accéder à davantage d’information ne conduit pas nécessairement à une amélioration du bien-être des consommateurs. De plus si la classification n’est pas autorisée, l’équilibre est maintenu, les risques faibles subventionnant les risques élevés.

Le second argument avancé pour justifier une segmentation est que cette dernière (et le fait, par conséquent, d’ajuster les primes au risque) serait juste et équitable. Mais cette vision de l’équité n’a pas toujours été de mise et semble portée par les développements techniques. Ainsi la classification est devenue de plus en plus fine, multipliant les classes de risque et conduisant à des tarifs « personnalisés ». En plus des avancées statistiques, des facteurs économiques pourraient justifier cette sophistication : la concurrence de plus en plus forte sur certaines branches.

Incertitude en assurance

Il y a plusieurs manières de caractériser l’incertitude en assurance. Comme bien souvent quand on fait des prévisions, il convient de distinguer l’incertitude associée à l’estimation des probabilités et l’incertitude réelle sur le résultat (aléa de l’évènement). Pour la seconde notion, Hacking (1975) parle de probabilité structurelle, et c’est celle qui est souvent utilisée pour introduire les concepts de probabilité, par exemple avec des dés ou des jeux de cartes : les probabilités sont connues, seule l’issue du jeu est incertaine. Par exemple je sais que la probabilité d’avoir 6 en lançant un dé est 1/6 (compte tenu de la géométrie du cube).

D’un point de vue statistique, la probabilité se mesure quand on peut observer une fréquence, c’est-à-dire une répétition de risques semblables. Les statisticiens ont ainsi défini une notion de probabilité empirique, basée sur la répétition[2]. Si, en lançant mille dés j’obtiens 173 fois la face 6, la probabilité empirique d’avoir 6 est de 17,3%. La loi des grands nombres nous assure que cette fréquence va tendre vers la vraie valeur en répétant l’expérience, et le théorème central limite permet d’en contrôler les fluctuations. C’est la première incertitude dont nous parlions au début de cette section, que nous appellerions l’erreur d’estimation.

On peut enfin mentionner deux notions supplémentaires ; tout d’abord, les probabilités conditionnelles. Cette idée est introduite en assurance par de Moivre, ou de Witt, lorsqu’ils notaient que pour estimer une probabilité de décès, il fallait considérer des personnes de même âge. C’est cette idée que l’on retrouve quand on considère une classification : on veut des risques homogènes, similaires, sans être pour autant identiques. La probabilité que l’on obtient est alors conditionnelle à ce facteur commun qui caractérise la classe observée. Dans notre exemple des dés, cela revient à dire qu’il ne faut pas lancer mille dés, mais mille fois le même dé – ou à défaut des dés semblables.

Enfin, les probabilités subjectives ont été formalisées par Bruno de Finetti et Leonard Savage (ainsi que plus philosophiquement par Frank P. Ramsey) pour comprendre et modéliser la prise de décision. Elles sont relativement populaires en économie de l’incertain, mais difficile à mettre en œuvre dans un contexte de valorisation de contrats d’assurance automobile ou habitation. Il s’agit d’un jugement, qui ne peut être confronté à la réalité, mais envisageable pour l’assurance de risques encore mal connus (McGrayne (2012) évoque ainsi les premiers contrats d’assurance aviation). Une approche bayésienne consiste alors à combiner cette probabilité subjective avec la probabilité comme fréquence observée d’un phénomène : partant d’une croyance a priori, on affine l’estimation par une mise à jour progressive en répétant les expériences. Classiquement, la probabilité d’avoir la face 6 sera une moyenne entre notre croyance (1 chance sur 6) et une probabilité dite historique, obtenue en faisant quelques lancés (3 sur 20 lancers, par exemple). Les poids attribués aux deux dépendant du nombre d’expériences effectuées : on donnera plus de crédit à l’expérience si on fait mille lancés que si on en fait soixante.

Incertitude sur le résultat, ou aléa fondamental

Les probabilités prédictives, utilisées pour calculer la prime d’un contrat d’assurance, sont la première étape d’un problème de classification. Un outil classique pour juger de la pertinence d’un classifieur est la courbe ROC, décrite dans Kuhn (2018)) : on compare la probabilité individuelle (a priori, telle que résultant du modèle de classification) à un seuil, compris entre 0 et 1; si la probabilité est inférieure au seuil, l’estimation est que la personne survit, sinon qu’elle décède.

On compare ensuite cette estimation aux réalisations (ex-post) de survie et de décès. Pour chaque seuil, on peut considérer la matrice classique dite matrice de confusion de théorie de la décision : elle consiste à répartir les observations suivant le résultat observé (en colonne) et l’estimation résultant du modèle en ligne (en fonction de la probabilité estimée pour l’individu et le seuil que l’on s’est fixé). On peut ainsi partager la population entre les classements corrects, et les erreurs (dont les « faux positifs » si la personne a survécu malgré une probabilité estimée de décès supérieure au seuil, et les « faux négatifs » si la personne décède malgré une probabilité estimée inférieure au seuil).

Figure 1: Courbe ROC et classification pour un seuil de probabilité valant 1.5%.

La courbe ROC est obtenue en faisant varier le seuil. Chaque seuil correspond à un point de la courbe, rapportant graphiquement les taux de faux positifs (en abscisse) et de vrais positifs (en ordonnée), comme sur la Figure 1.

Considérons un groupe de 1000 assurés, où 20 personnes sont décédées l’an passé. Supposons un modèle dans lequel on admet que la population est parfaitement homogène, la probabilité estimée de décès est de 2% pour tout le monde. Dans ce cas pour tout seuil supérieur à 2%, on estimera que la totalité de la population survit : on aura un taux de faux positifs de 0% et un taux de vrais positifs de 0%, d’où un point (0,0) sur le graphe. A l’inverse pour tout seuil inférieur à 2%, on estimera que la totalité de la population décède : on aura un taux de faux positifs de 100% et un taux de vrais positifs de 100%, d’où un point (1,1) sur le graphe. La courbe de ROC de ce modèle uniforme à 2% est donc la diagonale du carré sur la figure 1.

Mais on peut aussi imaginer qu’il existe un peu d’hétérogénéité avec, par exemple, une probabilité de décès de 1% pour une moitié de la population et de 3% pour l’autre moitié, ou encore que le modèle produit des probabilités comprises entre 1% et 3% de façon non dichotomique. Les données simulées pour construire la courbe noire sur la Figure 1 suppose que la population a des probabilités de décès variables, comprises entre 1% et 3%, obtenues par une régression logistique. Comme le montre le tableau de droite, on commet des erreurs, et comme le montre celle de gauche, la nature de celle-ci varie en fonction du seuil choisi, qui modifie les taux de faux positifs et de faux négatifs.

Le cas extrême serait celui où le modèle aurait correctement attribué une probabilité de 100% aux 20 personnes qui sont effectivement décédées. C’est la courbe rouge sur la Figure 1. Ce partage est possible ex-post, une fois réalisation de l’aléa : a posteriori, il y a une certitude de décès pour ceux qui sont effectivement morts. Mais cela n’a cependant pas grande réalité dans l’assurance, à moins d’imaginer que l’actuaire serait un oracle, qui saurait avec certitude qui va mourir, et qui va survivre. La réalité est plutôt celle de la situation intermédiaire entre la courbe rouge et la diagonale, avant d’arriver dans la région hachurée, où le taux d’erreur est faible, mais pas nul : on ne peut pas prédire, avec certitude, qui va décéder. L’assurance n’est possible que si cette borne supérieure n’est pas trop élevée.

Incertitude statistique, données et modèles

Une question fondamentale pour la survie de l’assurance est de savoir où se situe cette borne supérieure : jusqu’où peut-on aller, entre les deux cas extrêmes (population homogène avec une probabilité de 2% pour tous, et une population très discriminée, avec 2% de la population ayant 100% de chances de mourir, et l’autre 0%) ? Et de quoi cette borne dépend-elle ? En particulier, des modèles plus complexes, tels que les réseaux de neurones très profonds permettent-ils vraiment d’améliorer la prévision ? Et l’enrichissement de données, tel qu’on l’observe grâce aux objets connectés et la fusion avec toutes sortes d’informations externes, va-t-il déplacer la borne supérieure vers le haut ?

Si l’apprentissage profond – voir Goodfellow et al. (2018) – permet d’avoir des classifieurs d’images avec un taux d’erreur proche de 0%, il est difficile d’imaginer qu’il sera possible de prévoir, presque un an à l’avance (à la signature du contrat), qui décèdera dans l’année, qui aura la grippe, qui aura un dégât des eaux, etc. Les modèles plus complexes permettent d’améliorer les prévisions, en tenant compte de non-linéarités, d’effets croisés entre les variables tarifaires, mais pas au point de faire disparaître l’aléa. Et tant que l’assurance est envisagée ex-ante (la prime est fixée au début de la période de couverture), il est difficile d’imaginer que rajouter de l’information fera aussi disparaître l’aléa. C’est d’ailleurs le cas pour les tests génétiques qui n’expliquent qu’une (petite) partie du risque de cancer, par exemple. Et rajouter des données revient souvent à rajouter du bruit, ce qui rend le travail d’analyse plus complexe. Cependant, force est de constater que des modèles plus complexes et des données plus riches tendent effectivement à « améliorer » la prévision, en remontant la courbe ROC vers le haut. Mais se pose-t-on les bonnes questions ? Que signifie vraiment une borne très éloignée du cas homogène, sur la diagonale ?

Homogénéité, équité et causalité

Comme nous l’avons vu, la tarification en assurance repose sur une répartition des risques (des contrats) en catégories, au sein desquelles la distribution des pertes peut être estimée, afin de fixer un niveau de primes. La répartition se fait à partir des caractéristiques de l’assuré, et du bien assuré. En retraçant l’histoire de l’assurance, Ewald (1986) montre que les mécanismes de prévoyance se sont mis en place en déplaçant la charge des accidents du travail sur la société : on abandonne l’idée d’une responsabilité individuelle de l’accident en faveur de la solidarité. L’assurance distingue « entre le dommage que subit tel ou tel individu — c’est affaire de chance ou de malchance — et la perte liée au dommage dont l’attribution est, quant à elle, toujours collective et sociale ». Ce principe de solidarité sociale, de mutualisation des risques, fait que le risque (en assurance) est toujours pensé collectivement.

Aujourd’hui, les tarifs sont considérés comme « justes », ou « actuariellement équitables » si chaque prime correspond à la perte attendue (pour ne pas dire « espérée », au sens mathématique) pour chaque assuré. Dans cette perception de l’équité, une hypothèse essentielle est que les classes soient « homogènes ». En effet, dans l’hypothèse inverse, les personnes les moins risquées subventionnent les personnes les plus risquées, ce qui est perçu comme socialement injuste.

On peut décrire cette version de l’équité actuarielle à l’aide de la formule de décomposition de la variance. La variance globale se décompose en effet en deux termes, la variance inter-classes et la variance intra-classes : l’ « équité actuarielle » vise à ce que  les classes de risque soient relativement distinctes les unes des autres, donc une variance inter-classes forte, accompagnée d’une homogénéité des classes, donc une variance intra-classes faible. D’un point de vue statistique, chercher à augmenter l’une est équivalente à faire diminuer l’autre. Cette mécanique n’est pas toujours claire pour des observateurs non avertis ; ainsi dans l’affaire Manhart, un des cas les plus documentés sur la discrimination par le genre en assurance, le juge Stevens affirme : « we focus on fairness to individuals rather than on fairness to classes […] even a true generalization about a class is an insufficient reason for disqualifying an individual to whom the generalization does not apply» (cité dans Anzalone (2016)). Autrement dit, pour la justice, un critère statistique de type « true generalization » ne peut être opposable à un individu.

Une autre critique importante, que l’on retrouve dans la « gender directive », est le lien entre discrimination et causalité. En effet, statistiquement, les actuaires vont chercher des facteurs de classification fortement corrélés avec la sinistralité. Mais il est possible que ces facteurs ne soient qu’un proxy de la vraie variable causale, restée elle inobservée, conduisant à une mauvaise estimation du risque pour certains. Comme le notent Antonio et Charpentier (2017), le genre a ainsi été longtemps utilisé en assurance automobile car très corrélé avec des variables associées au style de conduite et à d’autres variables historiquement non observables (mais qui le sont aujourd’hui grâce aux objets connectés, comme le kilométrage, les heures de conduites, les types de routes utilisés, etc).

Ce lien avec les mécanismes causaux est d’ailleurs relativement profond, et Hacking (1975) y voit une connexion avec la « révolution probabiliste » : on peut assez facilement mettre en évidence des corrélations, mais les causes, si elles existent, nous restent plus opaques. Laplace au début du 19e siècle déclare ainsi que « la probabilité est relative en partie à nos connaissances, en partie à notre ignorance », liant les probabilités à la fois à une vision newtonienne déterministe du monde et à notre incapacité à le connaitre parfaitement. Cette dernière composante fait que l’on ne peut pas annoncer la date exacte du décès d’un individu, mais statistiquement, dans un groupe homogène, on peut prédire le nombre de décès au cours d’une année. Et pour revenir à la relation causale, le tabagisme par exemple ne cause pas forcément une mort prématurée mais fumer sera vu comme dangereux car il augmente la probabilité de décès pendant une période donnée. Ainsi nous montre Hacking (1975), la causalité est pensée aujourd’hui dans un contexte probabiliste, et non plus déterministe.

Abraham, K.S. (1985) Efficiency and fairness in insurance risk classification, Virginia Law Review 71: 403-451.

Anderson, A.W. (1978). A Critique of the Manhart Brief. The Actuary, 12:5.

Antonio, K. & Charpentier, A. (2017). La tarification par genre en assurance, corrélation ou causalité ? Risques, 109.

Anzalone, C.A. (2016). U.S. Supreme Court Cases on Gender and Sexual Equality. Routledge.

Bailey, H., Hutchison, T. & Narber, G. (1975) The regulatory challenge to life insurance classification, Drake Law Review Insurance Law Annual 4: 779-827

Barry L. (2019). Justice ou justesse? L’équité de l’assurance. Working paper #15, chaire PARI.

Charpentier, A. & Denuit, M. (2004). Mathématiques de l’Assurance Non-Vie : Principes Généraux de Théorie du Risque. Economica.

Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

Cummins, J.D., Smith, B.D., Vance, R.N. & VanDerhai, J.L. (1982). Risk Classification in Life Insurance. Kluwer-Nijhoff Publishing.

Ewald F. (1986). L’État providence. Grasset.

Fisher, R. A. (1936).  The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics. 7 (2): 179–188.

Frézal, S. & Barry, L. (2019). Fairness in Uncertainty: Some Limits and Misinterpretations of Actuarial Fairness, Journal of Business Ethics.

Goodfellow, U., Bengio, Y. & Courville, A. (2018) L’apprentissage profond. Massot Edition.

Hacking, I. (1975) The Emergence of Probability. Cambridge University Press

Hoy, M. (1982) Categorizing risks in the insurance industry, The Quarterly Journal of Economics 97: 321-336.

Kleindorfer, P. & Kunreuther, H. (1980) Misinformation and Equilibrium in Insurance Markets, in Economic Analysis of Regulated Markets, Jörg Finsinger Editor, Springer Verlag, 67-90

Kuhn, M. & Johnson, K. (2018). Applied Predictive Modeling. Springer Verlag.

Martin,G.D. (1976) Gender Discrimination in Pension Plans, Journal of Risk and Insurance 43.

McGrayne, S.B. (2012) The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy. Yale University Press.

Ramsey P.F. (1926). Truth and Probability.

Schauer, F. (2006) Profiles, Probabilities, and Stereotypes. Harvard University Press.

Von Mises, R. (1957). Probability, Statistics and Truth. Dover publications.

Wilson, C. (1977). A model of insurance markets with incomplete information. Journal of Economic Theory, 16:2, 167-207.

[1] Au sens statistique du mot, dans le sens introduit par Fisher (1936).

[2] Dans cette approche fréquentiste, et notamment pour Ronald Fisher et Richard von Mises, la probabilité d’un évènement unique (dit « one shot ») n’a pas de sens.

Insurance data science : use and value of unusual data #1

Next week, with , I will be at the Summer School of the Swiss Association of Actuaries, in Lausanne, with Jean-Philippe Boucher (UQAM) and Ewen Gallic (AMSE).

I will give an introductionary talk on Monday morning, and the slides are now available

There will be some hands-on applications, on R. I will share some codes in the slides.

La Ville, Conférence CartoStats

Vendredi, je participerai à la journée CartoStats organisée à Paris Diderot à l’occasion des 60 ans du Comité Français de Cartographie par Françoise Bahoken, Laurent Jégou et Hélène Mathian (entre autres). On m’a proposé une intervention de 20 minutes, et ça sera l’occasion de revenir sur des travaux récents avec Lucas Vernet (doctorant à Sciences Po), Alfred Galichon (professeur à Columbia) et Ewen Gallic (qui vient de défendre son doctorat, la semaine dernière, et qui commencera son post-doc avec moi à la rentrée sur la démographie participative). Et histoire d’avoir une cohérence dans mon exposé, j’ai pris un fil rouge qui sera la ville en économie. Les slides sont en ligne (ils ont \’et\’e mis \`a jour par rapport au billet, avec des r\’ef\’erences compl\’ementaires), mais ils sont purement visuels, et je me suis dit qu’un billet permettrait de les compléter.

Plus que la ville, c’est la notion d’agglomération autour d’un centre qui va m’intéresser – ce que les géographes appellent CBD (Central Business District) ou centre des affaires – et la localisation des vendeurs dans l’espace (souvent dans des “villes”) compte tenu d’une répartition d’acheteurs. On peut remonter aux travaux de von Thünen ou Kohl qui, les premiers, justifient l’existence de villes monocentriques, autour d’un centre.

Mais pour tous les économistes, quand on parle de “ville”, on pense instinctivement à la “ville linéaire” introduite par Hotelling en 1929

avec cette ville incroyablement simpliste, on peut justifier la concentration des firmes en un “centre”. On suppose  ici la population répartie (uniformément) dans cette ville, et les firmes cherchent où se placer, de manière optimale. Et on peut montrer assez facilement qu’avec deux firmes, un équilibre stable est que les villes s’installent toutes les deux au centre.

Mais on l’a dit, cette ville linéaire n’est pas très réaliste, et pas très représentative de la plupart des villes. Cela n’empêche pas certains auteurs de changer les hypothèses pour arriver à d’autres configurations d’équilibres, comme le montre Fujita & Thisse (1997), en reprenant Fujita & Ogawa (1982)

La première tentative pour passer en bi-dimensionnel a été la ville circulaire de Salop (1979), mais on le voit, c’est une fausse extension du cas précédant. On reste assez fondamentalement en une dimension,

Une autre extension pourrait être la ville “étoilée”, avec une population répartie en étoile, le long de routes. Mais là encore, si ce modèle est encore assez populaire, on se demande pourquoi ne pas envisager la création de route périphérique, autour du centre.

Autrement dit, l’économie spatiale et l’économie urbaine ne peuvent se penser sans faire intervenir la dimension des transports (comme le dit le théorème d’impossibilité de Starrett (1978)). L’extension logique est alors de voir la ville sous forme d’un réseau, comme l’avait fait Christaller (1933), afin de commencer une analyse microéconomique de la répartition de l’offre et de la demande.

Mais cette dimension rend le problème plus complexe, car il existe une infinité de réseaux, et tous ont leurs particularité. Et les réseaux ont une particularité : l’intuition que l’on a en terme de distance et de voisins tend à disparaître. Par exemple sur le graphique ci-dessous, 8 points sont répartis sur la région parisienne, et on y accède par métro. La figure du haut correspond au diagramme de Voronoï classique (c’est à dire aux zones d’affluence de magasins si les prix sont identiques) avec une norme euclidienne (à vol d’oiseau). Le graphique du dessous correspond au diagramme de Voronoï si on utilise la distance en métro. On voit que les zones d’affluence ne sont plus les mêmes.

Le problème économique que l’on peut envisager est alors celui d’un matching bipartite, où les “consommateurs” vont se rendre chez un “producteur” (et un seul). On peut imaginer des contraintes de capacité des producteurs, voire des problèmes de congestion sur le réseau, etc.

On a récemment étudié ce genre de modèles avec Alfred et Lucas, tout d’abord avec la localisation des maternités (le bien est ici un accouchement : on sait – par communes – le nombre de femme ayant accouché, et on a des informations sur les maternités) et les stations essence (le bien est ici un plein d’essence : on a la localisation des stations en Bretagne, et la localisation d’utilisateurs d’une api qui souhaitent connaître le prix de l’essence proche de leur localisation). Compte tenu de l’offre et de la demande, il est possible d’en déduire le “prix” de ces biens, et donc leur valeur. Ceci permet, par exemple, de voir où il serait nécessaire d’implanter une nouvelle maternité, par exemple.

Mais l’étude de réseaux ne peut se faire en statistique. L’important est d’ajouter une dimension temporelle. Dans l’exemple de la ville linéaire, on suppose les agents uniformément répartis dans la ville, puis les firmes se localisent, mais comment ne pas imaginer ensuite que cette localisation influencera a son tour la répartition des agents? Beckman (1976) note ainsi que les individus se répartiront au sein de la ville avec “a density distribution which is symmetric and decreasing from the center

C’est la notions d’externalités résultats des interactions entre firmes et/ou ménages.

Cette dimension est assez fondamentale, en témoigne la visualisation de l’épidémie de choléra à New York en 1832.

En effet, la prise en compte de la dimension spatiale permet de voir où est née l’épidémie de choléra, et de montrer l’importance des points d’eau, comme l’avait noté John Snow lors de l’épidémie de Londres en 1854 (raconté dans un article paru la semaine passée).

Les exposés de vendredi après midi étant autour des thèmes “nouvelles méthodes, nouveaux outils” et “nouveaux usages“, ça sera l’occasion de présenter des travaux récents avec Ewen sur la création de cartes évoluant dans le temps.

A titre indicatif, nous allons travailler cet été sur la création de telles cartes pour étudier le déplacement de descendants d’un couple, qui s’est marié en 1800 dans un village du Finistère, du bocage normand, ou du sud de la Bourgogne. A suivre donc….

Automatic Detection of the Language of a Tweet

Two days ago, in my post to extract automatically my own tweets, and to generate some html list, I mentioned that it would be great if there were a function that could be used to distinguish tweets in English, and tweets in French (usually, I tweet in one of those two languages). And one more time, @3wen came to rescue me! In my previous post, I spent some time to extract mentions to other twitter accounts, and urls. This time, I have to remove those two informations, that are not relevant to get the language of a tweet.

Continue reading Automatic Detection of the Language of a Tweet

Names in the U.S., from James Smith to Jose Rodriguez

Two weeks ago, @mona published an interesting post on her blog, about a difficult question, What’s The Most Common Name In America? There were stats about first names, in the U.S., and last names, too. Those informations are – somehow – easy to get. But usually, it is more complicated to get the first and the last name together. For confidentiality issues ! Datasets – the ones I deal with – are supposed to be anonymized, so I never see the first and the last names.  In a previous post, a few years ago, I did mention the so-called Social Security Death Master File. In that file, we have Social Security numbers, with the date of birth, the date of death as well as the first and the last name. So I did use those files to get stats about the first and the last names of American citizens. Of course, it is very restrictive. I have only U.S. citizens that have a Social Security number (which is not compulsary in the U.S. as far as I understood) and who passed away (as mentioned in the name of the dataset: the death master file). Another great thing about that dataset is that I have the date of birth, so I can look at some cohort effect (see opendata.stackexchange for an interesting discussion on that dataset).

Continue reading Names in the U.S., from James Smith to Jose Rodriguez

Shapefiles from Isodensity Curves

Recently, with @3wen, we wanted to play with isodensity curves. The problem is that it is difficult to get – numerically – the equation of the contour (even if we can easily plot it). Consider the following surface (just for fun, in order to illustrate the idea)

> f=function(x,y) x*y+(1-x)*(1-y)
> u=v=seq(0,1,length=21)
> v=seq(0,1,length=11)
> f=outer(u,v,f)
> persp(u,v,f,theta=angle,phi=10,box=TRUE,
+ shade=TRUE,ticktype="detailed",xlab="",
+ ylab="",zlab="",col="yellow")

For instance, assume that we want to locate areas where the density exceed 0.7 (here in the lower left corner, SW, and the upper right corner, NE)

> image(u,v,f)
> contour(u,v,f,add=TRUE,levels=.7)

Continue reading Shapefiles from Isodensity Curves

Extracting datasets from excel files in a zipped folder

The title of the post is a bit long, but that’s the problem I was facing this morning: importing dataset from files, online. I mean, it was not a “problem” (since I can always download, and extract manually the files), more a challenge (I should be able to do it in R, directly). The files are located on ressources-actuarielles.net, in a zip file. Those are mortality tables used in French speaking African countries, and I guess that one problem came from special characters, such as “é” or “è”… When you open the zip file, you see a folder

and in that folder, several files that I would like to import

Continue reading Extracting datasets from excel files in a zipped folder

Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction is now online, on hal.archives-ouvertes.fr/.

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and difficult to implement. We provide a simple technique — based of properties of Gaussian kernels — to efficiently compute weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. We illustrate the use of that technique to visualize hot spots of car accidents and campsite locations, as well as location of bike thefts.

There are new applications, and new graphs, too

Most of the codes can be found on https://github.com/ripleyCorr/Kernel_density_ripley (as well as datasets).

We are winter

Non, pas de billet pour annoncer la reprise de Game of Thrones dans les jours à venir. Mais un billet rapide pour rappeler qu’au Canada, l’hiver, c’est un sujet important, et sérieux.

C’est ce thème qu’avait repris le Comité Olympique Canadien avant et pendant les Jeux.

Et lors de mon excursion au Mexique la semaine passée (où nous avons retrouvé des collègues chiliens d’ailleurs), je me suis rendu compte que les Jeux Olympiques d’hiver n’intéressaient que peu de monde, en dehors des Canadiens, et des nations qui subissent l’hiver, comme les Russes, ou les Norvégiens (entre autres). Un autre point comment de toutes ces nations est d’avoir un intérêt très fort pour un endroit important d’un point de vue géopolitique, le pôle nord. Et on m’a demandé d’intervenir dans une conférence de Thierry Garcin, et Philippe Reka. J’ai toujours un plaisir fou à discuter avec Philippe, donc c’est avec plaisir que je viendrais, mais j’avoue (une nouvelle fois) avoir du mal à trouver ma place (oui, je suis un peu le sparadrap du capitaine Haddock)

Il y a tout juste un an, Philippe publiait dans le Monde Diplo un article inspirant sur la Cartographie Radicale,

Depuis des décennies, la cartographie traditionnelle revendique le statut de science exacte s’appuyant sur des données fiables. Elle se targue de fournir une image neutre et fidèle de la réalité. Mais une telle approche fait l’impasse sur l’utilisation politique et sociale de la carte, et sur son rôle tant de propagande que de contestation. Depuis le début des années 2000, et dans le désordre, émerge une pratique cartographique qui se dit « radicale » (on parle aussi de « cartographie critique » ou de « contre-cartographie »), riche combinaison revendiquée d’art, de sciences, de géographie, de politique et de militantisme social. Les cartographes orthodoxes ont vu naître avec beaucoup de méfiance les premiers projets exploratoires qui bousculaient, voire cassaient, les codes. « On peut “crobarder” des schémas réducteurs donnant l’illusion de la synthèse, et en faire des figures de style prétentieuses, écrivait un cartographe professionnel en écho à l’une des toutes premières expositions de cette mouvance. La carte géographique, ce fantastique outil de savoir et d’appréhension de l’espace et du monde, se voit ravalée au rang de simple outil de prestidigitation. »

C’est peu dire que le symposium « Art and Cartography — Cartography and Art », organisé par l’Académie des beaux-arts et l’Université technique de Vienne (Autriche) en février 2008, fut le théâtre d’intenses discussions, où s’exprimèrent des positions inconciliables. La cartographie engagée a trouvé la voie d’une forme de contestation au service de la justice sociale, dénonçant des pratiques politiques et économiques douteuses. Elle est un exercice libre de déconstruction de l’espace et des phénomènes sociaux, pour lequel les protagonistes se permettent de pervertir les conventions les plus classiques. Les premiers à s’en emparer furent des artistes, mais aussi des architectes, des urbanistes et des militants, rejoints plus tard par des géographes. Les projets, individuels ou collectifs, se sont surtout multipliés en Argentine, aux Etats-Unis, en Suisse, en Allemagne et en France, couvrant un champ thématique assez large (finance, surveillance et sécurité, marketing, environnement, organisation des territoires, pour ne donner que quelques exemples), et profitant largement de la puissance des outils de la cartographie participative et des réseaux sociaux.

Déchiffrer pour dénoncer, mettre en images des processus peu visibles qui concourent à confisquer l’espace public (voire des biens publics), à compromettre les libertés individuelles, à détourner des lois, tels sont les principaux objectifs des promoteurs de cette cartographie. Dans cet esprit, informer n’est qu’un début. L’étape suivante de la démarche reste l’action pour le changement sur le terrain. Largement informelle, cette initiative se reconnaît quelques affinités avec les mouvements de résistance tels que celui des « indignés » ou Occupy. Ainsi, le collectif d’artistes argentins Grupo de Arte Callejero ne s’est pas contenté de dénoncer les bourreaux de l’ancienne dictature en cartographiant les lieux de leur paisible retraite : il a aussi organisé des itinéraires « touristiques » qui passaient par leurs domiciles et mené campagne sur des panneaux publicitaires… Le projet « Duty Free Shop », en ce qu’il se propose de déconstruire, rendre lisible et donc fragiliser un processus d’aliénation des citoyens dans un espace public, s’inscrit de plain-pied dans ce mouvement pétillant.

En tant qu’amateur (plus qu’expert) de visualisation de données, je rejoins complètement Philippe dans sa lecture critique (je ne vois les cartes que comme un support graphique parmi beaucoup d’autres). Philippe a ainsi proposé une très jolie visualisation de l’arctique.

Et comme souvent, Philippe a raison. Pour comprendre l’arctique, il faut avoir une visualisation claire. Car comme beaucoup de monde, j’ai été éduqué avec des cartes qui ressemblaient à ça

Autrement dit, pour moi, l’arctique et le pôle nord ont longtemps été des concepts flous, un peu comme l’horizon… c’est au nord, tout là haut. Par contre, j’ai bien conscience que l’antarctique existe, par exemple, mais les cartes usuelles ne nous permettent de prendre conscience de sa superficie (un peu moins de deux fois la taille de l’Australie). D’ailleurs, si on va sur https://www.google.ca/maps/place/North+Pole/ on peut voir (en butant contre la partie supérieure de la carte) qu’on ne voit rien, si on cherche à visualiser ce pôle nord (même s’il est sous la glace)

Comme le notait ma fille, le pôle nord est une enjeu stratégique important car c’est là que se trouve le village du père Noël (c’est manifestement la croyance répandue en Amérique). Et si le père Noël n’existe pas (c’est ce que pense ma fille), peut être que le pôle nord n’existe pas non plus ? Le premier défi est donc d’avoir une visualisation aussi pertinente que possible du pôle. Je vais donc revenir un peu sur les problématique de projection. Car les Soviétiques ne se trompaient pas lorsqu’ils représentèrent l’arctique avec un point de vue qui donnait l’impression que l’arctique était russe,

Pour appréhender les projecteurs, deux ouvrages fabuleux expliquent cela très clairement. Flattening the Earth – Two Thousand Years of Map Projections de John Parr Snyder, mais aussi (et surtout) How to Lie with Maps de Mark Monmonier. Mark Monmonier a également publié un livre passionnant sur le sujet (en particulier le chapitre 6), intitulé Air Apparent – How Meteorologists Learned to Map, Predict, and Dramatize Weather

Pour visualiser le pôle nord, on peut utiliser quelques fonctions simples (basées sur un travail en cours avec @3wen), pour visualiser différentes projections. On va ici mettre l’accent sur trois pays importantes quand on parle du pôle nord, le Canada, la Russie, et le Groenland (qui est techniquement un pays consultatif), et sur trois deux régions, une située à moins de 1500km du pôle nord, et l’autre à 3000km du pôle (soit un peu en dessous du cercle arctique)

Si on veut changer de point de vue, on peut visualiser l’animation suivante

D’ailleurs, si on compare la proportion du disque occupé par chacun des pays (j’ai pris ici chacune des huit pays de l’arctique), on obtient, avec en abscisse la distance au pôle, et en ordonnée, la proportion du disque occupée par chacun des pays, avec le Canada (rouge), la Russie (bleu) et le Groenland (vert)

Techniquement, le pôle pose de réels soucis car la plupart des données sont représentées par leur latitude, et leur longitude. Les pôles sont les seuls points – dans cette représentation – à avoir une infinité de coordonnées !

File:Latitude and Longitude of the Earth.svg

La latitude est de 90°, mais toutes les longitudes sont valides. On a ainsi des soucis y compris sur les cartes construites à partir de données satellites, qui pourtant pourraient balayer uniformément la sphère.

Le pôle nord peut ainsi compter beaucoup de discontinuités, uniquement à cause de la représentation des données par la latitude et la longitude : en quelques kilomètres, on passe d’une zone rouge à une zone blanche, ce qui pose des soucis de lissage de la carte (les données n’étant manifestement pas considérées ici comme circulaire pour le lissage).

Maintenant, plus qu’un point sur la carte, on peut voir le pôle nord comme une région, et non pas comme un point sur la sphère, défini très arbitrairement. Bon, j’exagère un peu ici, car le pôle nord est défini par l’axe de rotation de la terre,

qui va induire une espère de symétrie, par exemple sur le climat. Mais le second problème, selon moi, qui va alors apparaître est que le pôle nord vu comme un territoire (un territoire de glace) est alors un territoire qui va changer avec les saisons.

Par exemple, en septembre dernier, la glace occupait la région suivante, ce qui montre que peu de pays peuvent effectivement rejoindre le pôle nord ‘à pieds’ quel que soit le moment de l’année.

On continue de travailler avec @3wen sur le sujet, mais le sujet et vaste, ambitieux… et complexe ! à suivre donc…

Visualizing densities of spatial processes

We recently uploaded on http://hal.archives-ouvertes.fr/hal-00725090 a revised version of our work, with Ewen Gallic (a.k.a. @3wen) on Visualizing spatial processes using Ripley’s correction: an application to bodily-injury car accident location

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and di#cult to implement. We provide a simple technique – based of properties of Gaussian kernels – to compute e#efficiently weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. An illustration on location of bodily-injury car accident (and hot spots) in the western part of France is discussed, where a lot of accident occur close to large cities, next to the sea.

Sketches of the R code can be found in the paper, to produce maps, an to describe the impact of our boundary correction. For instance, in Finistère, the distribution of car accident is the following (with a standard kernel on the left, and with correction on the right), with 186 claims (involving bodily injury)

and in Morbihan with 180 claims, observed in a specific year (2008 as far as I remember),

The code is the same as the one mentioned last year, except perhaps plotting functions. First, one needs to defi
ne a color scale and associated breaks

breaks <- seq( min( result $ZNA , na.rm = TRUE ) * 0.95 , max ( result$ZNA , na.rm = TRUE ) * 1.05 , length = 21)
col <- rev( heat . colors (20) )

to
finally plot the estimation

image . plot ( result $X, result $Y, result $ZNA , xlim = range (pol[,
1]) , ylim = range (pol[, 2]) , breaks = breaks , col = col ,
xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "n",
zlim = range ( breaks ), horizontal = TRUE )

It is possible to add a contour, the observations, and the border of the polygon

contour ( result $X, result $Y, result $ZNA , add = TRUE , col = "grey ")
points (X[, 1], X[, 2], pch = 19, cex = 0.2 , col = " dodger blue")
polygon (pol , lwd = 2)

Now, if one wants to improve the aesthetics of the map, by adding a Google Maps base map, the
first thing to do – after loading ggmap package – is to get the base map

theMap <- get_map( location =c( left =min (pol [ ,1]) , bottom =min (pol[ ,2]) , right =max (pol [ ,1]) , 
top =max (pol [ ,2])), source =" google ", messaging =F, color ="bw")

Of course, data need to be put in the right format

getMelt <- function ( smoothed ){
res <- melt ( smoothed $ZNA)
res [ ,1] <- smoothed $X[res [ ,1]]
res [ ,2] <- smoothed $Y[res [ ,2]]
names (res) <- list ("X","Y","ZNA")
return (res )
}
smCont <- getMelt ( result )

Breaks and labels should be prepared

theLabels <- round (breaks ,2)
indLabels <- floor (seq (1, length ( theLabels ),length .out =5)) 
indLabels [ length ( indLabels )] <- length ( theLabels ) 
theLabels <- as. character ( theLabels [ indLabels ])
theLabels [ theLabels =="0"] <- " 0.00 "

Now, the map can be built

P <- ggmap ( theMap )
P <- P + geom _ point (aes(x=X, y=Y, col=ZNA), alpha =.3 , data =
smCont [!is.na( smCont $ZNA ) ,], na.rm= TRUE )

It is possible to add a contour

P <- P + geom _ contour ( data = smCont [!is.na( smCont $ZNA) ,] ,aes(x=
X, y=Y, z=ZNA ), alpha =0.5 , colour =" white ")

and colors need to be updated

P <- P + scale _ colour _ gradient ( name ="", low=" yellow ", high ="
red", breaks = breaks [ indLabels ], limits = range ( breaks ),
labels = theLabels )

To remove the axis legends and labels, the theme should be updated

P <- P + theme ( panel . grid . minor = element _ line ( colour =NA), panel
. grid . minor = element _ line ( colour =NA), panel . background =
element _ rect ( fill =NA , colour =NA), axis . text .x= element _ blank() ,
axis . text .y= element _ blank () , axis . ticks .x= element _ blank() ,
axis . ticks .y= element _ blank () , axis . title = element _ blank() , rect = element _ blank ())

The
final step, in order to draw the border of the polygon

polDF <- data . frame (pol)
colnames ( polDF ) <- list ("lon","lat")
(P <- P + geom _ polygon ( data =polDF , mapping =( aes(x=lon , y=lat)), colour =" black ", fill =NA))

Then, we’ve applied that methodology to estimate the road network density in those two regions, in order to understand if high intensity means that it is a dangerous area, or if it simply because there is a lot of traffic (more traffic, more accident),

We have been using the dataset obtained from the Geofabrik website which provides
Open-StreetMap data. Each observation is a section of a road, and contains a few points identifi
ed by their geographical coordinates that allow to draw lines. We have use those points to estimate a proxy of road intensity, with weight going from 10 (highways) to 1 (service roads).

splitroad <- function ( listroad , h = 0.0025) {
pts = NULL
weights <- types . weights [ match ( unique ( listroad $ type ), types .
weights $ type ), " weight "]
for (i in 1:( length ( listroad ) - 1)) {
d = diag (as. matrix ( dist ( listroad [[i]]))[, 2: nrow ( listroad
[[i ]]) ]))
}}
return (pts )
}

See Ewen’s blog for more details on the code, http://editerna.free.fr/blog/…. Note that Ewen did publish a poster of the paper (in French), for the http://r2013-lyon.sciencesconf.org/ conference, that will be organized in Lyon on June 27th-28th, see

All comments are welcome…

Border bias and weighted kernels

With Ewen (aka @3wen), not only we have been playing on Twitter this month, we have also been working on kernel estimation for densities of spatial processes. Actually, it is only a part of what he was working on, but that part on kernel estimation has been the opportunity to write a short paper, that can now be downloaded on hal.

The problem with kernels is that kernel density estimators suffer a strong bias on borders. And with geographic data, it is not uncommon to have observations very close to the border (frontier, or ocean). With standard kernels, some weight is allocated outside the area: the density does not sum to one. And we should not look for a global correction, but for a local one. So we should use weighted kernel estimators (see on hal for more details). The problem that weights can be difficult to derive, when the shape of the support is a strange polygon. The idea is to use a property of product Gaussian kernels (with identical bandwidth) i.e. with the interpretation of having noisy observation, we can use the property of circular isodensity curve. And this can be related to Ripley (1977) circumferential correction. And the good point is that, with R, it is extremely simple to get the area of the intersection of two polygons. But we need to upload some R packages first,

require(maps)
require(sp)
require(snow)
require(ellipse)
require(ks)
require(gpclib)
require(rgeos)
require(fields)

To be more clear, let us illustrate that technique on a nice example. For instance, consider some bodiliy injury car accidents in France, in 2008 (that I cannot upload but I can upload a random sample),

base_cara=read.table(
"http://freakonometrics.blog.free.fr/public/base_fin_morb.txt",
sep=";",header=TRUE)

The border of the support of our distribution of car accidents will be the contour of the Finistère departement, that can be found in standard packages

geoloc=read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
header=TRUE,sep=",",comment.char="",check.names=FALSE,
colClasses=c(rep("character",5),rep("numeric",38)))
geoloc=geoloc[,c("dep","com","com_nom",
"long","lat","pop_2008")]
geoloc$id=paste(sprintf("%02s",geoloc$dep),
sprintf("%03s",geoloc$com),sep="")
geoloc=geoloc[,c("com_nom","long","lat","pop_2008")]
head(geoloc)
france=map('france',namesonly=TRUE,
plot=FALSE)
francemap=map('france', fill=TRUE, col="transparent",
plot=FALSE)
detpartement_bzh=france[which(france%in%
c("Finistere","Morbihan","Ille-et-Vilaine",
"Cotes-Darmor"))]
bretagne=map('france',regions=detpartement_bzh,
fill=TRUE, col="transparent", plot=FALSE,exact=TRUE)
finistere=cbind(bretagne$x[321:678],bretagne$y[321:678])
FINISTERE=map('france',regions="Finistere", fill=TRUE,
col="transparent", plot=FALSE,exact=TRUE)
monFINISTERE=cbind(FINISTERE$x[c(8:414)],FINISTERE$y[c(8:414)])

Now we need simple functions,

cercle=function(n=200,centre=c(0,0),rayon)
{theta=seq(0,2*pi,length=100)
m=cbind(cos(theta),sin(theta))*rayon
m[,1]=m[,1]+centre[1]
m[,2]=m[,2]+centre[2]
names(m)=c("x","y")
return(m)}
poids=function(x,h,POL)
{leCercle=cercle(centre=x,rayon=5/pi*h)
POLcercle=as(leCercle, "gpc.poly")
return(area.poly(intersect(POL,POLcercle))/
area.poly(POLcercle))}
lissage = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,
sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)

before defining our weights.

poidsU=function(i,U,h,POL)
{x=U[i,]
poids(x,h,POL)}
OMEGA=parLapply(cl,1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)
stopCluster(cl)
}else
{OMEGA=lapply(1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)}

Note that it is possible to parallelize if there are a lot of observations,

if(n>=500)
{cl <- makeCluster(4,type="SOCK")
worker.init <- function(packages)
{for(p in packages){library(p, character.only=T)}
NULL}
clusterCall(cl, worker.init, c("gpclib","sp"))
clusterExport(cl,c("cercle","poids"))

Then, we can use standard bivariate kernel smoothing functions, but with the weights we just calculated, using a simple technique that can be related to one suggested in Ripley (1977),

fhat=kde(U,H,w=1/OMEGA,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
fhat$estimate=fhat$estimate*sum(1/OMEGA)/n
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}
lissage_without_c = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)}
fhat=kde(U,H,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}

So, now we can play with those functions,

base_cara_FINISTERE=base_cara[which(point.in.polygon(
base_cara$long,base_cara$lat,monFINISTERE[,1],
monFINISTERE[,2])==1),]
coord=cbind(as.numeric(base_cara_FINISTERE$long),
as.numeric(base_cara_FINISTERE$lat))
nrow(coord)
map(francemap)
lissage_FIN_withoutc=lissage_without_c(coord,
monFINISTERE,optimal=TRUE)
lissage_FIN=lissage(coord,monFINISTERE,
optimal=TRUE)
lesBreaks_sans_pop=range(c(
range(lissage_FIN_withoutc$Z),
range(lissage_FIN$Z)))
lesBreaks_sans_pop=seq(min(lesBreaks_sans_pop)*.95,
max(lesBreaks_sans_pop)*1.05,length=21)

plot_article=function(lissage,breaks,
polygone,coord){
par(mar=c(3,1,3,1))
image.plot(lissage$X,lissage$Y,(lissage$ZNA),
xlim=range(polygone[,1]),ylim=range(polygone[,2]),
breaks=breaks, col=rev(heat.colors(20)),xlab="",
ylab="",xaxt="n",yaxt="n",bty="n",zlim=range(breaks),
horizontal=TRUE)
contour(lissage$X,lissage$Y,lissage$ZNA,add=TRUE,
col="grey")
points(coord[,1],coord[,2],pch=19,cex=.1,
col="dodger blue")
polygon(polygone,lwd=2,)}

plot_article(lissage_FIN_withoutc,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

plot_article(lissage_FIN,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

If we look at the graphs, we have the following densities of car accident, with a standard kernel on the left, and our proposal on the right (with local weight adjustment when the estimation is done next to the border of the region of interest),

Similarly, in Morbihan,

With those modified kernels, hot spots appear much more clearly. For more details, the paper is online on hal.