Category Archives: Statistics

Using “home made” statistics

Since I am still at the Fields Institute in Toronto, enjoying a workshop on Impacts of Climate Change on Economics, Finance, and Insurance, I wanted to share some experience, from this summer. After three years of lockdown because of the covid situation, the family has been able to travel, and we went to France, so that our kids could see their grand-parents (the last visit was a long time ago). And it was hot, very hot. While I was chating with my dad, about the weather, and told me that he had a lot of connected devices in the house, including measures of the temperature. One of the device was in a place where nothing did not really change over time. So I thought it could be sufficient to get robust data. My goal was to see how the popular IPCC graph was on real data

When I got the data, I did plot them, and did compare the distribution back in 2012, and in 2022 (or to be honest, half 2021-half 2022). As for the IPCC graph, I assume a Gaussian distribution.

As expected, there is a clear shift to the right (that is “climate change”). But the most scary part, was actually the linear trend,

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -637.30455   80.44650  -7.922 3.01e-15 ***
x              0.32273    0.03988   8.092 7.72e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

with a slope of 0.322, meaning that the average temperature is increasing by 0.322 degrees per year ! That is more than 3°C over the past ten years ! Let me write it again : in a house, +3°C on average over the past ten years.

I thought there were some issues with the data. So I tried to collect some official data, and since there were no official records in their village, I did use the data from Lyon (which is 80 kilometers from their house).

The shift on the right is confirmed here, but unfortuntely, I could not get data after 2020.  Now

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -567.27953   96.26577  -5.893 4.17e-09 ***
x              0.28803    0.04776   6.031 1.81e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

And here again, I have a slope close to 0.3. So again, mainland, about +3°C over the past 10 years was observed. You might not find that scary, but I do think that it is scary !

Statistique bayésienne, data sciences et nouveaux risques

En cette rentrée, l’Institut des Actuaires lance un cycle de conférences sur le thème Statistique bayésienne, data sciences et nouveaux risques. Comme ils m’ont fait le plaisir et l’honneur d’introduire ce cycle, je ferais un exposé introductif jeudi 22 septembre, en fin d’après midi. Mes slides sont en ligne (je déconseille de les imprimer, j’ai mis des animations qui s’étalent sur plusieurs slides, je conseille plutôt cette version).

Interprétabilité et explicabilité (formalisé) des modèles prédictifs

Dans Confessiones, Saint Augustin écrivait

quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio

que l’on traduit

Qu’est-ce-que le temps ? Si personne ne me le demande, je le sais. Si je veux l’expliquer à qui me le demande, je ne le sais plus.

Pour aller un peu plus loin (car souvent, si on nous demande d’expliquer, on a quelques idées), dans Une étude en rouge de Sir Arthur Conan Doyle, paru en 1887, on a l’échange suivant, entre Sherlock Holmes et le docteur Watson,

– Je me demande ce que cherche ce type là-bas, demandai-je,
désignant un grand individu habillé simplement qui suivait
l’autre côté de la rue, en examinant anxieusement les numéros.
Il tenait à la main une grande enveloppe bleue et, de toute
évidence, portait un message.
– Vous parlez de ce sergent d’infanterie de marine ? dit Sherlock Holmes.

puis, comme il s’avère que la personne est effectivement sergent dans la marine (tout comme un autre personnage de l’histoire, un certain Arthur Charpentier), le docteur Holmes lui demande une explication, il veut savoir comment il est arrivé à cette conclusion

« Comment diable avez-vous pu deviner cela ? demandai-je.
– Deviner quoi ? fit-il sans aménité.
– Eh bien, qu’il était un sergent de marine en retraite ?
– Je n’ai pas de temps à perdre en bagatelles ! répondit-il
avec brusquerie avant d’ajouter dans un sourire : excusez ma rudesse ! Vous avez rompu le fil de mes pensées. Mais c’est peut-être aussi bien. Ainsi donc vous ne voyiez pas que cet homme était un sergent de marine ?
Non, certainement pas !
– Décidément, l’explication de ma méthode me coûte plus que son application ! Si l’on vous demandait de prouver que deux et deux font quatre, vous seriez peut-être embarrassé ; et cependant, vous êtes sûr qu’il en est ainsi. Malgré la largeur de la rue, j’avais pu voir une grosse ancre bleue tatouée sur le dos de la main du gaillard. Cela sentait la mer. Il avait la démarche militaire et les favoris réglementaires ; c’était, à n’en pas douter, un marin. Il avait un certain air de commandement et d’importance. Rappelez-vous son port de tête et le balancement de sa canne ! En outre, son visage annonçait un homme d’âge moyen, sérieux, respectable. Tous ces détails m’ont amené à penser qu’il était sergent.
– C’est merveilleux ! m’écriai-je.
– Peuh ! L’enfance de l’art ! dit Holmes, mais d’un air qui me
parut trahir sa satisfaction devant ma surprise et mon admiration manifestes.

(pour être honnête, c’est Liu Cixin qui en parle dans Le problème à trois corps). Pour l’anecdote, c’est la première histoire du couple Holmes-Watson, qui introduit la méthode de travail de Sherlock Holmes. Pour ceux qui sont familier avec les nouvelles, cette approche narrative sera largement reprise par la suite: Sherlock Holmes énonce un fait, le docteur Watson est étonné et demande une explication, et Sherlock Holmes explique, point par point, comment il est arrivé à cette conclusion. C’est un peu cette approche qu’on tente de mettre en place quand on va construire un modèle prédictif : sur la base des données du Titanic, si on prédit que telle personne va mourir, et que telle autre va survivre, on veut comprendre pourquoi le modèle arrive à cette conclusion.

Continue reading Interprétabilité et explicabilité (formalisé) des modèles prédictifs

Interpretability and explainability of predictive models

In 400 AD, in his Confessiones, Augustine wrote

quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio

that can be translated as

What then is time? If no one asks me, I know what it is. If I wish to explain it to him who asks, I do not know.

To go a little further (because often, if we are asked to explain, we have some ideas), in A Study in Scarlet by Sir Arthur Conan Doyle, published in 1887, we have the following exchange, between Sherlock Holmes and Doctor Watson

– “I wonder what that fellow is looking for?” I asked, pointing to a stalwart, plainly-dressed individual who was walking slowly down the other side of the street, looking anxiously at the numbers. He had a large blue envelope in his hand, and was evidently the bearer of a message.
– “You mean the retired sergeant of Marines,” said Sherlock Holmes.

then, as it turns out that the person is indeed a sergeant in the navy (as is another character in the story, someone named Arthur Charpentier), Dr. Holmes asks him for an explanation, he wants to know how he arrived at this conclusion

– “How in the world did you deduce that?” I asked.
“Deduce what?” said he, petulantly.
“Why, that he was a retired sergeant of Marines.”
“I have no time for trifles,” he answered, brusquely; then with a smile, “Excuse my rudeness. You broke the thread of my thoughts; but perhaps it is as well. So you actually were not able to see that that man was a sergeant of Marines?”
“No, indeed.”
– “It was easier to know it than to explain why I knew it. If you were asked to prove that two and two made four, you might find some difficulty, and yet you are quite sure of the fact. Even across the street I could see a great blue anchor tattooed on the back of the fellow’s hand. That smacked of the sea. He had a military carriage, however, and regulation side whiskers. There we have the marine. He was a man with some amount of self-importance and a certain air of command. You must have observed the way in which he held his head and swung his cane. A steady, respectable, middle-aged man, too, on the face of him – all facts which led me to believe that he had been a sergeant.”

(to be honest, it is Liu Cixin who talks about it in The Three-Body Problem). For the record, this is the first story of the Holmes-Watson couple, which introduces Sherlock Holmes’ working method. For those who are familiar with the short stories, this narrative approach will be widely used thereafter: Sherlock Holmes states a fact, Dr. Watson is astonished and asks for an explanation, and Sherlock Holmes explains, point by point, how he arrived at this conclusion. This is a bit like the approach we try to implement when we build a predictive model: on the basis of the Titanic data, if we predict that such and such a person will die, and that such and such a person will survive, we want to understand why the model arrives at this conclusion.
Continue reading Interpretability and explainability of predictive models

Report des voix entre les deux tours aux élections présidentielles en France

Allez, on va profiter des élections présidentielles, à deux tours, en France, pour jouer un peu, maintenant que les données commencent à être disponibles. On peut avoir les données par département (une centaine) ou par commune (36000)

# département
loc_dpt_1 = "https://www.data.gouv.fr/fr/datasets/r/312b99d0-27a7-429e-be8a-b930975545e9"
loc_dpt_2 = "https://www.data.gouv.fr/fr/datasets/r/c7b72e78-5e78-4976-995f-30533f0a35d3"
# commune
loc_com_1 = "https://www.data.gouv.fr/fr/datasets/r/54782507-e795-4f9d-aa70-ed06feba22e3"
loc_com_2 = "https://www.data.gouv.fr/fr/datasets/r/1601b027-48d1-42b1-96f1-b0fb49f55aa3"
download.file(loc_dpt_1,"dpt1.csv")
download.file(loc_dpt_2,"dpt2.csv")
d1 = read.csv("dpt1.csv")
d2 = read.csv("dpt2.csv")
download.file(loc_com_1,"com1.csv")
download.file(loc_com_2,"com2.csv")
c1 = read.csv("com1.csv")
c2 = read.csv("com2.csv")

Commençons par les départements. On peut récupérer les pourcentages de voix obtenus par les différents candidats, au premier tour, et regarder les pourcentages obtenus par Marine Le Pen. Par exemple, on peut regarder, département par département, le score de Marine Le Pen au premier et au second tour, pour vérifier que le report des voix s’opère bien

D1 = d1[d1$cand_nom=="LE PEN",c("dep_code","cand_rapport_inscrits","cand_rapport_exprim")]
D2 = d2[,c("CodeDépartement","LE.PEN.ins","LE.PEN.exp")]
names(D1) = c("dep_code","ins1","exp1")
names(D2) = c("dep_code","ins2","exp2")
D = merge(D1,D2)
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

