Tag Archives: cartes

L’assurance, l’immobilier et la fabrique de l’ignorance

Cette fin de semaine, le New York Times rapportait que Zillow, premier site immobilier américain, a retiré discrètement de ses annonces les scores de risques climatiques (inondation, incendie, chaleur) produits par First Street Foundation. Le retrait suit les pressions d’agents immobiliers et le lobbying de plusieurs Multiple Listing Services (MLS), notamment celui de Californie, l’un des plus puissants du pays. Comme l’écrit la journaliste Claire Brown, Zillow a supprimé ces scores suite à une plainte,

“Zillow stopped displaying the scores after complaints from the California Regional Multiple Listing Service”

Mais que se passe-t-il lorsqu’une information sur les risques entre en collision avec des intérêts économiques solidement installés ? Mary Douglas le soulignait, en 1992, dans Risk and Blame

institutions themselves are the monitors which determine what is going to count as information

Le risque, autrement dit, ne se contente pas d’exister ; il est façonné, filtré, parfois instrumentalisé comme une technologie morale, c’est-à-dire un instrument utilisé pour orienter les comportements en fonction des intérêts d’un groupe. Le cas Zillow en est un exemple presque didactique : la connaissance produite par les modèles climatiques devient un coût pour l’industrie immobilière, donc elle doit être rendue silencieuse. Dans Risk and Culture, Mary Douglas écrivait :

“Blaming the system is a way of protecting the community from morally unacceptable risks. […] What is dangerous is not the event itself, but the threat it poses to an established social order.”

Les risques climatiques, dans cette perspective, menacent une “social order” : la fluidité des marchés immobiliers. Ils deviennent donc socialement inacceptables… non pas parce qu’ils sont faux, mais parce qu’ils sont vrais et dérangeants.

xL’objection immobilière : quand l’économie préfère l’ignorance

Le PDG du MLS californien déclare que

“When we saw entire neighborhoods with a 50 percent probability of the home flooding this year and a 99 percent probability of the home flooding in the next five years, especially in areas that haven’t flooded in the last 40 to 50 years, we grew very suspicious,”

L’argument est ancien et classique : un modèle est suspect lorsqu’il contredit l’expérience vécue. Or de nombreux travaux de climatologie et d’hydrologie rappellent que l’expérience passée n’est plus un guide fiable dans une distribution de risques non stationnaire. C’est précisément la thèse de tous les ouvrages sur les régimes extrêmes, la loi des queues épaisses ou le “stationarity is  dead” de Milly et al. (2008). Les hydrologues, les assureurs et les climatologues le répètent depuis 20 ans : le passé n’est plus un indicateur du futur. Dans The Physics of Climate Change, Lawrence Krauss écrit :

“The past climate is no longer a faithful record of the future—it is a fossil.”

Le commentaire du MLS illustre l’une des limites épistémologiques mises en avant par les philosophes des risques : la circularité de la validation. Si l’on exige qu’un risque futur soit déjà observé dans le passé, alors aucun risque émergent n’est légitime. Cette logique fut historiquement mobilisée pour nier les risques industriels dans les années 1950–1980 (asbestos, tabac, pollution). Sheila Jasanoff parle à ce propos de technologies of humility : la nécessité de reconnaître que l’absence d’événement n’est pas un argument.

La granularité excessive du risque

Le New York Times met en évidence, très justement, une tension fondamentale :  les modèles de risques climatiques ne sont pas conçus pour être lus au niveau de la maison individuelle. Madison Condon, chercheuse en droit et économie du climat, dit dans l’article

“With forward-looking models, we can’t check their outputs against a record, they are necessarily forecasting events that have not yet happened. The First Street models might provide a good enough answer for certain questions or certain hazards, But what level of accuracy is good enough changes substantially if the question is about one specific property you are about to spend your life savings on.”

Cette remarque rejoint deux axes centraux de la littérature scientifiques.  Tout d’abord, Les modèles de catastrophe sont fondamentalement des modèles de population, pas de micro-prédiction. Ils sont efficaces en agrégé, mais imprécis en micro-lieu. Mais il y a aussi l’illusion de précision (Lavine, Gigerenzer) : donner un chiffre très précis à un niveau spatial très fin crée une illusion d’exactitude, et donc une vulnérabilité cognitive.

Dans beaucoup d’ouvrages scientifiques sur le climat (notamment les rapports sur les risques de montée des eaux et les travaux de climat extrême), on retrouve cette idée : le risque est distribué spatialement, mais la vulnérabilité est locale et socio-économique. Le modèle ne peut donc être lu sans tenir compte du bâti, des protections locales, des pratiques d’assurance, et surtout des inégalités.

Les “cartographies interdites” 

De nombreux chercheurs ont montré que la publication d’une information sur les risques suscite presque toujours une controverse sociale. Comme le soulignent Reghezza-Zitt et Rufat dans Resilience Imperative,

“any representation of risk involves choices: what is shown, what is omitted, what is made salient”

Une carte n’est donc jamais un simple reflet du réel : elle sélectionne, ordonne et rend visibles certains dangers, tout en en invisibilisant d’autres, et contribue de ce fait à reconfigurer les rapports de force autour du territoire qu’elle représente. Pour Bankoff, Frerks et Hilhorst dans Mapping Vulnerability,

“vulnerability… defies the possibility of reading it from a general chart”

