Monthly Archives: February 2015

Publications, Research

Feldstein and Horioka, Mother of all puzzles

Leave a comment

The joint paper  The “Mother of all puzzles” at thirty: A Meta-Analysis is finally out.

This paper provides a meta-analysis of 1651 point estimates of Feldstein and Horioka saving retention coefficient from 49 peer-reviewed papers published over three decades. We get two main results. First, correcting for publication bias, we find a consistent underlying coefficient lying between 0.56 and 0.67 for studies using the original paper. Second, heterogeneity reported in the estimates of the Feldstein and Horioka can be explained by a few main factors. In particular, we find evidence that the saving retention coefficient is systematically underestimated with models written in first difference, models using the saving ratio or the current account ratio as the dependent variable instead of the investment ratio, and models including indicators of the public deficit or indicators of the country size as additional explanatory variables.

Maps

John Snow, and Google Maps

7 Comments

In my previous post, I discussed how to use OpenStreetMaps (and standard plotting functions of R) to visualize John Snow’s dataset. But it is also possible to use Google Maps (and ggplot2 types of graphs).

library(ggmap)
get_london <- get_map(c(-.137,51.513), zoom=17)
london <- ggmap(get_london)

Again, the tricky part comes from the fact that the coordinate representation system, here, is not the same as the one used on Robin Wilson’s blog.

> library(maptools)
> setwd("/cholera/")
> deaths <- readShapePoints("Cholera_Deaths")
> head(deaths@coords)
coords.x1 coords.x2
0  529308.7  181031.4
1  529312.2  181025.2
2  529314.4  181020.3
3  529317.4  181014.3
4  529320.7  181007.9
5  529336.7  181006.0
> X <- deaths@coords

or, use d X_deaths.RData. So now, we have to change it

df_deaths <- data.frame(X)
library(sp)
library(rgdal)
coordinates(df_deaths)=~coords.x1+coords.x2
proj4string(df_deaths)=CRS("+init=epsg:27700") 
df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))

Here, we have the same coordinate system as the one used in Google Maps. Now, we can add a layer, with the points,

london + geom_point(aes(x=coords.x1, y=coords.x2),data=data.frame(df_deaths@coords),col="red")

Again, it is possible to add the density, as an additional layer,

london + geom_point(aes(x=coords.x1, y=coords.x2), 
data=data.frame(df_deaths@coords),col="red")+
geom_density2d(data = data.frame(df_deaths@coords), 
aes(x = coords.x1, y=coords.x2), size = 0.3) + 
stat_density2d(data = data.frame(df_deaths@coords), 
aes(x = coords.x1, y=coords.x2,fill = ..level.., alpha = ..level..),size = 0.01, bins = 16, geom = "polygon") + scale_fill_gradient(low = "green", high = "red",guide = FALSE) + 
scale_alpha(range = c(0, 0.3), guide = FALSE)

 

Maps

John Snow, and OpenStreetMap

11 Comments

While I was working for a training on data visualization, I wanted to get a nice visual for John Snow’s cholera dataset. This dataset can actually be found in a great package of famous historical datasets.

library(HistData)
data(Snow.deaths)
data(Snow.streets)

One can easily visualize the deaths, on a simplified map, with the streets (here simple grey segments, see Vincent Arel-Bundock’s post)

plot(Snow.deaths[,c("x","y")], col="red", pch=19, cex=.7,xlab="", ylab="", xlim=c(3,20), ylim=c(3,20))
slist <- split(Snow.streets[,c("x","y")],as.factor(Snow.streets[,"street"]))
invisible(lapply(slist, lines, col="grey"))

Continue reading John Snow, and OpenStreetMap →

Datamining, M1-Data Science, Statistics

Visualizing Clusters

12 Comments

Consider the following dataset, with (only) ten points

x=c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
y=c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
plot(x,y,pch=19,cex=2)

We want to get – say – two clusters. Or more specifically, two sets of observations, each of them sharing some similarities.

Since the number of observations is rather small, it is actually possible to get an exhaustive list of all partitions, and to minimize some criteria, such as the within variance. Given a vector with clusters, we compute the within variance using

within_var = function(I){
I0=which(I==0)
I1=which(I==1)
xbar0=mean(x[I0])
xbar1=mean(x[I1])
ybar0=mean(y[I0])
ybar1=mean(y[I1])
w=sum(I0)*sum( (x[I0]-xbar0)^2+(y[I0]-ybar0)^2 )+
  sum(I1)*sum( (x[I1]-xbar1)^2+(y[I1]-ybar1)^2 )
return(c(I,w))
}

Then, to compute all possible partitions, use

base2=function(z,n=10){
  Base.b=rep(0,n)
  ndigits=(floor(logb(z, base=2))+1)
  for(i in 1:ndigits){
    Base.b[ n-i+1]=(z%%2)
    z=(z%/%2)}
  return(Base.b)}
L=function(x) within_var(base2(x))
S=sapply(1:(2^10),L)