Comme on pouvait s’y attendre, on a une belle droite de régression, avec les nombres de voix obtenus, relativement aux inscrits, en bleu, et aux exprimés, à droite. Pour Emmanuel Macron au premier et Marine Le Pen au second tour, on a une belle corrélation négative, mais moins prononcée qu’auparavant. Autrement dit, un département qui votait peu Macron au premier tour a beaucoup voté Le Pen au second

Mais là, on parlait des deux candidats qui se sont maintenus. Si on regarde l’autre candidat de l’extrême droite, les reports sont assez étranges… après, les taux étaient “relativement” faibles au premier tour pour Erik Zemmour.

Pour le candidat en tête à gauche, on a une corrélation négative, mais sur les votes exprimées, on observera des outliers en haut à droite, avec manifestement un vote d’opposition au président sortant dans certains départements

Pour le candidat écologiste, la corrélation est elle encore plus clair

autrement dit, quand on a voté écologiste au premier tour, on avait peu tendance à voter Le Pen (ou disons les départements où on votait fortement écologiste n’ont pas voté Le Pen au second #ecologicalfallacy, etc).

Mais on peut faire le même exercice par commune. Mon seul soucis a été que les codes communes dans les deux bases étaient différents, donc j’ai utilisé les noms des communes, ce que je n’aime pas faire (avec les accents, les tirets, etc, la fusion peut être assez mauvaises… mais là ça marche bien). Je me suis ensuite limité aux grosses communes, de plus de 10000 inscrits,

C1 = c1[c1$cand_nom=="LE PEN",c("commune_name","cand_rapport_inscrits","cand_rapport_exprim","votants_nb")]
C2 = c2[,c("Commune","LE.PEN.ins","LE.PEN.exp")]
names(C1) = c("com_code","ins1","exp1","taille")
names(C2) = c("com_code","ins2","exp2")
D = merge(C1,C2)
D = D[D$taille>5000,]
par(mfrow=c(1,2))
plot(D$ins1,D$ins2,xlab="inscrits (%), premier tour",ylab="inscrits (%), second tour",col=colr2[1])
plot(D$exp1,D$exp2,xlab="exprimés (%), premier tour",ylab="exprimés (%), second tour",col=colr2[2],main="Le Pen / Le Pen",cex.main=1.1)

On peut alors faire le même exercice, avec Le Pen au premier et second tour

ou Macron au premier, Le Pen au second

On observe des beaux nuages, sur la gauche, respectivement en haut et en bas… Beaucoup de communes qui avaient peu voté pour Le Pen semblent avoir eu peur, et on massivement voté au second pour elle. Alors que ce n’est pas trop le cas pour celles où Le Pen avait obtenu beaucoup de voix. Et pour Macron, c’est en bas à gauche, autrement dit, beaucoup de communes qui n’ont pas voté Macron au premier tour n’ont pas pour autant voté Le Pen au second. Pour ceux que ça amuse, les codes sont disponibles pour regarder qui sont ces communes dans les coins.

Pour Zemmour, c’est assez confus, j’avoue, le personnage est déroutant (ou dégoutant ? comme le suggérait le correcteur automatique)

ou alors Jean Luc Mélanchon

ou la relativement franche corrélation des écologistes, qui semblent résolument aller en sens inverse avec le vote d’extrême droite…

 

Du “Fichier des Victimes Indemnisées” à DataJust (et inversement)

Alors que je préparais un exposé sur la décision d’aller ou pas au tribunal, dans le cas d’accident de la route (pour une victime), j’ai voulu faire une digression autour de DataJust, en lien avec ce que j’ai pu voir en travaillant sur notre article, sur la base des données du Fichier des Victimes Indemnisées, mis en place dans le cadre de la loi Badinter, en 1985, car ça permet d’avoir une perspective intéressante, il me semble.

En commençant par la fin, DataJust, née du le décret nº 2020-
356 du 27 mars 2020 portant création d’un traitement automatisé de données à caractère personnel a été publié au Journal officiel du 29 mars 2020. L’article 1er du décret du 27 mars 2020 énonce que le ministre de la Justice « est autorisé à mettre en œuvre, pour
une durée de deux ans, un traitement automatisé de données à caractère personnel, dénommé « DataJust », ayant pour finalité le
développement d’un algorithme ». Un « algorithme »… le (gros) mot est lâché.  Car autant être clair: comme le rappellent Eugénie Petitprez, et Rodolphe Bigot (dans standard humain ou standardisation
algorithmique de l’évaluation du dommage corporel ?) « la Justice n’a pas la culture de la donnée, ni les réels moyens des ambitions qui lui sont prêtées. S’il n’y a pas de culture de la donnée, c’est parce que depuis toujours la Justice fonctionne fondamentalement à l’encontre de la donnée ». On ne pourrait être difficilement plus clair. Mais les deux premiers articles donnaient des informations intéressantes. Le premier article revient sur les finalités de cet « algorithme » qui sera créé : il doit « servir à : 1 º La réalisation d’évaluations rétrospectives et prospectives des politiques publiques en matière de responsabilité civile ou administrative ; 2º L’élaboration d’un référentiel indicatif d’indemnisation des préjudices corporels (…) ». Mais cet « algorithme » n’est pas tout, car un « algorithme » n’est rien sans données pour apprendre (dans la version moderne de l’apprentissage). Et le second article évoque les données utilisées. Plus précisément, dans les données, figureront « 2º les éléments d’identification des personnes physiques suivants : la date de naissance, le genre, le lien de parenté avec les victimes et le lieu de résidence ; 3º Les données et informations relatives aux préjudices subis, notamment : – la nature et l’ampleur des atteintes à l’intégrité, à la dignité et à l’intimité subies, en particulier la description et la localisation des lésions, les durées d’hospitalisation, les préjudices d’agrément, esthétique, d’établissement, d’impréparation ou sexuel, les souffrances physiques et morales endurées, le déficit fonctionnel, ainsi que le préjudice d’accompagnement et d’affection des proches de la victime directe ; – les différents types de dépenses de santé (notamment frais médicaux, paramédicaux, pharmaceutiques et d’hospitalisation) et d’aménagement (notamment frais de logement, d’équipement et de véhicule adaptés) ; – le coût et la durée d’intervention des personnes amenées à remplacer ou suppléer les victimes dans leurs activités professionnelles ou parentales durant leur période d’incapacité ; – les types et l’ampleur des besoins de la victime en assistance par tierce personne ; – les préjudices scolaires, universitaires ou de formation subis par la victime directe ; – l’état antérieur de la victime, ses prédispositions pathologiques et autres antécédents médicaux ; 4 º Les données relatives à la vie professionnelle et à la situation financière, notamment la profession, le statut, les perspectives d’évolution et droits à la retraite, le montant des gains et pertes de gains professionnels des victimes, ainsi que des responsables ou personnes tenues à réparation ; 5º Les avis des médecins et experts ayant examiné la victime et le montant de leurs honoraires ; 6º Les données relatives à des infractions et condamnations pénales ; 7º Les données relatives à des fautes civiles ; 8º Le numéro des décisions de justice. (….) ». Je vais m’arrêter là mais il y a beaucoup d’information.

En lisant un peu des articles de presse, de juristes (on pourra lire Le raisonnement par algorithme : le fantasme du juge-robot, par exemple de Frédéric Rouvière, ou Le droit isométrique : un nouveau paradigme juridique né de la justice prédictive, de Louis Larret-Chahine) il y a beaucoup de fantasmes autour de de l’intelligence artificielle, des « algorithmes », etc. Ce qui me frappe le plus, c’est que dans un cas très particulier ces données existent (plus ou moins) : c’est le cas des dommages corporels subits par les victimes d’accident de la circulation lors d’accidents automobiles

Il existait auparavant des banques de données des décisions des différentes juridictions, mais rien sur les règlements contractuels, c’est à dire, pour faire simple, lorsque l’assureur offre une indemnité à la victime (et que cette dernière l’accepte). Pour rendre comparables les montants, la loi (qui date de plus de 35 ans) a créé le “fichier national des personnes indemnisées“, où doivent figurer les montants d’indemnités résultant tant des décisions de justices que des transactions avec les assureurs. Cette obligation d’information a été inscrite dans la lettre de la loi. L’article 26 de la loi du 5 juillet 1985 (c’est la petite phrase de l’article L.211-23 du Code des Assurances), « sous le contrôle de l’autorité publique, une publication périodique rend compte des indemnités fixées par les jugements et les transactions ». Et depuis 1989, cette base est accessible (à l’époque, par…. Minitel – 3615 AGRIA). Comme le montre l’extrait que j’indique, il ne s’agit aucunement de garantir un niveau d’indemnisation à la victime, mais seulement de lui assurer la disponibilité de l’information, afin qu’elle puisse juger par elle-même de la valeur des règlements qui lui sont proposés. Je renvoie vers l’accident corporel de la circulation, entre transactionnel et juridictionnel de Evelyne Serverin, Boris Bertin, Simy Elbaz et Sophie Julliot pour un survol passionnant (le rapport date de 1997) !

Dans l’exposé que je donnerai lundi midi, on reviendra sur l’utilisation de ces données, car elles sont riches. Mais complexes à analyser. Un point intéressant, c’est qu’il y a un enregistrement par victime (et non pas par accident). Et pour chaque victime, on a le montant de l’indemnité perçue. Pour un règlement contractuel, c’est le montant offert par l’assureur, et pour un jugement, c’est le montant ordonné par le juge (le montant proposé par l’assureur est alors inconnu). Pour un économètre, c’est un joli problème, quand on essaye de comprendre qui décide d’aller au tribunal. Mais les conclusions de notre étude ont surtout souligné l’incroyable variabilité des indemnités perçues par les victimes. Y compris ceteris paribus sic stantibus. Et compte tenu du nombre d’information à disposition, cela signifie que l’on peut regarder à âge de la victime donné, à genre donné, à lieu de résidence donné, à gravité de l’accident donné (niveau d’IPP, nombre de jours d’hospitalisation, etc). Si cet « algorithme » est simplement un barème, une grille indicative d’indemnisation, j’imagine que l’idée est louable. Si notre modélisation donne des niveaux obtenus au tribunal et après négociation avec l’assureur très proches, la variabilité quand on compare aux montants réellement observés est colossale ! Dans notre étude, on observe une personne qui aurait eu un accident en 1996, dans la région d’Aix-en-Provence, passager, d’un peu plus de 50 ans, qui a passé un peu moins de 2 ans à l’hôpital. Alors que nos modèles prédisent un dommage corporel de l’ordre de 65,000 euros, elle a touché presque 100,000 euros (au tribunal). Un homme dans la région de Rennes, un peu moins de 50 ans, qui a passé moins de 2 ans à l’hôpital aurait du toucher selon nos modèles de l’ordre de 130,000 euros, et il a accepté les 100,000 euros proposés par l’assureur. Et ce ne sont que deux exemples pris au hasard.

La variabilité observé pour des victimes « semblables » (le fameux ceteris paribus sic stantibus) rend mal à l’aise, et laisse à penser qu’un barême national serait bénéfique à tous et à toutes (je ne parlerais pas ici des différences entre les hommes et les femmes que nous observons, là encore des victimes « semblables » … ou de la différence entre les lieux des accidents, en particulier entre Aix-Marseille et Rennes). Avoir des données pour mieux comprendre comment sont indemnisées les victimes d’accidents corporels est essentiel pour avoir un minimum de justice. Et DataJust n’apporte finalement rien de plus que ce que proposait la loi Badinter (tant que l’« algorithme » reste un modèle de régression – même si ici, on est dans cas de données dépendantes limitées pour reprendre la terminologie de Madalla (1983), qui restent des modèles classiques d’économétrie, mais qui pourraient ressembler à de la magie noire de l’intelligence artificielle à un juriste éloigné des techniques statistiques de base).

Bref, j’étais étonné de voir autant de plumes s’élever contre DataJust, surtout quand on voit que des choses très proches existent depuis 35 ans, avec des finalités louables

En l’occurrence, notre étude montre que les cours donnent des indemnités très différentes, et c’est problématique. La culture de la donnée, et des modèles, est indispensable si on veut un minimum de transparence et d’équité. C’est qu’évoque le retrait de DataJust, ou la fausse défaite des barèmes de Vincent Rivollier, Manon Viglino et Christophe Quézel-Ambrunaz.

Mais cela suppose d’avoir des historiques de données importants (dans notre cas, avec seulement 20 ans de données, on a du mal à avoir des estimations robustes, en grande partie à cause de la très grande variabilité des montants), et des données fiables, et bien renseignées. Car ensuite l’interprétation devient complexe: lorsque nous avions présenté nos résultats il y a quelques années, on nous avait fait remarqué que si les montants ne sont pas « semblables » alors que les victimes le sont, c’est peut-être qu’elles ne sont pas aussi « semblables » qu’elles semblent l’être. En particulier, les informations relativement subjectives (que je dirais sorties d’un « algorithme » si je voulais être taquin) comme le niveau d’I.P.P. – Incapacité Permanente Partielle – n’est peut être pas attribué de la même manière à tout le monde. Au lieu de monter sur ses grands chevaux quand on parle d’ « algorithme », il serait intéressant d’écouter les personnes qui font de l’économie de la justice, ou de statistique de la justice. Et je suis content, parce que je vais les rencontrer cette semaine lors de mon passage à Paris ! À suivre donc…

Could there be incentives to cycle through a red light?

This is of course a rhetorical question! Because cyclists must stop when the light is red! … But … there is always that moment, on a bicycle, when you stop, and  then you say to yourself

the worst part is that the lights are badly regulated, and I know that the next one will also be red once I will reach it … whereas if I had passed, I would have had the next one, and who knows, maybe a green wave afterwards?

Not having wanted to try the experiment with my bike, I wanted to try using simulations, on my computer.

Let us assume that in my city, there are red lights every 250 m, and they go from red to green in 1 minute, and similarly from green to red. Then, I can use a simple loop, where I compute the time from home, each time I pass a light: either the light is green, I go through, and I do not wait; or the light is red, then I wait, until the end of the minute. Here is the code. First, I define the constant parameter, with lights every 250 m, light turning either red or green every 1 minute, and I assume a constant speed here, of 15/60 km per min (or 15km per hour).

v = 15/60
d = .250
t = 1

Suppose that the office is at 10 km from home. The code to compute the time is simply (I do add a random noise on the time it takes to reach the next light, I will get back to that later on)

Dist = seq(0,10,by=.25)
Time1=rep(NA,length(D))
Time1[1]==0
for(k in 2:length(Dist)){
  noise=rnorm(1,sd=.05)
Time1[k]=Time1[k-1]+d/v+noise
if((floor(Time1[k])%%2==0)&(k<length(Dist))) Time1[k]=ceiling(Time1[k])
}

I can now visualize myself on my bike

T=seq(0,60,by=1/60)
colr = c("red","green")
plot(T,T*v,col="white")
for(k in 1:40){
points(T,rep(d*k,length(T)),cex=.3,pch=15,col=colr[rep(1:2,each=t*60)])
}
points(Time1,Dist,pch=19,cex=.5)

Here, it took me about 51 min to reach the office. And each time I have a red light, I stop and wait. Here is the distribution of the time it will take, as a function of my speed

simul = function(v=.250,sd=.05){
  Dist = seq(0,10,by=.25)
  Time1=Time2=rep(NA,length(D))
  Time1[1]=Time2[1]=0
  for(k in 2:length(Dist)){
    noise=rnorm(1,sd=sd)
    Time1[k]=Time1[k-1]+d/v+noise
    if((floor(Time1[k])%%2==0)&(k<length(Dist))) Time1[k]=ceiling(Time1[k])
  }
  max(Time1)
}
S=function(v=.250,sd=.05,n=1000){
  Vectorize(function(x) simul(v=x,sd=sd))(rep(v,n))
}
vit = seq(.2,.3,length=101)
MS = Vectorize(function(x) S(v=x,sd=.05))(vit)
qMSsup = apply(MS,2,function(x) quantile(x,.95))
qMSinf = apply(MS,2,function(x) quantile(x,.05))
MSmed = apply(MS,2,function(x) quantile(x,.5))
MSmean = apply(MS,2,mean)
par(mfrow = c(1,1))
plot(vit,MSmed,type="l")
polygon(c(vit,rev(vit)),c(qMSinf,rev(qMSsup)),col="light blue",border=NA)
lines(vit,MSmed,lwd=2)

Obviously, on average, the faster I cycle, the shorter the ride will be. Of course, 40 min is not a (real) lower bound : if I go much faster, it can be a rather quick ride

Now, consider a simple rhetorical alternative. What if I decide to go through a red light (after checking that there is no danger), say, with 1 chance of out 20 ?

simul_random = function(v=.250,sd=.05){
  Dist = seq(0,10,by=.25)
  Time1=rep(NA,length(D))
  Time1[1]=0
  for(k in 2:length(Dist)){
    noise=rnorm(1,sd=sd)
    Time1[k]=Time1[k-1]+d/v+noise
    red = sample(c(0,1),size=1,prob = c(.95,.05))
    if((floor(Time1[k])%%2==0)&(k<length(Dist))&(red == 0)) Time1[k]=ceiling(Time1[k])
  }
  max(Time1)
}

Here, 1 chance out of 20 could mean that on a 10 km ride, I wil always stop : first, I need a red light, and then 95% of the time, I stop. Here is the average distribution, with confidence bands

S_random=function(v=.250,sd=.05,n=1000){
  Vectorize(function(x) simul_random(v=x,sd=sd))(rep(v,n))
}
MS3 = Vectorize(function(x) S_random(v=x,sd=.05))(vit)
qMSsup3 = apply(MS3,2,function(x) quantile(x,.95))
qMSinf3 = apply(MS3,2,function(x) quantile(x,.05))
MSmed3 = apply(MS3,2,function(x) quantile(x,.5))

It is possible to compare the two actually (on average times)

plot(MSmed,MSmed3,type="l")
abline(a=0,b=1,col="red",lty=2)
plot(MSmed,(MSmed-MSmed3)/MSmed*100,type="l")

which means that, somehow I can save some time, maybe from 3% to 5% if I do not cycle to fast, otherwise probably less than 2%. I would not claim, here, that it is worth it.

Consider another rhetorical alternative. What if I decide to go through the red light only if it is during the very first second ?

simul_1sec = function(v=.250,sd=.05){
  Dist = seq(0,10,by=.25)
  Time1=rep(NA,length(D))
  Time1[1]=0
  for(k in 2:length(Dist)){
    noise=rnorm(1,sd=sd)
    Time1[k]=Time1[k-1]+d/v+noise
    if((floor(Time1[k])%%2==0)&(k<length(Dist))&(Time1[k]-floor(Time1[k])>1/60)) Time1[k]=ceiling(Time1[k])
  }
  max(Time1)
}

Here is the time it takes to go to the office

and here is the potential gain

We can now clearly see the impact of those nonlinearities : with my average 15 km per hour speed, I can save 8% of the time if I go through red during the very first second. Of course, I will never do such a think, but mathematically speaking, it is stricking. Here is the gain for the first 2 seconds (out of a full minute)

so the gain can be about 12%. And now, what if we remove the noise – or more precisely, what if the standard deviation become 0.001 instead of 0.05. Here is the first distribution of the time (when I stop at each red light, and wait)

If I do not go fast enough, I will stop at almost each light, and wait until the end of the minute: going at 20 km per hour or 24 km per hour is exactly the same. But if I can go slightly faster than 25km per hour, that is awesome, and I have my green wave. Now, her is the graph I get when I cycle through the red light only in the very first second

This is now the gain

I can save up to 40% of the time, and more realistically, my 55 min ride could now be a 40 min ride. Which is substantial.

But on that one, observe that another strategy is also possible : what if I do not cycle through red light, but I simply cycle +1% faster ?

simul_faster = function(v=.250,sd=.05){
  Dist = seq(0,10,by=.25)
  Time1=rep(NA,length(D))
  Time1[1]=0
  for(k in 2:length(Dist)){
    noise=rnorm(1,sd=sd)
    Time1[k]=Time1[k-1]+d/(v*1.01)+noise
    if((floor(Time1[k])%%2==0)&(k<length(Dist))&(Time1[k]-floor(Time1[k])>1/60)) Time1[k]=ceiling(Time1[k])
  }
  max(Time1)
}

We have a similar result here, since riding 1% faster can actually help me save up to 50% of my time ! and (more realistically) if I ride 1% faster with a large random noise in the time between two light, I get almost the same as previously (when passing through the red light)

Again, my point here is not that we should cycle through a red light, of course not. But there may be accumulation of (small) nonlinear effects that might have a major impact at the end. And I believe that the best way to avoid this is to offer a green wave to cyclists, assuming that they ride at a reasonnable speed…

Un double centenaire : Treatise on probabilities de John Maynard Keynes et Risk, Uncertainty and Profit de Frank Knight

Il y a tout juste cent ans, paraissait deux ouvrages fondamentaux essayant de lier l’économie, le risque, l’incertain et le calcul des probabilités. Et si l’occasion est rare de se plonger dans des textes aussi anciens, y compris pour la majorité des étudiantes et des étudiants, il est agréable de noter que ces ouvrages contiennent de nombreuses réflexions passionnantes, qui éclairent d’un œil (presque) nouveau ces domaines où nous avons parfois l’impression de tout savoir, depuis tout ce temps…

Cambridge, Keynes, et les probabilités

Commençons par Treatise on probabilities, le premier ouvrage publié par un des plus grands économistes du début du XXème siècle, et qui témoigne de l’importance des mathématiques dans le formalisme économique – John Maynard Keynes avait déjà publié des articles économiques depuis une dizaine d’année lorsque le traité est publié. Et étonnamment, cet ouvrage est souvent oublié par les économistes[i].

Mais commençons par le début… John Maynard Keynes est né à Cambridge en juin 1883, son père (John Neville Keynes) était professeur de logique et d’économie à l’université de Cambridge, et sa mère (Florence Ada Brown) deviendra maire de la ville (en 1932, une autre femme – Eva Hartree – avait déjà été élue maire dix ans auparavant[ii]). Au début du XXème siècle (en réalité entre 1880 et 1940, comme le montre MacLeod & Urquiola (2020)), Cambridge était probablement la plus prestigieuse université au monde. En 1903, Bertrand Russell y publie son Principles of Mathematics (dont le titre de la traduction en français – Écrits de logique philosophique – reflète peut-être mieux le contenu que le titre original). L’ouvrage de John Maynard Keynes est d’ailleurs dans cette lignée, car malgré un formalisme mathématique indéniable, le traité est peut-être avant tout un ouvrage de philosophie et de logique. C’est aussi l’époque où Srinivasa Ramanujan est sur le campus, à Cambridge, invité du Tamil Nadu, en Inde, pour travailler avec Godfrey Harold Hardy et John Littlewood, comme le raconte The Man Who Knew Infinity, de Robert Kanigel, qui décrit le Cambridge des années 1910. John Maynard Keynes y côtoiera aussi de nombreux artistes (en tant que membre du groupe, ou association, Bloomsbury) comme Virginia Woolf, ou E. M. Forster

En 1905, alors qu’il vient d’obtenir un Bachelor of Arts en mathématiques, il suit les cours d’Alfred Marshall, qui vont lui faire découvrir l’économie politique, et en 1908 (de retour de deux ans passés à Bombay au Bureau des Affaires Indiennes), il travaille sa thèse sur les probabilités. Il ne s’agit pas formellement d’une thèse de doctorat[iii], mais d’une « Fellowship Dissertation » afin de devenir Fellow du King’s College. Son manuscrit est presque fini lorsque la première guerre mondiale éclate, ce qui empêchera sa publication, et il retravaillera le document pendant encore 2 ans après la fin de la guerre, pour la voir publiée en 1921, sous le titre Treatise of Probabilities. Coïncidence, 1921 est aussi l’année ou Franck Ramsey arrive à Cambridge (il deviendra également Fellow de King’s College en 1924 avant de décéder en 1930, à tout juste 26 ans), dont la proximité intellectuelle avec John Maynard Keynes est importante, même s’ils auront beaucoup de divergence, en particulier sur le concept de probabilité. 1921, c’est enfin un peu moins de 10 ans avant la publication de l’axiomatique des probabilités par Kolmogorov, en 1929, qui placera les probabilités au sein de la théorie de la mesure, et qui proposera le cadre formel encore enseigné aujourd’hui.

Pour comprendre l’orientation prise par Keynes dans son livre, il faut avoir en mémoire les liens entre probabilités et logique. A cette époque, à Cambridge, William Ernest Johnson travaille sur son ouvrage Logic, dont les trois tomes[iv] paraitront en 1921, 1922 et 1924. Harold Jeffreys, à qui on devra de nombreux travaux en statistique Bayésienne, était aussi à Cambridge à cette époque, ainsi que Dorothy Wrinch, ce dernier ayant publié On Some Aspects of the Theory of Probability en 1919. Si la recherche était bouillonnante à cette époque, R. B. Braithwaite, qui signe la réédition de l’ouvrage de Keynes en 1972, note que malgré tout, rares étaient les ouvrages de référence en anglais sur le sujet, le précédant ayant été publié plus de 50 ans auparavant par John Venn (qui aura laissé à la postérité les fameux diagrammes) avec son Logic of Chance. Et si Harold Jeffreys a proposé une axiomatique des probabilités plus solide que le modeste chapitre 2 du Traité, son ouvrage ne sortira qu’en 1939.

Pour Keynes, une « probabilité » n’est pas un fait objectif, naturel, comme le laisse entendre l’interprétation fréquentiste.  Pour Keynes, une probabilité est avant tout affaire de logique, quantifiant un degré de confirmation, ou de croyance rationnelle, d’une proposition. Formellement,

« Let our premises consist of any set of propositions h, and let our conclusion consist of any set of propositions a, then, if a knowledge of h justifies a rational belief in a of degree a,we say that there is a probability-relation of degree a between a and h. »

Elle exprime le degré de croyance qu’il est raisonnable d’avoir envers un ensemble de propositions a, à la lumière d’un ensemble de propositions h. Dans cette mesure, donc, la probabilité peut être qualifiée de subjective, car il s’agit d’une croyance.  Mais dans le sens que lui donne Keynes, la probabilité n’est pas subjective puisqu’elle n’est pas soumise à un caprice humain : c’est une relation logique entre deux propositions (ou deux ensembles de propositions) qui se manifestent dans l’esprit des individus. Ce n’est pas une relation entre une déclaration et une réalité, « probability begins and ends with probability ».  Dans le livre de Keynes, une probabilité peut être exprimée par le symbole ‘a/h‘. La connaissance d’une situation dans laquelle ‘a/h = 1′ est certaine. Une situation dans laquelle ‘a/h = 0′ constitue une impossibilité. Dans la majorité des cas, ce chiffre est situé entre 0 et 1. Et, dans la majorité des cas, les nombres n’ont qu’une signification ordinale. Keynes veut dire par là qu’il est souvent impossible de comparer quantitativement deux probabilités et d’affirmer, par exemple, que les chances qu’un événement se produise sont x fois plus élevées que les chances qu’un autre événement se produise.  De plus, il est possible que deux probabilités soient incomparables : « it is not  always possible to say that the degree of our rational belief in one conclusion is either equal to, greater than, or less than the degree of our belief in another ». Et l’on peut penser que, dans la réalité humaine, les probabilités sont dans la plupart des cas non quantifiables et incommensurables, d’autant plus qu’elles sont contraintes par les limites de la raison humaine. Même lorsque les probabilités individuelles sont mesurables numériquement, nous ne pouvons pas aller très loin dans le raisonnement mathématique. Il est d’ailleurs difficile d’éliminer l’intuition et le jugement direct de la considération des probabilités. Cette difficulté restreint la portée des probabilités fréquentistes fondées sur la loi des grands nombres. Même dans le domaine des sciences naturelles, l’intuition et l’analogie jouent un rôle plus important que la manipulation des fréquences statistiques.

En fait, pour Keynes, la plupart des grands penseurs des probabilité, comme Condorcet, Bernoulli, Bentham, Laplace, ou Edgeworth, se sont trompés en pensant que l’on pouvait appliquer les principes issus de l’équiprobabilité aux sciences morales et ainsi quantifier, mesurer et formaliser mathématiquement la réalité sociale.  Keynes rappelle que certains de ces grands penseurs se sont aventuré dans le domaine de l’éthique, et ont a été amené à croire que les degrés de bonté étaient numériquement mesurables, et parfois additifs, alors que la question relève du jugement intuitif.  Keynes dénonce cette « mathematical charlatanry » qui a fait perdre du temps à la statistique théorique, et qui est soutenue par une confiance totale dans l’inférence statistique.  En économie, cette erreur est soutenue par l’illusion de la quantification selon laquelle des choses comme l’utilité peuvent être quantifiées, mesurées, additionnées et soustraites. Cette erreur conduit à une utilisation abusive des statistiques qui, en tant que moyen de décrire quantitativement la réalité, s’est transformée en un instrument de prédiction. Keynes poursuivra cette réflexion par la suite, avec en particulier des éléments de réflexions sur les probabilités dans sa Théorie Générale, avec le Chapitre 12. En particulier, c’est dans ce second ouvrage que l’on trouve la note de bas de page abondamment citée, où il affirme que « très incertain » et « très improbable » sont pour lui équivalent, laissant entendre que les probabilités serraient l’outils approprié pour modéliser les situations incertaines (récemment, Brady (2019) dédie un article à cette fameuse note 1 de la page 148). Pour une discussion plus poussée sur le contenu du Traité, je vais renvoyer à l’article de Christian Robert sur le sujet.

Pour revenir sur l’histoire du Traité, comme nous l’indiquions en introduction, il est passé inaperçu pour bon nombre d’économistes, tout du moins jusque dans les 1970, cinquante ans après sa parution. Lors de la réédition des travaux de Keynes, il a été décidé de publier ce livre non pas dans l’ordre chronologique comme cela avait été fait pour les sept autres monographies, mais comme un huitième titre. Au début des années 1980, les travaux de Robert Skidelsky, suivis par ceux de Bradley Bateman, Anna Carabelli ou Michael O’Donnell, vont inaugurer un ensemble de travaux étudiant les relations entre philosophie et économie chez Keynes ainsi que la continuité de sa vision philosophique. On peut noter cette réflexion de Peter Waley

« The first major effort to construct a theory of imprecise probability was made by Keynes (1921). Keynes aimed to develop an inductive logic, based on a logical interpretation of probability as a “degree of rational belief” […] Since Keynes, a large literature has grown up concerning the mathematics and interpretation of imprecise, epistemic probabilities »

Mais il ne serait pas honnête d’oublier la critique que Franck Ramsey a faite dans Truth and Probability (publié à titre posthume en 1931, mais basée sur une conférence donnée en 1926).  Pour ce dernier, Le calcul des probabilités consiste à établir un ensemble de règles permettant aux degrés de croyance de former un système cohérent. Les relations de probabilité, telles que Keynes les avait décrites, n’existent pas. Ce questionnement sur les probabilités, l’incertain et les croyances vont se retrouver dans un second ouvrage, publié lui aussi en 1921, mais dans un contexte assez différent.

Chicago, Knight et l’incertain

Le second ouvrage dont on célèbre le centenaire de la publication, Risk, Uncertainty and Profit de Frank Knight a été publié de l’autre côté de l’Atlantique, indépendamment de ce qui se passait en Angleterre. En 1871, Chicago est dévasté par un immense incendie, et sa reconstruction se fera grâce à l’école de Chicago, qui attirera les plus grands architectes, entre 1890 et 1910, avec William Le Baron Jenney, Henry Hobson Richardson, ou Frank Lloyd Wright qui est à Chicago pendant une bonne partie de cette période. Chicago atteint le million d’habitants et devient un moment une ville importante artistiquement, même si on retient souvent qu’elle était la capitale du crime organisé, avec Al Capone par exemple. Et si notre histoire ne commence pas à Chicago, c’est dans cette ville de Franck Knight s’imposera, en dirigeant, avec Jacob Viner, le département d’économie, entre 1920 et 1950. Même si c’est Milton Friedman qui rendit célèbre l’école de Chicago, Franck Knight a marqué par ses enseignement la génération rendue célèbre après-guerre, comme le rappelle Van Overtveldt (2007).

Dans les années 1890, et pendant une quinzaine d’années, les économistes américains essayent de répondre à la question de la justice sociale du partage et de la répartition des richesses. On peut retenir l’article de Frederick Hawley, The Risk Theory of Profits parmi les nombreux articles parus dans le Quarterly Journal of Economics, édité à Harvard, où les économistes se répondent les uns aux autres, pendant plusieurs années. Pour Frederick Hawley, la prise de risque était une composante inévitable de la dynamique de la production, et ceux qui prenaient des risques devaient avoir droit à une récompense, connue sous le nom de « profit ». Selon Frederick Hawley, le profit est le prix payé par la société pour assumer le risque. En particulier, il note que cette prime de risque devrait correspondre à une compensation supérieure à la valeur actuarielle, c’est-à-dire une prime sur le risque calculable. Autrement dit, le profit naît du risque non assurable, et voyant le profit comme un résidu, ce qui a donné le nom « residual theory of profit ». S’en suivra une longue polémique avec John Bates Clark, jusqu’à l’article de Francis Edgeworth en 1904, qui offrira une synthèse du sujet.

Maintenant que le contexte est posé, revenons un instant à Frank Knight. Frank Knight est né dans l’Illinois, en novembre 1885. Il était le fils aîné d’une famille de fermiers, et les contraintes de la vie à la ferme l’ont forcé à mettre ses études de côté un certain temps. Il a quitté sa famille alors qu’il avait plus de 20 ans pour commencer des études collégiales (évangéliques) dans le Tennessee, pour obtenir un Bachelor of Science en sciences naturelles, puis une maîtrise en allemand quelques années plus tard. Il part ensuite à Cornell, pour un doctorat en économie (discipline qu’il découvre à partir de 1913), et défendra sa thèse, intitulée « A Theory of Business Profit », qui avait remporté la deuxième place dans la catégorie générale du concours d’essais économiques Hart, Schaffner & Marx en 1916. A cette époque à Cornell, on retrouve Allyn Abbott Young qui y restera avant d’enseigner à Harvard, puis à la London School of Economics, qui présidera le jury de doctorat de Frank Knight, et qui lui offrira de rester comme stagiaire post-doctoral, jusqu’à son départ pour Harvard en 1917. Frank Knight de son côté rejoindra alors John Bates Clark à l’Université de Chicago, où il retravaillera le manuscrit de sa thèse, qui sera publiée en 1921 sous le titre Risk, uncertainty and profit, comme le raconte Emmet (2020) – qui revient longuement sur les changements entre le manuscrit de la thèse, et l’ouvrage publié. Les superviseurs de Knight – John Bates Clark et Allyn Abbott Young – ont attesté à l’éditeur que peu de choses avaient été modifiées, mais comme pour John Maynard Keynes, ces années de maturation ont permis de clarifier certains concepts.

Le point de départ de livre de Frank Knight, ce sont les travaux de John Clark, et en particulier Distribution of Wealth paru en 1899,

« Professor Clark’s work is taken as the basis and starting point of the present study. Some criticism of it will be found necessary, leading to disagreement on some rather essential points, but with its general aim and with many of its conclusions our own are in accord ».

Tout au long du livre, il revient sur des défauts des approches de Hawley et Clark, en notant par exemple que si les agents économiques connaissaient les distributions de probabilités, ils devraient alors se comporter exactement de la même manière. Même si Clark avait introduit une notion de « valeur actuarielle subjective », l’idée centrale à l’époque était que les probabilités (et donc la valeur actuarielle) était parfaitement connue par tous les agents

« It is this “true” uncertainty, and not risk, which forms the basis of a valid theory of profit and accounts for the divergence be-tween actual and theoretical competition ».

La première idée de Frank Knight est bien sûr la notion d’incertitude. L’ouvrage récent de John Kay et Mervyn King, Radical Uncertainty : Decision-Making Beyond the Numbers, revient sur cette idée centrale dans le livre de Knight. La deuxième idée est la théorie des prix, et plus précisément la théorie de la concurrence parfaite. La différence la plus importante, selon Knight, entre la théorie classique des prix et la concurrence imparfaite est la présence dans cette dernière de l’incertitude, qui émerge d’un changement dynamique, et donc imprévisible. C’est cette imprévisibilité dynamique qui a créé les conditions économiques pour l’action entrepreneuriale. Cette approche nécessitaient un jugement, une faculté humaine qui n’était pas nécessaire dans la concurrence parfaite, où tout était calculable, et calculé, parce que tout était « prévisible ». L’incertitude exigeait un jugement, non seulement sur les états du monde, mais aussi sur la capacité des autres à porter des jugements dans des situations incertaines. L’organisation et la coopération est donc apparue, tant au sein des entreprises que dans l’organisation économique sociale, comme un moyen de coordonner les jugements. Franck Knight a ainsi été le premier fait une distinction entre les problèmes de décision lors que les probabilités sont connues, et lorsqu’elles ne le sont pas, ce qui impose de clarifier (comme avait tenté de le faire Keynes) la notion de probabilité.

Dans Risk, Uncertainty and Profit, Knight a distingué trois types de probabilité différents, qu’il a appelées « probabilité a priori », « probabilité statistique » et « estimations ». Le premier type se situe sur le même plan logique que les propositions mathématiques, l’exemple classique étant la probabilité d’obtenir une face lors d’un lancer de dé (la symétrie du cube permet d’affirmer qu’il y a une chance sur six d’obtenir une des faces). La « probabilité statistique » dépend de l’évaluation empirique de la fréquence, dans des expériences répétées. Dans le langage bayésien contemporain, dans le premier cas, on retrouve la distribution de probabilité a priori. alors que dans le deuxième cas, on retrouverait une distribution a posteriori. Dans le troisième cas, les données existantes ne se prêtent pas à l’analyse statistique, lorsque des expériences répétées sont impossibles. Ce dernier cas est celui qui a le plus intéressé Knight, en tant qu’économiste essayant de comprendre et modéliser le monde des affaires, et la nature du profit dans ce contexte. Knight a identifié la confusion entre « le problème de l’estimation intuitive » et « la logique de la probabilité », qu’elle soit a priori ou statistique

Le livre n’a pas connu un grand succès, seulement deux articles publiés avant 1930 l’ayant cité. Mais en 1933, Lionel Robbins a fait en sorte que la London School of Economics (LSE) réimprime le livre afin que ses étudiants puissent y avoir accès. Et quarante ans plus tard, l’American Economic Association a décerné à Knight la médaille Walker, décernée tous les cinq ans pour honorer l’économiste américain qui a apporté la plus grande contribution à l’économie tout au long de sa carrière.

Arrous, J. (1982). Keynes et les probabilités: Un aspect du ” fondamentalisme ” keynésien. Revue économique, 33(5), 839-861

Bouvier-Patron, P. (1996). La question de la réduction d’incertitude chez F. Knight. L’Actualité économique, 72 (4), 397–415.

Brady, Michael Emmett, How Should the Post Keynesian School Define ‘Uncertainty’? 1 The Only Correct Answer Is to Use Keynes’s Own Definition Given in Footnote 1 on Page 148 of Chapter 12 of the General Theory: Uncertainty Is an Inverse Function of the Weight of the Argument (2016). SSRN:3438090

Combenale, P. (2010). Introduction à Keynes. La Découverte.

Diebold, F.X. Doherty, N.A. & Herring, R.J. (2010). The Known, the Unknown, and the Unknowable in Financial Risk Management: Measurement and Theory Advancing Practice. Pinceton University Press.

Emmet, R.B. (2020). The Writing and Reception of Risk, Uncertainty and Profit. SSRN: 3591596.

Kent, S. (2019). ‘A degree of a different character’: 100 years of the Cambridge PhD. Cambridge University Library Special Collections Blog

Leroy, S.F. & Singell, L.D. (1987). Knight on risk and uncertainty. Journal of Political Economy, XCV, 2, p. 294-406.

MacLeod, W.B. & Urquiola, M. (2021). Why Does the United States Have the Best Research Universities? Incentives, Resources, and Virtuous Circles. Journal of Economic Perspectives, 35:1, 185-206.

Muchlinski, E. (2003). Épistémologie et probabilité chez Keynes. L’Actualité économique, 79 (1-2), 57–70.

Pradier, PC. & Serrano, D.T. (2000) Frank H. Knight le risque comme critique de l’économie politique. Revue de synthèse 121, 79–116

Robert, C.P (2010). An attempt at reading Keynes’ Treatise on Probability. ArXiv: 1003.5544

Schmidt, C. (2003). Que reste-t-il du Treatise on Probability de Keynes? L’Actualité économique, 79 (1-2), 37–55.

Van Overtveldt, J. (2007). The Chicago school : how the university of Chicago assembled the thinkers who revolutionized economics and business, Chicago, Agate

Walley, Peter ( 1991), Statistical Reasoning with Imprecise Probabilities, London : Chapman and Hall.

Wrinch, D. & Jeffreys H. (1919) On Some Aspects of the Theory of Probability, Philosophical Magazine, no 38,‎ 1919, p. 715-731

Zabell, S.L. (1992) Predicting the unpredictable. Synthese 90, 205–232

[i] Cet ouvrage n’est même pas cité dans l’introduction à Keynes de Pascal Combemale, Combenale (2010)

[ii] En France, la première femme élue à un conseil municipal était Joséphine Pencalet, à Douarnenez en mai 1925 (suite à la grève des sardinières), mais son élection fut invalidée par le Conseil d’État, car à l’époque, les femmes ne sont ni électrices ni éligibles. Il faudra attendre 1945, et le droit de vote, pour que les femmes deviennent éligibles, en Frnace.

[iii] Le programme de Doctor of Philosophy de Cambridge n’existera qu’à partir de 1919, comme le raconte Kent (2019)

[iv] Zabell (1992) revient sur les apports de William Ernest Johnson, qui a en particulier introduit un des premiers relâchements de l’hypothèse d’indépendance, en probabilités, en définissant le concept d’échangeabilité.

From multinomial regression to binary classification on some Siamese data

There are two kinds of people in the world: people who think there are two kinds of people in the world and people who don’t

(borrowed from Menand (2018)). Because things are always simpler when we face only binary choice, aren’t they? But consider here the case were multiple options are possible, and let us see if we cannot get back to simpler binary choices. Consider a collection of observations (y_i,\boldsymbol{x}_i) where y_i is some categorical variable, y_i\in\mathcal{A} where \mathcal{A}=\lbrace A_1,\cdots,A_\kappa \rbrace, with \kappa possible categories. Let \mathcal{I}_k=\lbrace i:y_i\in A_k \rbrace.

In a classical multinomial logistic regression, suppose that A_1 is the reference, then \mathbb{P}[Y=A_j|\boldsymbol{X}=\boldsymbol{x}]=\frac{\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_j]}{1+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_2]+\cdots+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]}With a lot a categories, and a small number of observations, inference can be complicated, and non-robust.

  • the Siamese dataset