autrement dit, aucune carte ne peut prétendre offrir une représentation neutre ou complète du danger. Toute cartographie du risque simplifie, met en avant certains éléments, en omet d’autres, et participe ainsi à structurer les rapports sociaux autour de ce risque. Les débats sur Zillow s’inscrivent exactement dans cette dynamique. Rendre visible un risque climatique n’est jamais neutre, cela crée des gagnants et des perdants, et l’immobilier fait partie des secteurs où cette visibilité est la plus politiquement sensible. Trois précédents historiques viennent immédiatement à l’esprit.

L’histoire de la cartographie des inondations en France au XIXᵉ siècle montre déjà à quel point rendre visible un danger pouvait déclencher une crise politique locale. Après la crue de 1856, l’État et les ingénieurs des Ponts et Chaussées engagent un vaste travail de délimitation des zones inondables. Très vite, des communes protestent : être placées en “zone inondable”, c’est voir menacés la valeur foncière, l’attractivité de la commune, parfois même sa capacité à emprunter et à investir. Les travaux d’histoire des catastrophes et des politiques hydraulique que tu m’as fournis insistent sur ce point : la carte ne se contente pas d’enregistrer un risque naturel, elle reconfigure les rapports de force locaux en désignant des territoires comme “à problème”. Ce que l’on observe déjà à cette époque, c’est une grammaire politique très stable : la représentation savante du risque est perçue comme un acte performatif, qui fait exister publiquement une menace jusque-là socialement déniée.

Un siècle plus tard, les “redlining maps” de la Home Owners’ Loan Corporation aux États-Unis offrent un autre précédent, encore plus brutal, de la puissance performative des cartes. En classant les quartiers de A à D, avec les quartiers majoritairement noirs et pauvres systématiquement relégués en catégorie “hazardous”, l’administration ne se contente pas de décrire un risque de défaut de crédit : elle organise une véritable ségrégation financière et immobilière. La littérature historique et économique sur le redlining montre bien que ces cartes ont fixé durablement les trajectoires urbaines : un quartier coloré en rouge était privé d’accès au crédit, ce qui entraînait dégradation du bâti, désinvestissement public, puis renforcement des stigmates sociaux. Là encore, le parallèle avec les scores de risques climatiques est évident : dès qu’une représentation du risque est attachée à un lieu, elle pèse sur les prix, les flux d’investissements et, in fine, sur les populations qui y vivent.

Plus récemment, le cas des cartes de FEMA illustre une troisième configuration, où l’enjeu n’est plus tant la production de la carte que la gestion de son obsolescence. Plusieurs travaux que tu m’as transmis soulignent combien de municipalités américaines résistent à la mise à jour des zones inondables : une nouvelle carte signifie, très concrètement, hausse des primes d’assurance, contraintes de construction supplémentaires et baisse immédiate de la valeur immobilière. Il en résulte des cartes souvent en retard sur la réalité climatique, non par manque de capacité technique, mais parce que l’ajustement à la nouvelle distribution des risques est politiquement et économiquement explosif.

Vu depuis ces trois précédents historiques, la décision de Zillow de retirer les scores de risques des annonces s’inscrit dans une continuité troublante : plutôt que d’assumer les effets redistributifs d’une information climatique plus précise, on préfère en limiter la visibilité, au prix d’un retour à une asymétrie d’information défavorable aux acheteurs les plus exposés.

Les modèles First Street : science robuste, usages fragiles

Revenons un instant sur les données et les modèles mentionnés dans l’article du New York Times. Il existe un consensus large, dans la littérature scientifique et dans les évaluations institutionnelles américaines, sur le fait que les cartes FEMA (Flood Insurance Rate Maps — FIRMs) sont souvent obsolètes, parfois très anciennes, et méthodologiquement limitées. Le Department of Homeland Security lui-même, dont dépend FEMA, notait dès 2017 que plus de la moitié des cartes d’inondation américaines avaient plus de 15 ans, et qu’un nombre significatif n’avaient pas été mises à jour depuis 30 à 40 ans (DHS, National Flood Mapping Program Review, 2017). Le Government Accountability Office (GAO) confirme ce diagnostic : en 2022, il rapporte que

“FEMA’s flood maps remain inaccurate or outdated in many areas”

en raison de capacités limitées et de résistances locales à la révision (GAO-22-104, 2022). La littérature académique dit la même chose : dans Understanding Disaster Insurance, Kousky (2020) rappelle que

“FEMA’s flood maps remain incomplete and often outdated; they were never designed to provide a full accounting of flood risk, particularly under a changing climate.”

ou encore Smiley (2020)

“outdated and inaccurate, FEMA flood maps fail to fully capture risk”

Cela signifie que les cartes FEMA utilisent souvent des données hydrologiques, topographiques et pluviométriques qui ne reflètent plus les conditions environnementales actuelles, ni a fortiori celles induites par le changement climatique.

À l’inverse, la méthodologie de First Street Foundation s’appuie sur des modèles hydrologiques et hydrauliques continus, fondés sur des DEM haute résolution (USGS + LiDAR), l’intégration de l’imperméabilisation urbaine, des intensités de pluie mises à jour (NOAA Atlas 14), des données de débordement fluvial, des modèles côtiers (surcote, marées, élévation du niveau de la mer), et surtout des projections climatiques CMIP6. Matthew Eby, de First Street Foundation, est interrogé dans l’article du New York Times,  et comme il le dit,

“Our models are built on transparent, peer-reviewed science and the full methodologies are publicly available for anyone to review on our website.”

