# 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 &gt;= 39) data = subset(data, data$latitude &lt;= 53) data = subset(data, data$longitude &gt;= -5.3) data = subset(data, data$longitude &lt;= 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() %&gt;% addTiles() %&gt;% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57) france = france %&gt;% 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&lt;39, 'flecheverte360.png', ifelse(donnees_carte$ff&lt;79, 'flechejaune360.png', ifelse(donnees_carte$ff&lt;119, 'flecheorange360.png', ifelse(donnees_carte$ff&lt;159, 'flecherouge360.png', ifelse(donnees_carte$ff&gt;159, 'flechenoire360.png', 'flecheblanche360.png'))))), ifelse(donnees_carte$ff&lt;39, paste0('flecheverte',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff&lt;79, paste0('flechejaune',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff&lt;119, paste0('flecheorange',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff&lt;159, paste0('flecherouge',donnees_carte$dd,'.png'), ifelse(donnees_carte$ff&gt;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()) %&gt;% addTiles() %&gt;% 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 $P_i\sim\mathcal{B}(N_i,p_i)$ où $P_i$ est le nombre de procurations dans le bureau de vote $i$, et où $N_i$ est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que $p_i$, 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=";")
> 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 $p_i$ est fonction $X_i$ 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,

$p_i=h(X_i)=\frac{\exp[\beta_0+\beta_1 X_i]}{1+\exp[\beta_0+\beta_1 X_i]}$

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

$p_i=\tilde h(X_i)=\frac{\exp[s(X_i)]}{1+\exp[s(X_i)]}$

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 $\beta_1$ 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 $\frac{P_i}{N_i}=\beta_0+\beta_1 X_i+\varepsilon_i$(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 $N$ désigne le nombre total de paquets à acheter, on note $M_i$ le nombre de paquets à acheter, entre l’apparition de la $i-1$ème nouvelle carte, et la $i$ème (avec la convention que $M_1$vaut 1, i.e. à l’achat du premier paquet, on a notre première carte). On notera que
$M_i\sim\mathcal{G}((20-(i-1))/20)$
et à partir de là, les calculs sont simples, puisque

$\mathbb{E}(N)=\mathbb{E}(M_1+M_2+\cdots+M_{20})$

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

$\mathbb{E}(N)=\mathbb{E}(M_1)+\mathbb{E}(M_2)+\cdots+\mathbb{E}(M_{20})$

i.e.

$\mathbb{E}(N)=\left( 1+ \frac{20}{19} + \frac{20}{18} + \frac{20}{17} + ... + \frac{20}{2}+\frac{20}{1} \right)$

ou encore

$\mathbb{E}(N)=20\left( \frac{1}{20}+ \frac{1}{19} + \frac{1}{18} + \frac{1}{17} + ... + \frac{1}{2}+1 \right)$

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

$1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} \sim \ln(n)$
> log(20)*20 [1] 59.91465

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

$\gamma = \lim_{n \rightarrow \infty } \left( 1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} - \ln(n) \right)$

i.e.

$1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} \sim \ln(n)+\gamma$

Numériquement, on obtient, en prenant pour $\gamma$ 0.57721 (cf Google)

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