The name Siamese I use, here, comes from Siamese Networks. Or sort of… As we say in French, it is an « histoire de l’homme qui a vu l’homme qui a vu l’ours » (story of the man who saw the man who saw the bear). A few years ago, a student tried to explain to me the idea of Siamese Networks and this is what I understood. I might be completely wrong, but the idea I got from it did make sense, in my mind at least. That is the story of that blog post…

The idea of the siamese algorithm will be to consider all pairs of observations, (y_i,\boldsymbol{x}_i) and (y_j,\boldsymbol{x}_j) :

  1. \tilde y_{i,j}=\boldsymbol{1}(y_i=y_j) indicating if individuals i and j are in the same category
  2. \tilde{\boldsymbol{x}}_{i,j} is a collection of p-1 variables,
  • \tilde {x}_{k:i,j}={x}_{k:i}-{x}_{k:j} if x_k is continuous, or \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}| (we can use another metric, e.g. \tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|^2, and this is why I decided to use some GAM model in the logistic regression on the Siamese dataset)
  • \tilde {x}_{k:i,j}=({x}_{k:i},{x}_{k:j})\in\mathcal{X}_k\times\mathcal{X}_k if x_k is a categorical variables (taking values in the set \mathcal{X}_k), or \tilde {x}_{k:i,j}=\boldsymbol{1}({x}_{k:i}\neq{x}_{k:j})\in\{0,1\}