Effectivement First Street Foundation met réellement en ligne une méthodologie complète, notamment via le document FSF Flood Model – Technical Methodology, Version 3.0 (juillet 2023), accessible publiquement. On y trouve la description détaillée des données d’entrée : topographie haute résolution (USGS DEM, LiDAR quand disponible), hydrologie et hydraulique (débits USGS, marées et surcotes NOAA, intensités de pluie NOAA Atlas 14), occupation du sol (imperméabilisation, couvert végétal), ainsi que les projections climatiques issues de modèles CMIP6 et de données downscalées. Toutes ces sources sont listées dans les pages “Methodology” et “Data Sources” du site de FSF. Les données historiques d’inondation servent à la validation mais non au calcul direct du “Flood Factor” d’une propriété (FSF précise explicitement que les “historic flood events” ne déterminent pas la note d’un bien individuel).

Cette transparence documentaire est réelle, mais la reproductibilité intégrale du modèle reste limitée. D’une part, le modèle combine plusieurs couches (pluvial, fluvial, côtier, ruissellement urbain) à une résolution très fine (3 m), ce qui exige un traitement SIG massif et un accès à des données parfois difficiles à reconstituer dans leur état exact. D’autre part, les projections climatiques reposent sur des hypothèses (scénarios d’émissions, intensité des tempêtes, montée des eaux) que FSF reconnaît comme sources d’incertitude. Enfin, le modèle ne capture pas toutes les données locales de vulnérabilité : rehaussement individuel des maisons, drainage municipal récent, protections privées, etc. Certains analyses indépendantes (on peut penser à Review of a New Online Modeling Tool) recommandent donc de considérer les résultats comme indicateurs de risque probabiliste, utiles pour corriger l’asymétrie d’information, mais non comme diagnostics structurels à l’échelle de la maison.

Leur méthodologie est en réalité nettement plus avancée que les cartes FEMA, souvent vieilles de 20 à 40 ans.  Dans Climate Shock, Gernot Wagner et Martin Weitzman écrivent

“Inaction is often justified by questioning the models—especially when the newer models are more precise than the legacy systems they threaten.”

Les critiques adressées aux modèles First Street ne portent donc pas sur leur validité scientifique, mais sur leur impact économique.

Quand l’information change les comportements

Le New York Times mentionne un résultat majeur issu d’une expérience menée par Redfin et analysée dans Disaster on the horizon : rendre visible l’information sur le risque modifie de manière mesurable et significative les comportements d’achat. L’expérience, conduite auprès de 18 millions d’utilisateurs, consistait à afficher aléatoirement un score de risque d’inondation provenant de First Street Foundation sur les pages d’annonces. Les chercheurs montrent que, chez les utilisateurs exposés à l’information, les recherches se déplacent vers des zones à plus faible risque, et que les transactions portant sur des biens identifiés comme “à haut risque” enregistrent une baisse moyenne de 1 % du prix de vente (Bernstein, Gustafson & Lewis, 2024). L’étude conclut que

“the disclosure of flood risk information leads to a statistically significant reduction in demand for high-risk properties”

confirmant une intuition déjà présente dans les travaux de psychologie du risque : lorsque l’information devient saillante, elle modifie l’arbitrage des acheteurs. Kunreuther & Meyer le disaient dans The Ostrich Paradox, (2017)

“Markets underreact to low-probability high-consequence risks when information is not salient.”

Le New York Times résume correctement cet effet en notant

“The Redfin experiment, which lasted for three months, affected sales of 8,150 properties that were listed as being at high risk of flooding and reduced their total sales prices by about 1 percent of their value.”

Ce résultat n’est pas anodin : il établit empiriquement que rendre visible le risque climatique crée un effet de marché immédiat, ce qui explique pourquoi la visibilité même de cette information suscite des résistances politiques et économiques, exactement comme dans les précédents historiques que nous avons rappelés. Le problème est donc clair : soit on donne l’information, et le marché ajuste ; soit on la retire, et le marché fonctionne “normalement”… mais au prix de désinformer les acheteurs.

Qui a le droit de savoir ?

Ce contraste entre le faible devoir de transparence des vendeurs immobiliers et les attentes bien plus rigides souvent imposées aux assurés mérite qu’on s’y arrête. Comme je le rappelle dans un billet récent, L’impossible droit à l’erreur, l’impossible droit à l’oubli ?, les systèmes américains d’information publique sont d’une sévérité remarquable vis-à-vis des individus : dans de nombreux États, un assuré peut être lourdement sanctionné s’il omet de déclarer un élément de son passé, notamment des condamnations, des accusations encore pendantes, ou même des arrestations récentes, pourtant non suivies de jugement. Les bases de données publiques — Missouri Sunshine Law (§ 610.100), New York Criminal History Record Search, et les répertoires agrégés par des courtiers privés — transforment chaque trace pénale en stigmate durable, librement accessible à tout employeur, assureur ou bailleur. Comme je l’écrivais

“une simple arrestation (même classée sans suite) peut se muer en exclusion durable”

Autrement dit, l’individu doit “tout dire”, faute de quoi il risque d’être exclu d’un prêt, d’une assurance, ou d’un emploi.

L’asymétrie est frappante puisque l’acheteur immobilier n’a aucun moyen légal d’exiger du vendeur qu’il déclare un antécédent d’inondation ou une vulnérabilité structurelle aux feux. Le New York Times le dit explicitement

“In many states, sellers are not required to disclose information about whether a house has flooded recently or if it is vulnerable to wildfires”