The cluster indices at the mimimum is here

I=S[1:n,which.min(S[n+1,])]

To visualize those clusters, use

cluster_viz = function(indices){
library(RColorBrewer)
CL2palette=rev(brewer.pal(n = 9, name = "RdYlBu"))
CL2f=CL2palette[c(1,9)]
plot(x,y,pch=19,xlab="",ylab="",xlim=0:1,ylim=0:1,cex=2,col=CL2f[1+I])
CL2c=CL2palette[c(3,7)]
I0=which(indices==0)
I1=which(indices==1)
xbar0=mean(x[I0])
xbar1=mean(x[I1])
ybar0=mean(y[I0])
ybar1=mean(y[I1])
segments(x[I0],y[I0],xbar0,ybar0,col=CL2c[1])
segments(x[I1],y[I1],xbar1,ybar1,col=CL2c[2])
points(xbar0,ybar0,pch=19,cex=1.5,col=CL2c[1])
points(xbar1,ybar1,pch=19,cex=1.5,col=CL2c[2])}

and then, simply

cluster_viz(I)

But that was possible only because https://latex.codecogs.com/gif.latex?n is not to large (since the total number of scenarios – with only 2 clusters – is https://latex.codecogs.com/gif.latex?2^n, or https://latex.codecogs.com/gif.latex?2^{n-1} if we changes zeros in ones).

Continue reading Visualizing Clusters →

Computer, M1-Data Science

k-means clustering and Voronoi sets

5 Comments

In the context of https://latex.codecogs.com/gif.latex?k-means, we want to partition the space of our observations into https://latex.codecogs.com/gif.latex?k classes. each observation belongs to the cluster with the nearest mean. Here “nearest” is in the sense of some norm, usually the https://latex.codecogs.com/gif.latex?\ell_2 (Euclidean) norm.

Consider the case where we have 2 classes. The means being respectively the 2 black dots. If we partition based on the nearest mean, with the https://latex.codecogs.com/gif.latex?\ell_2 (Euclidean) norm we get the graph on the left, and with the https://latex.codecogs.com/gif.latex?\ell_1 (Manhattan) norm, the one on the right,

Points in the red region are closer to the mean in the upper part, while points in the blue region are closer to the mean in the lower part. Here, we will always use the standard https://latex.codecogs.com/gif.latex?\ell_2 (Euclidean) norm. Note that the graph above is related to Voronoi diagrams (or Voronoy, from Вороний in Ukrainian, or Вороно́й in Russian) with 2 points, the 2 means.

Continue reading k-means clustering and Voronoi sets →

M2-inequalities

Readings on Inequalities

Leave a comment

Here are some reading for the course on inequalities,

M2-inequalities

Inequalities and Quantile Regression

2 Comments

In the course on inequality measure, we’ve seen how to compute various (standard) inequality indices, based on some sample of incomes (that can be binned, in various categories). On Thursday, we discussed the fact that incomes can be related to different variables (e.g. experience), and that comparing income inequalities between coutries can be biased, if they have very different age structures.

So we’ve seen quantile regressions. I can mention some old slides (used in a crash course at McGill three years ago)., as well as a more technical discussion on ties, and non-unicity of the regression line.

In order to illustrate, consider  the following dataset

> salary <- read.table("http://data.princeton.edu/wws509/datasets/salary.dat",header=TRUE)
> plot(salary$yd,salary$sl)
> abline(lm(sl~yd,data=salary),col="blue")

We have here the stndard regression line, obtained using ordinary least squares. Here we have the expected income given the experience. But we can also use a quantile regression,

https://latex.codecogs.com/gif.latex?Q_\tau(Y\vert\boldsymbol{X})=\boldsymbol{X}^{\text{\sffamily%20T}}\boldsymbol{\beta}

> library(quantreg)
> Q10 <- rq(sl~yd,data=salary,tau=.1)
> Q90 <- rq(sl~yd,data=salary,tau=.9)
> abline(Q10,col="red")
> abline(Q90,col="purple")