The original dataset was a n\times p matrix, and (if there are no categorical variable), it becomes a n(n-1)/2\times p matrix. The key point is that if the original variable y_i was multinomial, y_{i,j} is now binomial. For instance, if our initial dataset was the following, with two covariates, one continuous and one categorical

its siamese counterpart is the following

  • Classification step

On the dataset (\tilde y_{i,j},\tilde {x}_{k:i,j})_{i,j}, fit a logistic regression, \mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}]=\frac{\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}(or any classification model – CART, random forest, etc). But that is the easy part (unless n is large, because the siamese dataset has (roughly) n^2/2 rows). The difficult task is the prediction

  • Prediction step

Consider a new input variable \boldsymbol{x}_{\cdot}, and define its siamese version, \tilde{\boldsymbol{x}}_{\cdot}=(\tilde{\boldsymbol{x}}_{\cdot,j})_j, i.e. a database with n rows. Then compute
p_{\cdot,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{\cdot,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]} where p_{\cdot,j} is the probability that (y_j,\boldsymbol{x}_{j}) and (y_{\cdot},\boldsymbol{x}_{\cdot}) are in the same category, as well as p_{i,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{i,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}
Let \boldsymbol{p}_{\cdot}=(p_{\cdot,j}), and similarly \boldsymbol{p}_{i}=(p_{i,j}). Then several techniques can be used to predict y_{\cdot}.

  1. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j=1,\cdots,n}{\text{argmax}}\{p_{\cdot,j}\}: the predicted class is the one of the observation the most likely to be other same class
  2. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{\ell=1,\cdots,k}{\text{argmax}}\{\overline{p}_{\ell}\}, where\overline{p}_{\ell} = \frac{1}{n_{\ell}}\sum_{j\in\mathcal{I}_s} \boldsymbol{1}(y_j =y_{\ell}),\text{ where }\mathcal{I}_s=\lbrace i:p_{\cdot,i}>s\rbraceconsider only probabilities sufficiently high to be considered, and predict the most important class (majority rule)
  3. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{\theta_i\} where \theta_i=\cos(\boldsymbol{p}_{\cdot},\boldsymbol{p}_{i})=\displaystyle{\frac{\boldsymbol{p}_{\cdot}\cdot\boldsymbol{p}_{i}}{\|\boldsymbol{p}_{\cdot}\|\|\boldsymbol{p}_{i}\|}}
  4. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{KL_{\cdot|i}\} and KL_{\cdot|i}=\displaystyle{\sum_{j=1}^n p_{\cdot,j}\log\frac{p_{\cdot,j}}{p_{i,j}}} (but one can select another metric)
  5. \widehat{y}_{\cdot}=y_{j^\star} where {j^\star}=\underset{j\in\mathcal{J}}{\text{argmax}}\{p_{\cdot,j}\} and \mathcal{J} is a sample of k observations, chosen randomly, one in each group (one-shot procedure): the predicted class is the one of the observation the most likely to be other same class