Cette dissymétrie dans le rapport de force est d’autant plus choquante que, dans un cas, l’enjeu est un risque structurel, durable, affectant potentiellement la sécurité future d’un foyer ; dans l’autre, il s’agit souvent de micro-événements individuels dont la pertinence assurantielle ou actuarielle est discutable. La transparence exigée du particulier contraste violemment avec l’opacité tolérée du vendeur immobilier : le premier ne peut invoquer ni droit à l’erreur ni droit à l’oubli, tandis que le second est libre de taire des informations essentielles sur l’exposition réelle du bien au changement climatique.

On touche ici à une injustice structurelle que révèle l’affaire Zillow. L’information climatique, parce qu’elle affecte immédiatement la valeur immobilière, est politiquement neutralisée ; l’information pénale ou personnelle, parce qu’elle permet d’ajuster finement le risque individuel, est, elle, amplifiée, agrégée, rendue indélébile. Dans un environnement où, comme le montre l’expérience Redfin, la seule visibilité d’un risque suffit à modifier le marché, l’effacement volontaire de ces informations dans l’immobilier apparaît comme un privilège que n’ont pas les assurés ordinaires. Ils doivent supporter la mémoire totale du système, alors même que les biens immobiliers, eux, bénéficient d’une forme de droit à l’oubli institutionnel. Cette tension éclaire d’un jour cru le retrait des scores climatiques par Zillow : ce n’est pas seulement un débat technique sur la précision des modèles, mais l’expression d’un régime de responsabilité profondément asymétrique. Dans Climate Change and Cities, Jesse Keenan écrit

“Low-income households are more likely to be harmed by insufficient access to information about climate-related hazards and adaptation options.”

et il est difficile de ne pas être d’accord…

La politique de la visibilité

Cela dit, il convient d’être précis : Zillow n’a pas entièrement supprimé l’accès aux scores, il renvoie désormais vers le site de First Street via un lien discret. Comme le dit le New York Times,

“Real estate listings on Zillow now display hyperlinks to First Street’s website, and users can click through to view climate risk scores for a specific property”

C’est un geste politique très classique : on rend l’information disponible, mais invisible. C’est l’équivalent numérique de mettre un rapport en ligne dans un onglet “archives”. La sociologie des sciences parle de visibility management : gérer la visibilité d’un risque, non pour le corriger, mais pour contrôler son effet social.

La dépolitisation du risque climatique

La décision de Zillow révèle une tension fondamentale dans les sociétés du XXIᵉ siècle. Les modèles climatiques deviennent assez précis pour produire des risques localisés. Mais l’économie et la politique ne sont pas prêtes à intégrer les conséquences de ces risques. Donc l’information est atténuée, déplacée, rendue non saillante. Toute la littérature scientifique récent converge vers une même idée (qu’il s’agisse d’épistémologie du risque, d’anthropologie des catastrophes, d’économie du climat ou d’histoire des assurances) : le risque n’est jamais seulement un fait naturel, c’est un fait politique.

Le cas Zillow n’est pas un détail du marché immobilier. C’est un laboratoire de la gouvernance du risque climatique, avec cette question essentielle : pouvons-nous réellement entrer dans l’ère de l’adaptation climatique si l’information sur les risques devient politiquement indésirable ? Comme l’écrivait Ulrich Beck dans Risikogesellschaft

“Dabei wird exemplarisch deutlich, worum es im öffentlichen Definitionsstreit um Risiken eigentlich geht: nicht nur um gesundheitliche Folgeprobleme … sondern um die sozialen, wirtschaftlichen und politischen Nebenfolgen dieser Nebenfolgen.”

autrement dit,

“Le cœur du conflit public autour des risques n’est pas seulement leurs effets sanitaires, mais les conséquences sociales, économiques et politiques qui en découlent.”

Le risque devient “réel” socialement lorsqu’il devient un objet de controverse publique, c’est-à-dire quand son existence a des conséquences politiques. Zillow, en retirant les scores, contribue à rendre le risque… moins réel. Mais seulement en apparence.

Références

Bankoff, G., Frerks, G., & Hilhorst, D. (Eds.). (2004).
Mapping Vulnerability: Disasters, Development and People. London: Earthscan.

Beck, U. (1992). Risikogesellschaft, Suhrkamp.

Bernstein, A., Gustafson, M. T., & Lewis, R. (2024). Flood Risk Disclosure and Housing Market Behavior. NBER Working Paper No. 31908.

Brown, C. (2025, November 30). Zillow Removes Climate Risk Scores From Home Listings. The New York Times.

Douglas, M., & Wildavsky, A. (1982). Risk and Culture: An Essay on the Selection of Technical and Environmental Dangers. University of California Press.

First Street Foundation. (2020–2025). Methodologies & Technical Documentation. First Street.

Jasanoff, S. (2003). Technologies of Humility. Nature, 450, 33.

Keenan, J. M. (2018). Sea Level Rise and Housing.
In C. Rosenzweig, W. Solecki, P. Romero-Lankao, et al. (Eds.),
Climate Change and Cities: Second Assessment Report of the Urban Climate Change Research Network (ARC3.2)
(pp. 343–385). Cambridge University Press.

Kousky, C. (2023). Understanding Disaster Insurance: New Tools for a More Resilient Future. Routledge / RFF Press.

Krauss, L. (2021). The Physics of Climate Change. Post Hill Press.

Kunreuther, H., & Meyer, R. (2017). The Ostrich Paradox: Why We Underprepare for Disasters. Wharton School Press.