A classical tool to describe inequalities is the ratio of the 90% quantile over the 10% quantile (among so many others,

> ratio9010 = function(age){
+   predict(Q90,newdata=data.frame(yd=age))/
+   predict(Q10,newdata=data.frame(yd=age))
+ }

For instance, among people with 5 years of experience, there is an inequality index of

> ratio9010(5)
1.401749

while for people with 30 years of experience, it would be

> ratio9010(30)
1.9488

If we plot the evolution of this 90-10 ratio, as a function of the experience, we get the following increasing trend

> A=0:30
> plot(A,Vectorize(ratio9010)(A),type="l",ylab="90-10 quantile ratio")

So clearly, comparing inequalitis ceteris paribus between two groups, should be performed carefully, and probably including some covariates.

M1-Data Science

Analyse des Données et Cartes

Leave a comment

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.

M2-inequalities

Inequality, Poverty and Welfare

Leave a comment

For the fourth cours on Inequalities, we will get back on the quantile regression, and discuss welfare functions as well as poverty indices. Slides are now online

To illustrate, we will use the following datasets

uk88 <- read.csv("http://www.vcharite.univ-mrs.fr/pp/lubrano/cours/fes88.csv",sep=";",header=FALSE)$V1
uk92 <- read.csv("http://www.vcharite.univ-mrs.fr/pp/lubrano/cours/fes92.csv",sep=";",header=FALSE)$V1
uk96 <- read.csv("http://www.vcharite.univ-mrs.fr/pp/lubrano/cours/fes96.csv",sep=";",header=FALSE)$V1
cpi <- c(421.7, 546.4, 602.4)
y88 <- uk88/cpi[1]
y92 <- uk92/cpi[2]
y96 <- uk96/cpi[3]

and for the part on applications of quantile regression

salary <- read.table("http://data.princeton.edu/wws509/datasets/salary.dat",header=TRUE)
Actuarial science, Big Data, Publications

Big Data : passer d’une analyse de corrélation à une interprétation causale

2 Comments

Le rôle d’un actuaire dans une compagnie d’assurance est assez souvent d’estimer la probabilité qu’un événement survienne, ou ses possibles conséquences financières, et fonction de variables dites « explicatives ». On voit en effet que certaines variables sont « statistiquement corrélées » avec la survenance d’un accident dans l’année, mais prétendre que l’on dispose d’une « explication » est un peu dangereux. Ce sont pourtant les interrogations qui étaient soulevées lors des débats sur l’utilisation du sexe en tarification : si les femmes ont moins d’accidents, en moyenne, que les hommes, pourquoi ne pas utiliser cette variable en tarification ? Le problème est que s’arrêter à une étude des corrélations ne permet pas de vraiment comprendre ce qui se cache vraiment derrière un phénomène. Dans un précédant billet, on notait que de telles études pouvaient conduire à des interprétations paradoxales, et erronées.
Comme le notait Dubuisson (2008), « comme le reconnaissent les actuaires, la mise en évidence d’un lien causal entre le critère choisi et la variation de la sinistralité s’apparente à la quête du Graal ». Les données massives, big data, permettent peut être d’avoir accès à davantage d’information, et de mieux comprendre ce qui peut causer un risque. Un exemple classique est un paradoxe qui a été long à comprendre par les épidémiologistes, sur la mortalité infantile, le poids des bébés et le tabagisme de la mère. Nous allons reprendre cet exemple ici et analyser comment l’utilisation de données massives a permis de mieux comprendre ce qui pouvait réellement causer une surmortalité infantile.

Continue reading Big Data : passer d’une analyse de corrélation à une interprétation causale →

M1-Time Series

Mensualiser une série hebdomadaire

Leave a comment

Pour le devoir de série temporelles, les données fournies par Google trend sont hebdomadaires, ce qui peut rendre la modélisation compliquée. Comme on l’a évoqué en cours, il peut être plus simple de travailler sur des données mensuelles. La petite fonction suivante permet de transformer les données pour avoir des données mensuelles (avec des moyennes par mois, ce qui fait qu’un mois de 28 jours et un mois de 31 jours sont comparables). Histoire d’illustrer, prenons un mot au hasard, disons…. épilation. On commence par sauver le fichier csv, et on le lit sous R,

EPILATION=read.table("/home/charpentier/report-epilation.csv", skip=4,header=TRUE,sep=",",nrows=426)

La petite fonction suivant va nous aider à convertir la base en une série temporelle mensuelle,

H2M=function(BASE){ 
X=BASE[,2] 
date1=substr(as.character(BASE$Semaine),1,10) 
date2=substr(as.character(BASE$Semaine),14,23) 
D1=as.Date(date1,"%Y-%m-%d") 
D2=as.Date(date2,"%Y-%m-%d") 
vm=vy=N=NA 
for(t in 1:length(D1)){ 
mois=seq(D1[t],D2[t],length=7) 
vm=c(vm,as.POSIXlt(mois)$mon+1) 
vy=c(vy,as.POSIXlt(mois)$year+1900) 
N=c(N,rep(X[t],7))} N=N[-1]; 
vm=vm[-1]; 
vy=vy[-1] 
YM=vy*100+vm 
Z=tapply(N,as.factor(YM),mean) 
Zts=ts(Z,start=c(2004,1),frequency=12) 
return(Zts)}

Si on travaille sur la série hebdomadaire, on a la série suivante

hebdo=ts(EPILATION$épilation,start=2004, frequency=52)

La fonction d’autocorrélation est alors

acf(hebdo,160)

Maintenant, on peut travailler sur les données mensuelles. On utilise alors

mensuel=H2M(EPILATION)

Le graphique est alors le suivant

La fonction d’autocorrélation est cette fois la suivante

acf(mensuel,40)

On retrouve le comportement cyclique, avec une saisonnalité annuelle.