Heuristically, it can be related to some k nearest neighbors strategy: we give the attribute that most neighbors have. The total distance is a weighted sum of the componentwise distances (for the logistic regression).

  • Simulation study

In order to test that technique, let us generate some multinomial model where y has 10 possible labels, with 6 (independent) covariates x_1,\cdots,x_6, and \mathbb{P}[Y=A_k|\boldsymbol{X}=\boldsymbol{x}]\propto \exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k] (where coefficients \boldsymbol{\beta}_k where generated randomly) for k\in\{1,2,\cdots,10\} (there were 10 categories) and with n=700 observations.

n=700
X1=rnorm(n)
X2=rnorm(n)
X3=rnorm(n)
X4=rnorm(n)
X5=rnorm(n)
X6=rnorm(n)
X=cbind(1,X1,X2,X3,sqrt(abs(X4)),X5*X1,X6)
k = 10
 PARAM = matrix(rnorm(k*6),k,6)
 PARAM[,1]=PARAM[,1]-1
 PARAM=cbind(PARAM,0)
 P=matrix(NA,n,k-1)
 for(j in 1:(k-1)) P[,j] = X %*% (PARAM[j,])+rnorm(n)
 P=cbind(P,0)
S=apply(exp(P),1,sum)
Pb = exp(P)/S
tirage = function(i){
      sample(1:10,size=1,prob = Pb[i,])
}
Y = LETTERS[Vectorize(tirage)(1:n)]
dbase = data.frame(Y=as.factor(Y),X1,X2,X3,X4,X5)