Milly, P. C. D., Betancourt, J., Falkenmark, M., Hirsch, R., Kundzewicz, Z., Lettenmaier, D., & Stouffer, R. (2008). Stationarity Is Dead: Whither Water Management? Science, 319(5863), 573–574.

Redfin & National Bureau of Economic Research (NBER). (2024). Flood Risk Disclosure and Homebuyer Behavior (Working Paper).

Reghezza-Zitt, M., & Rufat, S. (2015).
Resilience Imperative: Uncertainty, Risks and Disasters. London: ISTE Press / Elsevier.

Roux, J., November, V., & Reghezza-Zitt, M. (2013). Les territoires du risque. Armand Colin.

Smiley, S. (2020). Outdated and inaccurate, FEMA flood maps fail to fully capture risk. Kinder Institute for Urban Research.

U.S. Government Accountability Office. (2022). FEMA: Flood Maps—Actions Needed to Improve Oversight and Management of Mapping Program. GAO-22-104079.

Wagner, G., & Weitzman, M. L. (2015). Climate Shock: The Economic Consequences of a Hotter Planet. Princeton University Press.

Faire ses cartes météo

Allez, un peu de scraping aujourd’hui (promis, il y aura d’autres billets d’ici la fin de la semaine sur le sujet). Je vais m’inspirer des codes de Romain de la formation Data Science pour l’Actuariat. Le but est de récupérer des données de températures, de précipitations, de vitesse de vent, journalières, en France. On commence par charger quelques librairies,

library(plyr)
library(stringr)
library(OpenStreetMap)
library(leaflet)
library(shiny)
library(rsconnect)
library(mapview)
library(png)
library(magick)
library(yaml)

On va récupérer les deux dernières années, mais comme toujours, le code peut facilement s’adapter.

annee = c(2016:2017)
mois = c('01','02','03','04','05','06','07','08','09','10','11','12')
aaaamm = sort(sub(pattern=' ',replacement = '', x=outer(annee,mois,paste)))
myCols = c(NA,NA,"NULL","NULL","NULL",NA,NA,"factor","NULL","NULL","NULL",NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL",NA,NA,NA,NA,NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL")
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (!file.exists(fichier)) {
    page = paste0('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/Archive/synop.',mois,'.csv.gz')
    download.file(page, fichier, quiet = TRUE, cacheOK = FALSE)
  }
}  
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (mois=='201601') {
    data = read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols)
  } else {
    data = rbind(data,read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols))
  }
}

On supprime une des colonnes qui ne sert a rien

data[,"X"] = NULL

et on bricole un peu (sinon on a des soucis avec minuit, qui est décalé de 24 heures)

data$temp = as.numeric(substr(data$date,1,14))
decalage_minuit = function(dateheure){
  r = dateheure
  if (dateheure-trunc(dateheure/1000000)*1000000==0){
    r = dateheure + 10000
  }
  return(r)
}
data$temp = sapply(data$temp,decalage_minuit)
data$jourheure_obs = strptime(data$temp, format = '%Y%m%d%H%M%S', 'UTC')
data$date = as.POSIXct((strptime(data$jourheure_obs,format = "%Y-%m-%d")))
data$heure = as.numeric(substr(data$jourheure_obs, start=12, stop=13))
data[,"temp"] = NULL
data$t = as.numeric(as.character(data$t))-273.15

oui, on doit convertir en degrés Celsius. On peut aussi travailler sur les vitesses de vent

data$dd = as.numeric(as.character(data$dd))
data$ff = as.numeric(as.character(data$ff))/1000*3600
data$ff = round(data$ff/5)*5

que l’on va convertir en km/h, et enfin les précipitations

data$rr1 = pmax(0,as.numeric(as.character(data$rr1)))
data$rr3 = pmax(0,as.numeric(as.character(data$rr3)))
data$rr6 = pmax(0,as.numeric(as.character(data$rr6)))
data$rr12 = pmax(0,as.numeric(as.character(data$rr12)))
data$rr24 = pmax(0,as.numeric(as.character(data$rr24)))

On peut d’ailleurs commencer par une simple visualisation de séries de températures

data_temperature = subset(data, !is.na(data$t))
temp_min = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=min)
names(temp_min)[2] = "temp_min"
temp_max = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=max)
names(temp_max)[2] = "temp_max"
temp_moy = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=mean)
names(temp_moy)[2] = "temp_moy"
temp_stat = cbind(temp_min,temp_max[],temp_moy)
rm(temp_min,temp_max,temp_moy)
temp_stat[,3] = NULL
temp_stat[,4] = NULL

avec le maximum observé, par jour, sur toutes les stations, le minimum, et la moyenne des stations

plot(x=temp_stat$date,y=temp_stat$temp_moy,'l',ylim=c(-15,40),xlab='Date',ylab='Température (°C)')
title('Températures françaises métropolitaines de 01/2016 à 12/2017')
points(x=temp_stat$date,y=temp_stat$temp_min,col="blue",'l')
points(x=temp_stat$date,y=temp_stat$temp_max,col="red",'l')

Pour faire des jolies cartes, il nous faut des informations sur les stations météo.

stations = read.table('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv',header=TRUE, sep=";",quote='"',fileEncoding="UTF-8")
stations = rename(stations, replace=c("ID"="numer_sta"))
stations = rename(stations, replace=c("Nom"="nom"))
stations = rename(stations, replace=c("Latitude"="latitude"))
stations = rename(stations, replace=c("Longitude"="longitude"))
stations = rename(stations, replace=c("Altitude"="altitude"))
data = merge(x=data,y=stations)