In the paragraph previously, I suggested to take the most likely one. Being wrong means that it was not the first choice. But perhaps being the second or the third choice is not that bad, actually. So in my simulations, I look at the proportion of predictions were our prediction is the good one (top 1), if the true one is either the most likely or the second most likely (top 2), or in the top 3. That will be on my x-axis. I draw some line, but we simply have three points (top 1, top 2 and top 3). I compute the proportion of good prediction, using cross-validation techniques (10-fold). The black lines are one of the methods described above. The red one is the standard multinomial model (with a logistic link function). For the Siamese model, I tried several models. I tried a logistic regression, and some smooth version (GAM) on top

and a classification tree, on the left, as well as some random forest on the right, below.

It looks like the multinomial approach performs always better than any Siamese one… and to be honest, I am disappointed.

Here is the code I did use when I considered a logistic regression on the Siamese dataset,

set.seed(1)
kfold = sample(rep(1:10,n/10))
 
KFOLDglm = function(i){
i_test=which(kfold==i)
i_calibration=which(kfold!=i)
y=credit[i_calibration,"Y"]
tirage = function(){
v=c(sample(i_calibration[y==levels(y)[1]],size=1),
    sample(i_calibration[y==levels(y)[2]],size=1),
    sample(i_calibration[y==levels(y)[3]],size=1),
    sample(i_calibration[y==levels(y)[4]],size=1),
    sample(i_calibration[y==levels(y)[5]],size=1),
    sample(i_calibration[y==levels(y)[6]],size=1),
    sample(i_calibration[y==levels(y)[7]],size=1),
    sample(i_calibration[y==levels(y)[8]],size=1),
    sample(i_calibration[y==levels(y)[9]],size=1),
    sample(i_calibration[y==levels(y)[10]],size=1))
names(v)=levels(y)
return(v)
}
 
LogisticModel <- multinom(Y ~ ., data = credit[i_calibration,], trace=FALSE)
 
comparaisonx = function(base,x=base[1,]){
  mix_base = base
  for(j in 1:ncol(base)){
    xj = as.numeric(x[j])-base[,j]
    mix_base[,j] = (xj)
  }
  mix_base
}
comparaisony = function(base,y=base[1]){
  as.factor(base == y)
}
creditx = credit[,-which(names(credit) == "Y")]
nc=length(i_calibration)
B=comparaisonx(base = creditx[i_calibration[2:nc],],x=creditx[i_calibration[1],])
B$Y=comparaisony(base = credit[i_calibration[2:nc],"Y"],y=credit[i_calibration[1],"Y"])
for(i in 2:(nc-1)){
  B0=comparaisonx(base = creditx[i_calibration[(i+1):nc],],x=creditx[i_calibration[i],])
  B0$Y=comparaisony(base = credit[i_calibration[(i+1):nc],"Y"],y=credit[i_calibration[i],"Y"])
  B=rbind(B,B0)
}
credit_mix = B
 
OneShotLogisticModel <- glm(Y ~ ., data = credit_mix, family=binomial)
A_ref = table(credit[i_calibration,"Y"])/length(i_calibration)
 
vect_oneshot = function(i){
  B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
  predict(OneShotLogisticModel,type="response",newdata=B2)
}
 
prediction_oneshot = function(i,type=1){
B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
y=credit[i_calibration,"Y"]
base = data.frame(p,y)
base = base[rev(order(base$p)),]
if(type==1){T = table(base$y[1:11])
return(names(which.max(T)))}
if(type==2){return(base$y[1])}
if(type==3){A=table(base$y[1:10])/10
T=A/A_ref
return(names(which.max(T)))}
if(type==4){
  costheta = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    costheta[j] = sum(vecteur_proba*p)/(sqrt(sum(vecteur_proba^2))*sqrt(sum(p^2)))
  }
  return(y[which.max(costheta)])}
if(type==5){
  kl = rep(NA,length(i_calibration))
  for(j in 1:length(i_calibration)){
    vecteur_proba = vect_oneshot(i_calibration[j])
    kl[j] = sum(p*log(vecteur_proba/p))
  }
  return(y[which.max(as.vector(kl))])}
if(type==6){ ## one shot : tirer au hasard un de chaque, et dire lequel est plus credible !
 
y=credit[i_calibration,"Y"]
tirage = function(){
    v=c(sample(i_calibration[y==levels(y)[1]],size=1),
        sample(i_calibration[y==levels(y)[2]],size=1),
        sample(i_calibration[y==levels(y)[3]],size=1),
        sample(i_calibration[y==levels(y)[4]],size=1),
        sample(i_calibration[y==levels(y)[5]],size=1),
        sample(i_calibration[y==levels(y)[6]],size=1),
        sample(i_calibration[y==levels(y)[7]],size=1),
        sample(i_calibration[y==levels(y)[8]],size=1),
        sample(i_calibration[y==levels(y)[9]],size=1),
        sample(i_calibration[y==levels(y)[10]],size=1))
    names(v)=levels(y)
    return(v)
  }
  pd=rep(NA,101)
for(ix in 1:101){
ids = tirage()
B2=comparaisonx(base = creditx[ids,],x=creditx[i,])
p=predict(OneShotLogisticModel,type="response",newdata=B2)
pd[ix]=levels(y)[which.max(p)]
}
levels(y)[which.max(table(pd))]
}
}
PRED0=as.character(credit[i_test,"Y"])
PRED1=as.character(predict(LogisticModel,type = "class",
                           newdata=credit[i_test,]))
PRED21=as.character(Vectorize(function(i) prediction_oneshot(i,type=1))
                    (i_test))
PRED22=as.character(Vectorize(function(i) prediction_oneshot(i,type=2))(i_test))
PRED23=as.character(Vectorize(function(i) prediction_oneshot(i,type=3))(i_test))
PRED24=as.character(Vectorize(function(i) prediction_oneshot(i,type=4))(i_test))
PRED25=as.character(Vectorize(function(i) prediction_oneshot(i,type=5))(i_test))
PRED26=as.character(Vectorize(function(i) prediction_oneshot(i,type=6))(i_test))
B=data.frame(PRED0,PRED1,PRED21,PRED22,PRED23,PRED24,PRED25,PRED26)
B
}
 
s=1/100;setTxtProgressBar(pb, s*2)
for(i in 2:10){
  PREDICTION = rbind(PREDICTION,KFOLDglm(i))
s=s+1/100;setTxtProgressBar(pb, s*2)}
for(j in 1:5) PREDICTION[,j]=as.character(PREDICTION[,j])
L=list()
v=mean(PREDICTION[,1]!=PREDICTION[,2])
names(v)="logistic"
L[["logistic"]]=v
v=c(mean(PREDICTION[,1]!=PREDICTION[,3]),
  mean(PREDICTION[,1]!=PREDICTION[,4]),
  mean(PREDICTION[,1]!=PREDICTION[,5]),
  mean(PREDICTION[,1]!=PREDICTION[,6]),
  mean(PREDICTION[,1]!=PREDICTION[,7]),
  mean(PREDICTION[,1]!=PREDICTION[,8]))
names(v)=c("top10","max","10norm","cos","KL","OS")
L[["glm"]]=v

Autocalibration for Insurance Pricing with Machine Learning

With Michel Denuit and Julien Trufin, we recently uploaded a joint paper on ArXiv, entitled Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning.

Boosting techniques and neural networks are particularly effective machine learning methods for insurance pricing.
Often in practice, there are nevertheless endless debates about the choice of the right loss function to be used to train the machine learning model, as well as about the appropriate metric to assess the performances of competing models. Also,
the sum of fitted values can depart from the observed totals to a large extent and this often confuses actuarial analysts.
The lack of balance inherent to training models by minimizing deviance outside the familiar GLM with canonical
link setting has been empirically documented in Wüthrich (2019, 2020) who attributes it to the early stopping rule in gradient descent methods for model fitting.
The present paper aims to further study this phenomenon when learning proceeds by minimizing Tweedie deviance.
It is shown that minimizing deviance involves a trade-off between the integral of weighted differences of lower partial moments and the bias measured on a specific scale.
Autocalibration is then proposed as a remedy. This new method to correct for bias adds an extra local GLM step to the analysis.
Theoretically, it is shown that it implements the autocalibration concept in pure premium calculation and ensures that balance also holds on a local scale, not only at portfolio level as with existing bias-correction techniques. The convex order appears to be the natural tool to compare competing models, putting a new light on the diagnostic graphs and associated metrics proposed by Denuit et al. (2019).

In this paper, we started with a simple observation : with GLM (like a Poisson regression) we have unbiased predictions, in the sense that \widehat{y}_1+...+\widehat{y}_n=y_1+...+y_n (on the training dataset, we should expect to have \widehat{y}_1+...+\widehat{y}_n\approx y_1+...+y_n on the validation dataset). But with any machin learning algorithm, this is no longer the case. For instance, with a boosting algorithm, we can end up with a significant bias. This is an application on an insurance dataset from the CASdataset package, with a GLM (Poisson regression) on the left, and additive smooth version, GAM in the middle and some boosting algorithm on the right. Here is the dispersion of the predictions, \widehat{y}_i‘s,

Let \widehat{\pi} denote a fited model. If \widehat{\pi} was close to the true parameter of the Poisson model, \mu where \mu(\boldsymbol{x})=\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}], then function  s\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=s] should be very close to the first diagonal (since \mathbb{E}[Y|\mu(\boldsymbol{X})=s]=s). And this is was we observe below, for the GLM. But the boosting algorithm has a significant bias.

Furthermore, the bias is not the same everywhere. We can see it more precisely on a quantile version of the x-axis

that is u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]. Here is a multiplicative version of that bias u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)]/u

The idea of autocalibration is to use that multiplicative factor to correct from the bias. Thus, if \widehat{\pi}(\boldsymbol{x}) is close to the 75% upper quantile, then the true value should be 20% larger for the bootsting algorithm. Here is the new distribution of the predictions, with that correction

and we can observe that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] is now much closer to the first diagonal (below is the multiplicative bias)

We discuss in the paper this correction for local bias, so that u\mapsto\mathbb{E}[Y|\widehat{\pi}(\boldsymbol{X})=F_{\widehat{\pi}}^{-1}(u)] becomes as close as possible to the first diagonal… see https://arxiv.org/abs/2103.03635 for a complete version of the paper, and https://github.com/freakonometrics/autocalibration/ for the complete version of the R codes

Principal Component Analysis: A Generalized Gini Approach

Our paper, with Stéphane Mussard and Téa Ouraga, entitle Principal Component Analysis: A Generalized Gini Approach is finally out in the European Journal of Operations Research.

A principal component analysis based on the generalized Gini correlation index is proposed (Gini PCA). The Gini PCA generalizes the standard PCA based on the variance. It is shown, in the Gaussian case, that the standard PCA is equivalent to the Gini PCA. It is also proven that the dimensionality reduction based on the generalized Gini correlation matrix, that relies on city-block distances, is robust to outliers. Monte Carlo simulations and an application on cars data (with outliers) show the robustness of the Gini PCA and provide different interpretations of the results compared with the variance PCA.

Some general thoughts on Partial Dependence Plots with correlated covariates

The partial dependence plot is a nice tool to analyse the impact of some explanatory variables when using nonlinear models, such as a random forest, or some gradient boosting.The idea (in dimension 2), given a model m(x_1,x_2) for \mathbb{E}[Y|X_1=x_1,X_2=x_2]. The partial dependence plot for variable x_1 is model m is function p_1 defined as x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)]. This can be approximated, using some dataset using \widehat{p}_1(x_1)=\frac{1}{n}\sum_{i=1}^n m(x_1,x_{2,i})My concern here what the interpretation of that plot when there are some (strongly) correlated covariates. Let us generate some dataset to start with