Pour la carte, on va se restreindre aux données en France métropolitaine

data = subset(data, data$latitude >= 39)
data = subset(data, data$latitude <= 53)
data = subset(data, data$longitude >= -5.3)
data = subset(data, data$longitude <= 9.8)
coordonnees_stations = unique(data[, c("nom","latitude","longitude")])

Commencons par visualiser le maximum observee le 19 juillet 2016

date_sel = "2016-07-19"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
temp = aggregate(x=donnees_carte_0$t,by=list(numer_sta=donnees_carte_0$numer_sta),FUN="max")
donnees_carte_0 = merge(x=donnees_carte_0,y=temp)
donnees_carte_0 = rename(donnees_carte_0, replace=c("x"="TEMP"))
donnees_carte_0$TEMP = round(donnees_carte_0$TEMP)
donnees_carte = unique(donnees_carte_0[,c('numer_sta','nom','longitude','latitude','TEMP')])
get_france = get_map(c(lon=2.25,lat=46), zoom=5, col='bw')
colfunc = colorRampPalette(c("darkblue","royalblue","cyan","lightblue","orange","red","darkred"))
liste_couleurs = data.frame(TEMP=as.numeric(seq(-49,50,1)),couleur=as.character(colfunc(100)))
donnees_carte = merge(x=donnees_carte,y=liste_couleurs)
france_tmax = ggmap(get_france) + 
  ggtitle(paste0('Températures maximales du ', date_sel_OK)) +
  scale_x_continuous(limits = c(-5, 10), expand = c(0, 0)) +
  scale_y_continuous(limits = c(41.3, 51.1), expand = c(0, 0))
france_tmax = france_tmax + geom_label(data=donnees_carte,aes(x=longitude,y=latitude,label=TEMP),col=as.character(donnees_carte$couleur))
france_tmax

On peut tenter aussi une visualisation des précipitations. On se place le 13 octobre 2016 (oui, il y a eu des inondations dans la region de Montpelier ce jour la)

date_sel = "2016-10-13"
colfunc = colorRampPalette(c('white','darkblue'))
data$precip_en_mm = data$rr3
data_precipitation = subset(data, !is.na(data$precip_en_mm))
precipitations = aggregate(x=data_precipitation$precip_en_mm,by=list(date=data_precipitation$date,station=data_precipitation$nom),FUN=sum)
names(precipitations)[3] = "pp"
precipitations$pp = pmax(precipitations$pp,0)
donnees_carte = subset(precipitations, precipitations$date==date_sel)
qt_prec = quantile(round(donnees_carte$pp),seq(0, 1, 0.05))
donnees_carte$Couleur_prec = colfunc(21)[(findInterval(round(donnees_carte$pp), qt_prec, all.inside=TRUE))]
donnees_carte = rename(donnees_carte, replace=c("station"="nom"))
donnees_carte = merge(x=donnees_carte,y=coordonnees_stations)

On peut tenter un leaflet ici (ça sera plus joli)

france = leaflet() %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = france %>% addCircles(lng=donnees_carte$longitude, lat=donnees_carte$latitude, color=donnees_carte$Couleur_prec, opacity = 1, fillColor=donnees_carte$Couleur_prec, fillOpacity = 1, radius=donnees_carte$pp*500)
france

Comme souvent, j’ai des soucis pour integrer du leaflet en wordpress, alors on se contentera dans le billet d’une copie d’ecran (mais un copier/coller du code suffit).

Et pour finir, la meme chose avec les vitesses de vent. On se place le 6 mars 2017 (la encore, la date n’est pas choisie au hasard, c’était la tempete Zeus)

date_sel = "2017-03-06"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
heure_sel = 18
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
donnees_carte_0 = subset(donnees_carte_0, donnees_carte_0$heure==heure_sel)
donnees_carte = donnees_carte_0

Pour les precipitations, des cercles plus ou moins grands suffisaient. La, on va mettre des fleches