n=1000
library(mnormt)
r=.7
set.seed(1234)
X = rmnorm(n,mean = c(0,0),varcov = matrix(c(1,r,r,1),2,2))
Y = 1+X[,1]-2*X[,2]+rnorm(n)/2
df = data.frame(Y=Y,X1=X[,1],X2=X[,2])

As we can see, the true model is here is y_i=\beta_0+\beta_1 x_{1,i}+\beta_2x_{2,i}+\varepsilon_i where \beta_1 =1 but the two variables are positively correlated, and the second one has a strong negative impact. Note that here

reg = lm(Y~.,data=df)
summary(reg)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.01414    0.01601   63.35   <2e-16 ***
X1           1.02268    0.02305   44.37   <2e-16 ***
X2          -2.03248    0.02342  -86.80   <2e-16 ***

If we estimate a wrongly specified model y_i=b_0+b_1 x_{1,i}+\eta_i, we would get

reg1 = lm(Y~X1,data=df)
summary(reg1)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03522    0.04680  22.121   <2e-16 ***
X1          -0.44148    0.04591  -9.616   <2e-16 ***

Thus, on the proper model, \widehat{\beta}_1\sim+1.02 while \widehat{b}_1\sim-0.44  on the mispecified model.

Now, let us look at the parial dependence plot of the good model, using standard R dedicated packages,

library(pdp) 
pdp::partial(reg, pred.var = "X1", plot = TRUE,
              plot.engine = "ggplot2")

which is the linear line y=1+x, that corresponds to y=\beta_0+\beta_1x.

library(DALEX)
plot(DALEX::single_variable(DALEX::explain(reg,
data=df),variable = "X1",type = "pdp"))

which corresponds to the previous graph. Here, it is also possible to creaste our own function to compute that partial dependence plot,

pdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  mean(predict(reg,newdata=nd))
}

that will be the straight line below (the dotted line is the theoretical one y=1+x,

vx=seq(-3.5,3.5,length=101)
vpdp1 = Vectorize(pdp1)(vx)
plot(vx,vpdp1,type="l")
abline(a=1,b=1,lty=2)

which is very different from the univariate regression on x_1

abline(reg1,col="red")

Actually, the later is very consistent with a local regression, only on x_1

library(locfit)
lines(locfit(Y~X1,data=df),col="blue")

Now, to get back to the definition of the partial dependence plot, x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)], in the context of correlated variable, I was wondering if it would not make more sense to consider some local version actually, something like x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2|X_1}}[m(x_1,X_2)]. My intuition was that, somehow, it did not make any sense to consider any X_2 while X_1 was fixed (and equal to x_1). But it would make more sense actually to look at more valid X_2‘s given the value of X_1. And a natural estimate could be some k neareast-neighbors, i.e. \tilde{p}_1(x_1)=\frac{1}{k}\sum_{i\in\mathcal{V}_k(x)}^n m(x_1,x_{2,i}) where \mathcal{V}_k(x) is the set of indices of the k x_i‘s that are the closest to x, i.e.

lpdp1 = function(x1){
  nd = data.frame(X1=x1,X2=df$X2)
  idx = rank(abs(df$X1-x1))
  mean(predict(reg,newdata=nd[idx<50,]))
}
vlpdp1 = Vectorize(lpdp1)(vx)
lines(vx,vlpdp1,col="darkgreen",lwd=2)

Surprisingly (?), this local partial dependence plot gives a curve that corresponds to the simple regression…

Régression sur des variables corrélées, un peu de géométrie

Depuis quelques jours, je mets en ligne des billets sur la régression sur des variables corrélés, en commençant par du deuxième effet kiss-cool suivi d’un court billet sur des régressions en cascade. Oui, depuis le premier, mon collègue Olivier ne cesse de demander des compléments. Alors ce soir, on va faire des dessins…

On a ici 3 vecteurs, \color{red}{\vec{x_1}}, \color{green}{\vec{x_2}} et \color{blue}{\vec{y}}. Pour rappel, la longueur de mes vecteurs, c’est la variance de la variable associée, \|{\color{red}{\vec{x_1}}}\|=\text{Var}({\color{red}{x_1}}). Et commençons par supposer que les variables \color{red}{x_1} et \color{green}{x_2} sont indépendantes (ou plus simplement non-corrélées) soit {\color{red}{\vec{x_1}}}\perp\color{green}{\vec{x_2}}. C’est le dessin ci-dessous,

La régression linéaire multiple est ici équivalente aux régressions simples, comme on en parlait dans les billets précédant. Géométriquement, on voit que les coordonnées de la projection de \color{blue}{\vec{y}} sur l’espace engendré par \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} (c’est à dire le plan (x,y)) correpond au point noir. C’est la régression multiple. Les deux régressions simples sont obtenus en projetant respectivement sur sur l’espace engendré par \color{red}{\vec{x_1}} (c’est à dire la droite (x)) et sur \color{green}{\vec{x_2}} (c’est à dire la droite (y)).

Supposons maintenant que \color{red}{\vec{x_1}} et \color{green}{\vec{x_2}} ne sont plus orthogonaux. Plus précisément, on va garder \color{red}{\vec{x_1}} et \color{blue}{\vec{y}} inchangés, et on va changer , \color{green}{\vec{x_2}}, tout en gardant la norme \|\color{green}{\vec{x_2}}\| inchangée. C’est ce qu’on a sur le schéma ci-dessous

On avait vu qu’on ne pouvait plus faire deux régressions indépendantes, mais qu’on pouvait garder la première, et en faire ensuite une seconde, mais en faisant un peu attention. Pour la première régression, on va projeter \color{blue}{\vec{y}} sur \color{red}{\vec{x_1}} et ca c’est facile, rien n’a changé,

Pour la seconde étape, rappelons qu’on doit projeter à la fois \color{blue}{\vec{y}} et \color{green}{\vec{x_2}} sur l’orthogonal de \color{red}{\vec{x_1}}, c’est à dire ici le plan (y,z), en rouge ci-dessous

(je renvois au précédant billet pour les formules algébriques). On obtient alors les projections suivantes, qu’on pourrait noter \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} et \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}}

On le voit, si on regarde sur le mur du fond, en rouge, la projection de \Pi_{\color{red}{x_1^{\perp}}}\color{blue}{\vec{y}} sur \Pi_{\color{red}{x_1^{\perp}}}\color{green}{\vec{x_2}} est un peu plus grande, autrement dit, le coefficient de la régression simple sera un peu plus grand. En fait, si on veut faire des calculs, on voit que le sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle va intervenir ici… Plus l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle sera petit, plus les variables seront corrélées, et plus la correction sera importante. Graphiquement, la correction dont je parle, c’est le petit segment, à gauche sur le dessin

Mais simulons quelques données pour visualiser tout ça…

library(mnormt)
n=10000
r=0
set.seed(1)
Z=cbind(rnorm(n)*sqrt(2),rnorm(n)*sqrt(3))
df=data.frame(y=Z[,1]+Z[,2]+rnorm(n)*sqrt(3),x1=Z[,1],x2=Z[,2])

On simule ici des variables explicatives indépendamment. Et on peut vérifier que les distances sont bien celles de notre dessin

var(df$x1)
[1] 2.049731
var(df$x2)
[1] 2.810419
var(df$y)
[1] 8.066665

et en plus, les variables sont quasiment orthogonales

cor(df$x1,df$x2)
[1] 0.004845078

Si on fait la régression multiple

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  0.99904    0.01219   81.94   2e-16 ***
x2  1.00090    0.01017   98.40   2e-16 ***

On a des coefficients qui sont identiques (à très peu de choses près, mais le soucis vient du fait que j’ai enlevé la constante du modèle pour coller au mieux avec mon dessin)

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

et pour la seconde régression aussi

reg2=lm(y~0+x2,data=df)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x2  1.00496    0.01315   76.42   2e-16 ***

On a déjà parlé de ce résultat longuement (je renvoie au précédant billet, des régressions en cascade). Regardons maintenant ce qui se passe si on garde nos observations \color{blue}{y_i}, \color{red}{x_{1,i}} mais que maintenant, la variable \color{red}{x_{2,i}} se retrouve corrélée avec \color{red}{x_{1,i}} , disons avec une corrélation de l’ordre de 0.4

r=.4
df$x2=r*df$x1+sqrt(1-r^2)*df$x2
cor(df$x1,df$x2)
[1] 0.3461496

Faire la régression multiple  donne ici

reg=lm(y~0+x1+x2,data=df)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1   0.5622     0.0130   43.26   2e-16 ***
x2   1.0921     0.0111   98.40   2e-16 ***

ce qui est différent des régressions simples, comme la première

reg1=lm(y~0+x1,data=df)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
x1  1.00489    0.01711   58.75   2e-16 ***

Ce que nous dit le théorème de Frisch-Waugh, c’est qu’on peut faire des régressions en cascades: on peut garder cette première régression, mais pour la seconde, on va régresser les projections sur l’orthogonal à \color{red}{\vec{x_1}} (le mur en rouge au fond dans notre dessin). Ici, on obtient

reg2bis=lm(residuals(lm(y~x1,data=df))~residuals(lm(x2~x1,data=df)))
summary(reg2bis)
 
Coefficients:
           Estimate Std. Error t value Pr(|t|)      
residuals(2) 1.0921     0.0111    98.4   2e-16 ***

qui donne effectivement un coefficient de régression un peu supérieur à celui que nous avions auparavant, comme nous l’observions sur le dessin. En fait on peut remarquer que la différence relative s’obtient très facilement avec la corrélation

1/sqrt(1-r^2)
[1] 1.091089

et pour rappel, d’un point de vue géométrique, la corrélation est un cosinus d’un angle, r=\cos(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle). Autrement dit \sqrt{1-r^2}=\sin(\langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle), c’est à dire que la correction est directement liée au sinus de l’angle \langle\color{red}{\vec{x_1}},\color{green}{\vec{x_2}}\rangle, comme mentionné auparavant… amusant, non?