logo = image_trim(image_read("https://image.freepik.com/icones-gratuites/fleche-noire-vers-le-haut_318-30934.jpg"))
if (file.exists('flecheblanche90.png')==FALSE){
  for (angle in seq(10,360, by=10)) {
    temp = image_rotate(logo,angle)
    image_write(temp,paste0('flecheblanche',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='green')
    image_write(temp,paste0('flecheverte',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='yellow')
    image_write(temp,paste0('flechejaune',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='orange')
    image_write(temp,paste0('flecheorange',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='red')
    image_write(temp,paste0('flecherouge',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='black')
    image_write(temp,paste0('flechenoire',angle,'.png'))
  }
}
arrowIcons = icons(
  iconUrl = ifelse(donnees_carte$dd == 0,
            ifelse(donnees_carte$ff<39,
                    'flecheverte360.png',
            ifelse(donnees_carte$ff<79,
                    'flechejaune360.png',
            ifelse(donnees_carte$ff<119,
                    'flecheorange360.png',
            ifelse(donnees_carte$ff<159,
                    'flecherouge360.png',
            ifelse(donnees_carte$ff>159,
                    'flechenoire360.png',
                    'flecheblanche360.png'))))),
            ifelse(donnees_carte$ff<39,
                   paste0('flecheverte',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<79,
                    paste0('flechejaune',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<119,
                    paste0('flecheorange',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<159,
                    paste0('flecherouge',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff>159,
                    paste0('flechenoire',donnees_carte$dd,'.png'),
           paste0('flecheblanche',donnees_carte$dd,'.png'))))))),
  iconWidth = 30, iconHeight = 30)

Cette fois, on est bon ! On peut faire la carte

france = leaflet(options = leafletOptions()) %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = addMarkers(map = france,lng=donnees_carte$longitude,lat=donnees_carte$latitude,icon=arrowIcons,popup=as.character(donnees_carte$ff))
france

Accidents de voiture et de deux roues à Paris

Il y a quelques mois, j’avais commencé à reprendre des codes que j’avais vu passer dans les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. J’avais commencé avec un billet sur sur cartographier le vote pour le Brexit (qui avait été repris – et amélioré – sur le site des voisins de rgeomatic). Aujourd’hui, en rangeant un peu mon ordinateur, je suis tombé sur un code de Quentin sur les acccidents de la route à Paris.

Comme toujours on charge quelques packages

require(downloader)
require(ggmap)
require(stringr)
options(encoding = "UTF-8")

et on récupère les données, en ligne sur https://opendata.paris.fr

path="https://opendata.paris.fr/explore/dataset/accidentologie/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true"
dest="test"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
obj <- read.csv(file = loc, stringsAsFactor = FALSE,sep=";")
file.remove(dest)

On peut commencer par regarder quand on eu lieu les accidents, quelle jour et à quelle heure. On va alors agréger les données

t=data.frame(obj,counter=1)
parJour=aggregate(t$counter, by=list(format(as.Date(t$Date),"%A")), sum)

et créer une base qui contient toutes les informations

vec_dat=data.frame();
vec_dat[ (1:3),1]=parJour$Group.1[3:5]
vec_dat[ (1:3),2]=parJour$x[3:5]
vec_dat[ 4,1]=parJour$Group.1[2]
vec_dat[ 4,2]=parJour$x[2]
vec_dat[ 5,1]=parJour$Group.1[7]
vec_dat[ 5,2]=parJour$x[7]
vec_dat[ 6,1]=parJour$Group.1[6]
vec_dat[ 6,2]=parJour$x[6]
vec_dat[ 7,1]=parJour$Group.1[1]
vec_dat[ 7,2]=parJour$x[1]

on va rajouter un nom, pour que ce soit plus clair

names(vec_dat)[1]="Jour";names(vec_dat)[2]="Nbsinistre"

On choisit éventuellement une palette de couleurs

Palette<-colorRampPalette(c("lightblue","darkblue"))

et on fait un barplot, pour compter le nombre d’accidents, par jour de la semaine, sur Paris

barplot(vec_dat$Nbsinistre,names.arg=vec_dat$Jour,main="Répartition du nombre d'accidents par jour de la semaine ", col=Palette(10), xlab="jour de la semaine", ylab="nombre d'accidents")

On peut ensuite regarder par heure de la journée (en données brutes, sans lissage)

parHeure=aggregate(t$counter, by=list(format(strptime(t$Heure,format="%H:%M"),"%H")), sum)
plot(parHeure,type="l", col=c("blue"), main=" Répartition du nombre d'accidents par heure de la journée", xlab="heure de la journee", ylab="nombre d'accidents")

Passons maintenant à de la géolocalisation

t2 = subset(t, Coordonnées!= "")
coord=data.frame(t(matrix(as.numeric(unlist(strsplit(t2$Coordonnées,"[,]"))),nrow=2)))
t2$lat=coord[,1]
t2$lon=coord[,2]
t2 = subset(t2, (
(48.81 <= lat) & (lat<=48.90) & (2.25<=lon) &(lon<=2.43) ) )

On peut regarder sur une carte où ont eu lieu les accidents

map<-get_map(location = "Paris", zoom=12, maptype="roadmap", color="bw")
vis <- ggmap(map) + geom_point(data = t2, aes(x=lon, y=lat, color=1)) + theme( legend.position="none") + labs(title="Cartographie des accidents à Paris")
vis

On peut tenter une densité (brute, non corrigée d’une densité de routes, ou de population)

overlay<-stat_density2d(data=t2, aes(x=lon, y=lat, fill=..level.., alpha=..level..), contour=T, n=100, geom="polygon")
densi<- ggmap(map) + overlay + scale_fill_gradient("Accident Density") +scale_alpha(range=c(0.4,0.75),guide=FALSE) + guides(fill=guide_colorbar(barwidth=1.5, barheight = 10))
densi = densi + labs(title="Densite des accidents à Paris")
densi

et non-corrigée ici d’effets de bords (ce qui donne l’illusion qu’il y a plus d’accidents intra-muros que sur le boulevard périphérique par exemple). On peut également regarder en fonction de variables explicatives, comme le type de véhicules impliqués, ou la gravité

densi_byGravUs<- densi + facet_wrap(~ Usager.1.Grav, nrow=2)
densi_byGravUs = densi_byGravUs + labs(title="Densité des accidents à Paris par gravité d'accident")
densi_byGravUs

Bref, on peut faire facilement plein de chose à partir de ces données.

Analyse des Données et Cartes

Vendredi, on continuera en cours la classification (non supervisée) et en particulier les méthodes hiérarchiques, et les arbres. On va utiliser la base suivante, avec les résultats des élections présidentielles (premier tour) de 2012,

elections2012 <- read.csv("http://komodo.regardscitoyens.org/public/presidentielles2012/resultats/resultats_departements_final_T1.csv",
sep=";",header=TRUE,dec=",")

La base qu’on garde contient juste les voix exprimées, en pourcentage

voix <- which(substr(names(elections2012),1,12)=="X..Voix.Exp.")
X <- as.matrix(elections2012[,voix])
colnames(X) <- substr(names(elections2012)[voix],13,nchar(names(elections2012)[voix]))
rownames(X) <- elections2012[,3]

Parfois, il convient de normaliser les données. Ici, sur les département, ça n’est pas nécessaire. Par contre, sur les candidats, ça le serait (comme pour l’ACP en fait).  Pour visualiser les distances (entre départements), on peut utiliser

heatmap(X)

On peut faire une CAH, sur notre matrice de distance entre rangs

cah <- hclust(dist(X))
plot(cah,cex=.6)

et si on souhaite garder 5 groupes (par exemple) on utilise

rect.hclust(cah,k=5)
groups.5 <- cutree(cah,5)

On peut aussi visualiser les classes avec la fonction

library(dendroextras)
plot(colour_clusters(cah,k=5))

Pour comprendre qui se trouve dans nos groupes, on peut utiliser

aggregate(X,list(groups.5),mean)

qui va nous renvoyer les “votes  moyens” pour chaque candidat, dans chaque groupe. Maintenant, pour en finir avec la visualisation, on peut utiliser le code suivant (un peu long, certes) qui permettra de visualiser sur un carte les différents groupes

carte_classe <- function(groupes){
library(stringr)
elections2012$dep <- elections2012$Libellé.du.département
elections2012$dep <- tolower(elections2012$dep)
elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
library(maps)
france<-map(database="france")
france$dep <- france$names
france$dep <- tolower(france$dep)
france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
corresp_noms <- elections2012[, c("Libellé.du.département", "dep")]
corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
col2001<-groupes+1
names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms$Libellé.du.département)]
color <- col2001[match(france$dep, names(col2001))]
map(database="france", fill=TRUE, col=color)
}

carte_classe(groups.5)

Avec ces fonctions, on devrait pouvoir tester différentes méthodes pour constituer des groupes.

Le Vote par Procuration en France

La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.

Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.

  • Modèle de régression, et recherche de variables explications

Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que  où  est le nombre de procurations dans le bureau de vote , et où  est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que , la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

Si on suppose que  est fonction  le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

il est légitime de tenter une régression logistique,

voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,

Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

et, pour la régression logistique standard (linéaire)

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration).

On peut faire la même chose sur le taux de logement HLM, dans le quartier,

Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de  laisse à désirer (toute suggestion est la bienvenue !).

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

 

On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire

(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre)

Cette fois, on a une sortie de régression plus classique,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions).

  • Visualiser les bureaux de vote à Paris

Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net)

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

Reste à conclure, en faisant une visualisation graphique de nos données

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

Dans ce cas, la visualisation des résidus donne

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.

ACT2121, éléments de correction

La plupart des calculs pouvant se faire sans calculs (trop) complexes avec une calculatrice, je vais revenir sur un exercice (l’exercice 11 de la seconde série) pour proposer une correction.

Le petit Nestor collectionne les cartes de joueurs de Baseball dans les paquets de gommes à mâcher. Il y a en tout 20 cartes différentes (réparties aléatoirement, une par paquet). Combien de paquets de gommes Nestor devrait-il s’attendre à avoir à acheter pour obtenir la collection complète ? 

La bonne stratégie était d’ordonner les cartes par ordre d’apparition, pour la première fois, et de modéliser le nombre de paquets entre deux premières apparitions, comme indiqué sur le dessin ci-dessous


Si  désigne le nombre total de paquets à acheter, on note  le nombre de paquets à acheter, entre l’apparition de la ème nouvelle carte, et la ème (avec la convention que vaut 1, i.e. à l’achat du premier paquet, on a notre première carte). On notera que
https://latex.codecogs.com/gif.latex?M_i\sim\mathcal{G}((20-(i-1))/20)
et à partir de là, les calculs sont simples, puisque

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=mathbb{E}(M_1+M_2+cdots+M_{20})

soit (par linéarité de l’espérance)

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=mathbb{E}(M_1)+mathbb{E}(M_2)+cdots+mathbb{E}(M_{20})

i.e.

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=left(%201+%20frac{20}{19}%20+%20frac{20}{18}%20+%20frac{20}{17}%20+%20...%20+%20frac{20}{2}+frac{20}{1}%20<br /> ight)

ou encore

https://latex.codecogs.com/gif.latex?mathbb{E}(N)=20left(%20frac{1}{20}+%20frac{1}{19}%20+%20frac{1}{18}%20+%20frac{1}{17}%20+%20...%20+%20frac{1}{2}+1%20<br /> ight)

La bonne réponse était alors (il faut sommer les vingt termes)

> sum(20/(1:20)) [1] 71.95479

Mais cette sommation de vingt termes n’est pas triviale, aussi, en cours, j’avais suggéré que

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)
> log(20)*20 [1] 59.91465

qui diffère du résultat numérique car il manquait la constante d’Euler

https://latex.codecogs.com/gif.latex?gamma%20=%20lim_{n%20<br /> ightarrow%20infty%20}%20left(%201+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20-%20ln(n)%20<br /> ight)

i.e.

https://latex.codecogs.com/gif.latex?1+%20frac{1}{2}%20+%20frac{1}{3}%20+%20frac{1}{4}%20+%20...%20+%20frac{1}{n}%20sim%20ln(n)+gamma

Numériquement, on obtient, en prenant pour https://latex.codecogs.com/gif.latex?gamma 0.57721 (cf Google)

> (log(20)+ 0.57721)*20 [1] 71.45885