Category Archives: Statistics

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?

Des régressions en cascade

Cette fin de semaine, je mettais en ligne un court billet du deuxième effet kiss-cool où je rappelais que quand on fait une régression sur plusieurs variables explicatives corrélées, ce n’est pas équivalent à faire plusieurs régressions simple. C’est ce que disait le théorème de Frisch-Waugh, publié il y a près de 90 ans. Mais je glissais qu’il était possible de faire des régressions en cascade, sans aller beaucoup plus loin. Il faut relire l’article de Michael Lovell qui donner une interprétation géométrique de ce dessin. Mais programmons le, pour le voir… Considérons une base avec 3 variables explicatives, correspondant à des nombres d’incendies à Chicago

chicago=read.table("http://freakonometrics.free.fr/chicago.txt", header=TRUE,sep=";")

La régression multiple s’écrit iciy_i=\beta_0+\beta_1 x_{1,i}+\beta_2 x_{2,i}+\beta_3 x_{3,i}+\varepsilon_iOn peut estimer ces paramètres par moindre carrés,

reg=lm(Fire~.,data=chicago)
summary(reg)
 
Coefficients:
            Estimate Std. Error t value  Pr(|t|)    
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X_1         -0.62764    5.28130  -0.119 0.905953    
X_2          0.22378    0.06161   3.632 0.000744 ***
X_3         -1.55059    0.38195  -4.060 0.000204 ***

et on a alors la prévision suivante
\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 x_{1,i}+\widehat{\beta}_2 x_{2,i}+\widehat{\beta}_3 x_{3,i}En fait, je peux même visualiser cette prévision.

BC = cbind(x0=coefficients(reg)[1],
           x1=coefficients(reg)[2]*chicago[,2],
           x2=coefficients(reg)[3]*chicago[,3],
           x3=coefficients(reg)[4]*chicago[,4],
           yp=predict(reg),y=chicago$Fire)
colr = 1:4
dessin_fire = function(i,BC=BC){
plot(0:3,0:3,xlim=c(0,35),col="white",ylim=c(-1,3),xlab="",ylab="",axes=FALSE)
abline(v=BC[i,1],lty=2)
rect(BC[i,1],3,BC[i,1]+BC[i,2],2,col=colr[1],border=NA)
arrows(BC[i,1],2.5,BC[i,1]+BC[i,2],2.5,lwd=2,length=.1,col="white")
rect(BC[i,1]+BC[i,2],2,BC[i,1]+BC[i,2]+BC[i,3],1,col=colr[2],border=NA)
arrows(BC[i,1]+BC[i,2],1.5,BC[i,1]+BC[i,2]+BC[i,3],1.5,lwd=2,length=.1,col="white") 
rect(BC[i,1]+BC[i,2]+BC[i,3],1,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],0,col=colr[3],border=NA)
segments(BC[i,5],0,BC[i,5],1,lwd=3)
arrows(BC[i,1]+BC[i,2]+BC[i,3],.5,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],.5,lwd=2,length=.1,col="white") 
abline(v=BC[i,1],lty=2)
rect(BC[i,5],-1,BC[i,6],0,col=colr[4],density=20,border=NA)
segments(BC[i,6],0,BC[i,6],-1,lwd=3,col=colr[4])
text(26,2.5,expression(X[1]),pos=4,col=colr[1])
text(26,1.5,expression(X[2]),pos=4,col=colr[2])
text(26,.5,expression(X[3]),pos=4,col=colr[3])
text(26,-.5,expression(epsilon),pos=4,col=colr[4])
axis(1)}

Par exemple, pour la douzième observation de ma base

dessin_fire(12)

on obtient

On lit de haut en bas : on commence par constante, vers 22. Puis on rajoute trois fois rien, parce que x_1 était non-significative dans notre régression. Puis on rajoute un petit quelque chose à cause de x_2, pour passer à 25, et x_3 nous fait plonger de plus de 20 points, pour finir vers 3.67.Le trait noir est la prévision que l’on obtient, à partir des trois variables. En bleu, on peut même voir l’erreur.

On peut noter que x_1 est non-significative, tout en étant relativement corrélée avec y

cor(chicago[,1:2])
          Fire       X_1
Fire 1.0000000 0.3773486
X_1  0.3773486 1.0000000

Car si x_1 est non-significative, c’est parce que cette variable est très corrélée avec une autre variable explicative. En fait, on peut faire une série de régressions en cascade, en commençant par x_1, et plus spécifiquement, on peut écrire\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_1 x_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,2}+\widehat{b}_2 \tilde{x}_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}où le premier terme (1) est obtenu en régression simplement y sur x_1 (oui, on fait juste une régression simple),  y_i=b_0+b_1 x_{1,i}+\eta_ipuis on apporte une petite correction, pour tenir compte de ce que nous apprend x_2 une fois la première régression effectuée. C’est exactement ce que raconte le théorème de Frish-Waugh : on va projeter y sur x_1^{\perp} (car c’est justement ce que x_1 n’explique pas), et x_2 sur x_1^{\perp}  et faire la régression de ces deux projections, soit, pour (2) \Pi_{x_1^{\perp}}y_i=b_{0,2}+b_2 \Pi_{x_1^{\perp}}x_{2,i}+\eta_{2,i}puis pour (3) \Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}En terme informatique, cela donne

reg1=lm(Fire~X_1,data=chicago)
reg2=lm(residuals(lm(Fire~X_1,data=chicago))~residuals(lm(X_2~X_1,data=chicago)))
reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))

soit ici

BC123=cbind(x0=coefficients(reg1)[1],
            x1=coefficients(reg1)[2]*chicago[,2],
            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_2~X_1,data=chicago)),
            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),
            yp=predict(reg),y=chicago$Fire)

En adaptant la fonction de dessin précédante on obtient

dessin_fire_123(12)

Autrement dit, on commence avec une régression simple, en rouge: on part d’une constante vers 3.5, puis on augemente notre prévision, en tenant compte de x_1. Puis on corrige avec x_2 (ou plutôt la partie de x_2 non expliquée par x_1). sur ce qui n’avait pas été expliqué auparavant, ce qui pousse à revoir un peu à la baisse notre prévision, puis on tient compte de x_3 (ou là encore, sa projection sur l’orthogonal de x_1 et x_2). On note que la prévision est très exactement la même qu’auparavant.

Mais on peut aller plus loin… pourquoi commencer par x_1? On peut aussi commencer par x_2. On considère alors\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_2 x_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,1}+\widehat{b}_1 \tilde{x}_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}où le premier terme (2) est obtenu en régression simplement y sur x_2y_i=b_0+b_2 x_{2,i}+\eta_ipuis on apporte une petite correction, pour tenir compte de ce que nous apprend x_1 une fois la première régression effectuée, (1) \Pi_{x_2^{\perp}}y_i=b_{0,1}+b_1\Pi_{x_2^{\perp}}x_{1,i}+\eta_{1,i}puis pour (3) \Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i} Soit

reg1=lm(Fire~X_2,data=chicago)
reg2=lm(residuals(lm(Fire~X_2,data=chicago))~residuals(lm(X_1~X_2,data=chicago)))
reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago)))
BC213=cbind(x0=coefficients(reg1)[1],
            x1=coefficients(reg1)[2]*chicago[,3],
            x2=coefficients(reg2)[1]+coefficients(reg2)[2]*residuals(lm(X_1~X_2,data=chicago)),
            x3=coefficients(reg3)[1]+coefficients(reg3)[2]*residuals(lm(X_3~X_1+X_2,data=chicago)),
            yp=predict(reg),y=chicago$Fire)

et visuellement, on obtient

dessin_fire_213(12)

Cette fois, la constante est un peu plus élevée, et on commence par faire une régression seulement avec x_2, pour arriver un peu en dessous de 10. Puis on corrige. On note que cette seconde correction nous ramène… très exactement comme dans le cas précédant. Mais si on pense en projection successives, on ne doit pas être surpris.

On peut récapituler avec une autre prévision, soit avec le modèle de régresion multiple (mais là encore, x_1 n’explique pas grand chose)

dessin_fire(15)

soit en commencant par une régression simple sur x_1 puis en faisant des régression en cascade,

dessin_fire_123(15)

soit, alternativement, en commencant par une régression simple sur x_2

dessin_fire_213(15)

Comme annoncé, cette trois approches sont équivalentes, et donnent très exactement la même prévision.

Regression discontinuity model for TV series

In September, we are usually happy to see our favorite TV series back on air… Or not? Because, admit it, if we are happy to see those characters back, most of the time, we are disappointed, too. So why not look at the data, to confirm this feeling? Nazareno Andrade shared some nice codes to get IMDB ratings in a nice csv file (you can either use the large csv file, or run your own codes)

download.file("https://github.com/nazareno/imdb-series/raw/master/data/series_from_imdb.csv",
destfile="series_from_imdb.csv")
base = read.csv("series_from_imdb.csv")

It is a large dataset, with more than 64,000 episodes of almost 890 TV series,

str(base)
'data.frame':	64018 obs. of  18 variables:
 $ series_name: Factor w/ 889 levels "'Allo 'Allo!",..: 137 137 137 137 137 137 137 137 137 137 ...
 $ episode    : Factor w/ 54090 levels "-30-","¡Viva los muertos!",..: 32314 7446 16 7176 17748 9562 1379 36218 17845 5553 ...
 $ series_ep  : int  1 2 3 4 5 6 7 8 9 10 ...
 $ season     : int  1 1 1 1 1 1 1 2 2 2 ...
 $ season_ep  : int  1 2 3 4 5 6 7 1 2 3 ...
 $ user_rating: num  8.9 8.7 8.7 8.2 8.3 9.2 8.8 8.7 9.2 8.3 ...

Just pick a TV series, for instance Dan Harmon’s Community,

sbase = base[base$series_name=="Community",]

We can plot the evolution of the rating over the 110 episodes.

sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),]
sbase$series_ep=1:nrow(sbase)

()since there could be some problem with the data (such as duplicates, let us clean it quickly)

plot(sbase$series_ep,sbase$UserRating,xlab=sbase$series_name[1])
idx=c(0,which(diff(sbase$season)!=0),nrow(sbase))
abline(v=idx+.5,lty=2,col=colr[2])
a = unique(sbase$season)
for(u in a){
  ssbase = sbase[sbase$season==u,]
  reg = lm(UserRating~series_ep,data=ssbase)
  lines(ssbase$series_ep,predict(reg),col=colr[3],lwd=2)
}

The vertical lines are here to visualize the seasons. On issue is that the lenght can vary with time. Consider Linwood Boomer’s Malcom in The Middle,

sbase = base[base$series_name=="Malcolm in the Middle",]

or Craig Thomas and Carter Bays’s How I Met Your Mother,

sbase = base[base$series_name=="How I Met Your Mother",]

On those two, the evolution is rather stable. Look at AMC’s The Walking Dead,

sbase = base[base$series_name=="The Walking Dead",]

Now, look at Howard Gordon and Alex Gansa’s Homeland,

sbase = base[base$series_name=="Homeland",]

There is an issue here with the last episode of season4, “Long Time Coming“, that has a very poor rating. If we remove that point, we get the thin line. Note that the regression line is always increasing. For Michael Hirst’s Vickings, we have

sbase = base[base$series_name=="Vicking",]

If we look more carefully on the previous graph, for five seasons (out of six), we have a positive slope. Well, to be honest, it is not significantly positive most of the time, but still. Out of 80 shows, and a total of 583 seasons, the slope is postive 75% of the time (433) and negative 25% of the time (150).

BASE = NULL
L80 = unique(base$series_name)
for(j in 1:length(L)){
sbase=base[base$series_name==L[j],]
sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),]
sbase=sbase[sbase$season>0,]
sbase$series_ep=1:nrow(sbase)
a=unique(sbase$season)
a=a[!is.na(a)]
for(u in a){
  ssbase=sbase[sbase$season==u,]
  reg=lm(UserRating~series_ep,data=ssbase)
  pente = NA
  if((!is.na(coefficients(reg)[2]))&(!is.na((summary(reg)$coefficients[2,4])))){
  if((summary(reg)$coefficients[2,4]<.05)&(coefficients(reg)[2]>0)) pente="positive"
  if((summary(reg)$coefficients[2,4]<.05)&(coefficients(reg)[2]<0)) pente="negative" sdf=data.frame(nom=sbase$series_name[1],season=u,slope=coefficients(reg)[2],inf=confint(reg)[2,1],sup=confint(reg)[2,2],signe=pente) BASE=rbind(BASE,sdf)} }} str(BASE) 'data.frame': 583 obs. of 6 variables: $ nom : Factor w/ 80 levels "Friends","Game of Thrones",..: 1 1 1 1 1 1 1 1 1 1 ... 
 
mean(BASE$slope>0)
[1] 0.7427101
table(BASE$signe)
negative positive 
      15      144

Most of the time, the slope is not significant. To be more specific, 72% of the time, the slope is not significant. But when it is, 90% of the time, it is positive (144 seasons). Let us look at other TV series, for instance Joel Surnow and Robert Cochran’s 24,

sbase = base[base$series_name=="24",]

Álex Pina’s La Casa de Papel,

sbase = base[base$series_name=="La Casa de Papel",]

Steven Knight’s Peaky Blinders,

sbase = base[base$series_name=="Peaky Blinders",]

or David Simon’s The Wire,

sbase = base[base$series_name=="The Wire",]

The slope is increasing over almost all seasons. But a major drawback is that when we get back to our show, for a new season, we usually get disapointed. More specifically, we can quantify the difference in red below

that can be estimated using

sbase12 = sbase[sbase$season%in%c(a[ij],a[ij+1]),]
seuil = sbase12$series_ep[which(diff(sbase12$season)!=0)]+.5
s = function(x) (x-seuil)*(x>seuil)
reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep>seuil),data=sbase12)

Here we have

summary(reg)
Coefficients:
                         Estimate Std. Error t value  Pr(|t|)    
(Intercept)               8.45000    0.16338  51.719    2e-16 ***
series_ep                 0.10000    0.03235   3.091 0.008598 ** 
s(series_ep)              0.02000    0.04218   0.474 0.643291    
I(series_ep)TRUE.        -1.01778    0.20486  -4.968 0.000257 ***

so the drop of 1 point (out of 10) cannot be claimed as being significant. That is the idea of regression discontinuity.

If we loop again over all our series, we have 485 pairs of consecutive seasons. As expected, in 75% of the casse, from season t-1 to season t, we observe a negative rupture. As previously, in 70% of the cases, it is not significat (with linear models before and after), and when it is significant, it is negative in 96% of the cases ! But an alternative can be to use nonparametric models, on both sides.

To illustrate, consider David Benioff and D. B. Weiss’s Game of Thrones,

sbase = base[base$series_name=="Game of Thrones",]

but let us remove the last season (no spoiler here, but clearly not worst watching)

Consider for instance the drop between season 1 and season 2,

library(rdd)
sbase12=sbase[sbase$season%in%c(1,2),]
lmr=RDestimate(UserRating~series_ep,data=sbase12,cutpoint=mean(range(sbase12$series_ep)))
plot(lmr)

This is very consistent with what we observed with our linear regressions actually,

seuil=10.5
s = function(x) (x-seuil)*(x>seuil)
reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep>seuil),data=sbase12)
summary(reg)
 
Coefficients:
                         Estimate Std. Error t value  Pr(|t|)    
(Intercept)               8.70000    0.15458  56.281    2e-16 ***
series_ep                 0.07273    0.02491   2.919  0.01003 *  
s(series_ep)              0.01455    0.03523   0.413  0.68520    
I(series_ep)TRUE         -0.94000    0.20316  -4.627  0.00028 ***

Here, the drop of one point is significant…

So, your favorite show had an outstanding finale ? and you can’t wait to watch the new season… Well, statistically, it’s very likely that you will be disapointed by the first episode of the forthcoming season…

Du deuxième effet kiss-cool (régression multiple, scoring et évaluation)

Lorsque j’étais petit (il y a fort longtemps, à une époque où je regardais pas mal la télévision) il y avait une publicité pour les pastilles kiss cool,

Et quand je présente la régression multiple à mes étudiants, je ne peux m’empêcher d’y penser… Mais avant d’aller plus loin sur le parallèle, faisons un peu de mathématiques.

Les techniques de régression permettent d’avoir des jolis théorèmes, souvent d’une portée incroyablement générale (moyennant quelques petites hypothèses techniques). Par exemple le théorème de Frisch-Waugh, en régression multiple, dont j’ai déjà parlé dans des vieux billets. Un des corolaires est que lorsque les variables explicatives dans un modèle de régression sont orthogonales, la régression multiple correspond à une collection de régressions simples (autrement dit, les estimateurs par moindres carrés coïncident). Formellement, si on considère le modèley_i=\beta_0+\beta_1x_{1,i}+\beta_2x_{2,i}+\varepsilon_i(avec les hypothèses usuelles des modèles de régression) alors, si les variables x_1 et x_2 sont non-corrélées, \widehat{\beta}_1 coïncide avec \widehat{b}_1 dans le modèley_i=b_0+b_1x_{1,i}+\eta_iOn peut faire une petite simulation pour confirmer (sinon, bien entendu, on peut regarder la démonstration qui se trouve dans tous les livres d’économétrie, qui est d’ailleurs un simple résultat d’algèbre linéaire, ou de géométrie, avec des projections successives sur des sous-espaces orthogonaux – même si c’est à Michael Lovell que l’on doit l’approche géométrique).

library(mnormt)
r = 0
S = matrix(c(1,r,r,1),2,2)
n = 1000
set.seed(1)
X = rmnorm(n,c(0,0),S)
E = rnorm(n,0,.3)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))

Petite note technique: je vais centrer les variables, histoire de ne pas avoir à garder la constante dans mon modèle (qui va compliquer les notations, et potentiellement embrouiller un peu le billet). La constante est nulle ici, on le voit,

reg = lm(Y~X1+X2,data=base)
summary(reg)
 
Coefficients:
              Estimate Std. Error t value  Pr(|t|)    
(Intercept)  7.449e-18  9.777e-03     0.0        1    
X1           1.012e+00  9.171e-03   110.3    2e-16 ***
X2          -1.988e+00  9.719e-03  -204.6    2e-16 ***

Bref, je régresse sans constante

reg = lm(Y~0+X1+X2,data=base)
summary(reg)
 
Coefficients:
    Estimate Std. Error t value  Pr(|t|)    
X1  1.011520   0.009166   110.3   2e-16 ***
X2 -1.988321   0.009714  -204.7   2e-16 ***

Maintenant, on va regarder les deux régressions simple

reg1 = lm(Y~0+X1,data=base)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1  1.01300    0.06006   16.86   2e-16 ***

quand on régresse juste sur la première variable, et pour la seconde, on obtient

reg2 = lm(Y~0+X2,data=base)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X2 -1.98916    0.03528  -56.39   2e-16 ***

Autrement dit, nos estimateurs sont très proches (si on avait laissé la constante, ils coïncideraient).

Maintenant, le gros soucis est que ce résultat n’est plus valide lorsque les variables explicatives sont corrélées. Le théorème de Frisch-Waugh nous explique comment ces estimateurs divergent, mais le point ici est que si les variables sont corrélées, utiliser les régressions simples donne deux estimateurs biaisés des vrais paramètres (du modèle multiple). Recommençons l’exercice précédant

r=.9
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~X1+X2,data=base)
summary(reg)
reg = lm(Y~0+X1+X2,data=base)
summary(reg)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1  0.98740    0.02205   44.79   2e-16 ***
X2 -1.97321    0.02229  -88.54   2e-16 ***

(on retrouve des valeurs proches de celles utilisées pour simuler nos données, donc tout va bien). En revanche, pour les régressions simples, on obtient

reg1 = lm(Y~0+X1,data=base)
summary(reg1)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X1 -0.78784    0.02726   -28.9   2e-16 ***

et

reg2 = lm(Y~0+X2,data=base)
summary(reg2)
 
Coefficients:
   Estimate Std. Error t value Pr(|t|)    
X2 -1.06543    0.01607  -66.31   2e-16 ***

Autrement dit, \widehat{b}_1\neq \widehat{\beta}_1 (et pareil pour le second). Ce qui signifie que si on construit une prévision à partir de ce modèle, \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, on sera potentiellement très loin de la “bonne” prévision \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i} (qui sera sans biais, etc, je renvoie ici vers n’importe quel cours d’économétrie). On peut le voir sur un dessin,

Yp=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
plot(base$Y,predict(reg),ylim=range(Yp),col=rgb(0,0,1,.5),cex=.7,xlab="observé",ylab="prédit")
abline(a=0,b=1,lty=2)
points(base$Y,Yp,col=rgb(1,0,0,.5),cex=.7)
abline(lm(predict(reg)~base$Y),col="blue")
abline(lm(Yp~base$Y),col="red")

En bleu, on a les prévisions avec le modèle linéaire multiple, et en rouge, en faisant deux régressions indépendantes… Si on regarde sur la droite, le modèle rouge sur-valorise, ou disons sur-estime largement, alors qu’il sous-valorise à gauche. La différence entre les deux droites s’interprète ici comme un biais. Sur le graphique ci-dessous, on peut visualiser la distribution des \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, en rouge, et les \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i} , en bleu. Cet excès de dispersion, de variance, qu’on observe sur les points rouges, j’interprète ça comme de la polarisation

plot(density(Yp),col="red")lwd=2)
lines(density(predict(reg)),col="blue",lwd=2)

En fait, ce que raconte le théorème de Frisch-Waugh (et je renvoie à mon précédant billet pour plus de détails), c’est qu’on a le droit de faire plusieurs régressions, mais en cascade ! et surtout pas indépendamment : je peux expliquer y avec la première variable x_1, et ensuite régresser le résidu (ce qu’on n’a pas pu expliquer) sur la seconde x_2. Cette méthode donnera la même prévision que le modèle multiple.

On peut aller un peu plus loin: on peut jouer sur la valeur de la corrélation, pour mesure l’écart entre les deux prévisions (ici je prévois pour une observation au hasard, la 78ème)

comp=function(r){
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~0+X1+X2,data=base)
reg1 = lm(Y~0+X1,data=base)
reg2 = lm(Y~0+X2,data=base)
y1=predict(reg)
y2=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
c(y1[78],y2[78],(y2[78]-y1[78])/y1[78])}
vR=seq(0,.98,by=.02)
vc=Vectorize(comp)(vR)
plot(vR,vc[3,]*100,ylab="Différence relative (%)",xlab="Corrélation",type="l")

On prédit ici pour une observation avec un large y_i, ce qui correspondait à la partie de droite du graphique précédant (avec les points rouges et bleus).

C’est ce que j’appelle le deuxième effet kiss-cool. Quand on est dans le premier cas, avec les variables indépendantes (corrélation nulle, c’est à dire à gauche sur ma figure), on explique ce qu’on peut avec la première variable, et on rajoute l’impact de la seconde. Et \widetilde{y}_i\sim\widehat{y}_i. Le soucis ici est que je n’ai pas le droit de considérer deux modèles indépendants lors que les variables sont très corrélées. Une partie de l’explication fournie par la seconde variable était déjà inclue dans la première. Par exemple avec deux variables très (positivement) corrélées, la prévision qu’on obtient en ajoutant les deux effets estimés indépendamment avec deux régressions simple \widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}, on sur-estime de 50% à 70% la “vraie” prévision \widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i}.

Tous les chercheurs savent savent ça… on parle ici de résultats du tout premier cours de modèles linéaires. Et malgré tout, en pratique, on continue à utiliser cette seconde méthode. Un exemple bien connu est celui de l’évaluation (des étudiants, des chercheurs, peu importe). Par exemple, quand on évalue un dossier de financement pour un chercheur, on nous demande de mettre un score

  • pour les publications scientifiques (nombre, qualité, etc)
  • pour l’encadrement d’étudiants (nombre, niveau, etc)
  • pour la qualité de l’environnement (prestige du labo, etc)
  • etc

Et à la fin, on somme tout. Mais on le voit, ces variables sont très très corrélées: si vous êtes dans un labo prestigieux, vous attirez beaucoup de candidatures d’étudiants (et des bons), et avoir beaucoup d’étudiants va permettre d’avoir plus de publications (si on ajoute son nom comme co-auteur). Bref, on est typiquement dans un modèle à double (voire triple) effet kiss-cool. Quelqu’un dans un bon labo aura un bon score sur le troisième item, mais aussi un bon sur le nombre d’étudiants, et aussi sur les publications. Ajouter ces scores est stupide, car on a une spirale infernale (les bons sont sur-évalués, et les moins bon, sous-évalués), c’est ce que racontait mon premier dessin, avec les points rouges et bleus. C’est un effet clivant de polarisation forte.

Si on voulait faire les choses proprement, ce que dit le le théorème de Frisch-Waugh, c’est que les scores devraient être attribués en corrigeant de la corrélation entre les variables

  • on peut commencer par calculer un score pour la qualité de l’environnement (prestige du labo, etc)
  • à environnement donné, calculer un score pour l’encadrement d’étudiants (nombre, niveau, etc)
  • à environnement donné, et à encadrement d’étudiants donné, calculer un score pour les publications
  • etc

C’est comme la situation que je voyais en France, où on pouvait avoir une variable qui tenait compte d’un prestige du chercheur (par exemple, être chercheur CNRS donnait un bonus) et une autre sur le dossier de publications. Sauf que les deux sont corrélés. Et la plupart des classements d’universités sont construits à partir de scores qui sont loin d’être indépendants.

Bref, tant que l’évaluation se fera en sommant des scores qui sont construits sur des critères souvent très corrélés, on polarise fortement la population.

Est-ce gênant? A priori oui. Car le message que cela envoie est qu’il existe deux classes, les bons, et les mauvais, alors qu’en réalité, le niveau est beaucoup plus homogène qu’il n’y paraît.  Un petit effet positif se retrouve démultiplié par le fait qu’il va se répercuter (positivement) sur plein d’autres variables. C’est mon effet kiss-cool. Mais on pourrait se dire que c’est un problème de distribution des notes finales. Si l’ordre est préservé, on pourrait se dire que ce n’est pas très grave. Malheureusement, ce n’est pas le cas.

Si on quitte un instant le cas de la corrélation très forte, les rangs des prédictions (c’est à dire les rangs des chercheurs une fois donnés les notes \widehat{y}_i ou \widetilde{y}_i ) sont moins corrélés si la corrélation sous-jacente est importante (mais pas trop)

comp=function(r){
S=matrix(c(1,r,r,1),2,2)
set.seed(1)
X=rmnorm(n,c(0,0),S)
Y = 2+X[,1]-2*X[,2]+E
base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))
reg = lm(Y~0+X1+X2,data=base)
reg1 = lm(Y~0+X1,data=base)
reg2 = lm(Y~0+X2,data=base)
y1=predict(reg)
y2=reg1$coefficients[1]*base$X1+reg2$coefficients[1]*base$X2
cor(y1,y2,method="spearman")}
vR=seq(0,.98,by=.02)
vc=Vectorize(comp)(vR)
plot(vR,vc,ylab="Corrélation de rangs",xlab="Corrélation",type="l")

Autrement dit, si les variables x_1 et x_2 sont très peu corrélées, on a les mêmes rangs  (globalement). En revanche, si la corrélation entre les variables x_1 et x_2 augmente, le rang des \widetilde{y}_i est de moins en moins cohérent avec celui entre les \widehat{y}_i  (qui devrait être celui que l’on recherche).

Bref, il serait temps de comprendre enfin sérieusement les conséquences de ce joli papier, publié il y a presque 90 ans

Prédire le nombre de morts, suite

Il y a quelques jours, un tweet de Baptiste Coulmont m’inspirait un billet sur l’estimation de la surmortalité, en France. Ce soir, c’est un tweet d’Alexandre Blanchet, qui parlait de surmortalité, au Québec.

J’étais surpris, parce que la prévision est bien au dessus des années passés contrairement à ce qu’on pouvait observer en France. En repartant des codes d’Alexandre (data_prep.R puis predictions.R) j’ai remis un peu en forme le graphique, pour avoir la vision de long-terme, pour comprendre la tendance linéaire utilisée par Alexandre (qui utilise un modèle à la Buys-Ballot)

Je me disais qu’on pouvait tenter des modèles de lissages exponentiel, à la Holt-Winters,

library(forecast)
x=ts(df.l_training$deces,start = 2010, frequency = 12)
x.hw <- ets(x, model = "AAA")
autoplot(forecast(x.hw,24))

en version additive, ou multiplicative

x.hw2 <- ets(x, model = "MAM")
autoplot(forecast(x.hw2,24))

Je laisse à d’autre le soin de comparer aux nombres de morts observés depuis 2 mois, je reste fasciné pour ma part par la tendance à la hausse de la mortalité, depuis 10 ans… Et j’essayerais de comparer ces trois modèles dès que j’aurais finis mes corrections de copies de la session d’hiver.

Testing for Covid-19 in the U.S.

For almost a month, on a daily basis, we are working with colleagues (Romuald, Chi and Mathieu) on modeling the dynamics of the recent pandemic. I learn of lot of things discussing with them, but we keep struggling with the tests. Paul, in Montréal, helped me a little bit, but I think we will still have to more to get a better understand. To but honest, we stuggle with two very simple questions

  • how many people are tested on a daily basis ?

Recently, I discovered Modelling COVID-19 exit strategies for policy makers in the United Kingdom, which is very close to what we try to do… and in the document two interesting scenarios are discussed, with, for the first one, “1 million ‘reliable’ daily tests are deployed” (in the U.K.) and “5 million ‘useless’ daily tests are deployed”. There are about 65 millions unhabitants in the U.K. so we talk here about 1.5% people tested, on a daily basis, or 7.69% people ! It could make sense, but our question was, at some point, is that realistic ? where are we today with testing ? In the U.S. https://covidtracking.com/ collects interesting data, on a daily basis, per state.

url = "https://raw.githubusercontent.com/COVID19Tracking/covid-tracking-data/master/data/states_daily_4pm_et.csv"
download.file(url,destfile="covid.csv")
base = read.csv("covid.csv")

Unfortunately, there is no information about the population. That we can find on wikipedia. But in that table, the state is given by its full name (and the symbol in the previous dataset). So we new also to match the two datasets properly,

url="https://en.wikipedia.org/wiki/List_of_states_and_territories_of_the_United_States_by_population"
download.file(url,destfile = "popUS.html")
#pas contaminé 2/3 R=3
library(XML)
tables=readHTMLTable("popUS.html")
T=tables[[1]][3:54,c("V3","V4")]
names(T)=c("state","pop")
url="https://en.wikipedia.org/wiki/List_of_U.S._state_abbreviations"
download.file(url,destfile = "nameUS.html")
tables=readHTMLTable("nameUS.html")
T2=tables[[1]][13:63,c(1,4)]
names(T2)=c("state","symbol")
T=merge(T,T2)
T$population = as.numeric(gsub(",", "", T$pop, fixed = TRUE))
names(base)[2]="symbol"
base = merge(base,T[,c("symbol","population")])

Now our dataset is fine… and we can get a function to plot the number of people tested in the U.S. (cumulated). Here, we distinguish between the positive and the negative,

drawing = function(st ="NY"){
sbase=base[base$symbol==st,c("date","positive","negative","population")]
sbase$DATE = as.Date(as.character(sbase$date),"%Y%m%d")
sbase=sbase[order(sbase$DATE),]
par(mfrow=c(1,2))
plot(sbase$DATE,(sbase$positive+sbase$negative)/sbase$population,ylab="Proportion Test (/population of state)",type="l",xlab="",col="blue",lwd=3)
lines(sbase$DATE,sbase$positive/sbase$population,col="red",lwd=2)
legend("topleft",c("negative","positive"),lwd=2,col=c("blue","red"),bty="n")
title(st)
plot(sbase$DATE,sbase$positive/(sbase$positive+sbase$negative),ylab="Ratio of positive tests",ylim=c(0,1),type="l",xlab="",col="black",lwd=3)
title(st)}

Let us start with New York

drawing("NY")

As at now, 4% of the entiere population got tested… over 6 weeks…. The graph on the right is the proportion of people who tested positive… I won’t get back on that one here today, I keep it for our work. In New Jersey, we got about 2.5% of the entiere population tested, overall,

drawing("NJ")

Let us try a last one, Florida

drawing("FL")

As at today, it is 1.5% of the population, over 6 weeks. Overall, in the U.S. less than 0.1% people are tested, on a daily basis. Which is far from the 1.5% in the U.K. scenarios. Now, here come the second question,

  • what are we actually testing for ?

On that one, my experience in biology is… very limited, and Paul helped me. He mentioned this morning a nice report, from a lab in UC Berkeley

One of my question was for instance, if you get tested positive, and you do it again, can you test negative ? Or, in the context of our data, do we test different people ? are some people tested on a regular basis (perhaps every week) ? For instance, with antigen tests (Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) – also called molecular or PCR – Polymerase Chain Reaction – test) we test if someone is infectious, while with antibody test (using serological immunoassays that detect viral-specific antibodies — Immunoglobin M (IgM) and G (IgG) — also called serology test), we test for immunity. Which is rather different…

I have no idea what we have in our database, to be honest… and for the past six weeks, I have seen a lot of databases, and most of the time, I don’t know how to interpret, I don’t know what is measured… and it is scary. So, so far, we try to do some maths, to test dynamics by tuning parameters “the best we can” (and not estimate them). But if anyone has good references on testing, in the context of Covid-19 (for instance on specificity, sensitivity of all those tests) I would love to hear about it !

Estimer la surmortalité

Ce matin, Baptiste Coulmont m’envoyait un tweet avec un joli graphique, présentant le nombre quotidien de décès en France.

Comme les données sont en ligne, je me suis dit que je pourrais jouer un peu avec. Pour les plus curieux, on a la liste de tous les décès depuis… longtemps ! (on a plus de 25 millions de décès). La seule chose qui nous intéresse, c’est la date, alors on la récupère. Puis on compte, par jour, combien de décès on a eu,

D = read.csv("insee_deces.csv",header=TRUE)
Vecteur_Dates = as.character(D[,8])
TV = table(Vecteur_Dates)
B = data.frame(dateC=names(TV),nb=as.numeric(TV))
B$date=as.Date(B$dateC,"%Y-%m-%d")
B$year = format(B$date,"%Y")
B$month = format(B$date,"%m")
B$day = format(B$date,"%d")
B$deb=as.Date(paste(B$year,"01-01",sep="-"),"%Y-%m-%d")
B$fin=as.Date(paste(B$year,"12-31",sep="-"),"%Y-%m-%d")
B$dif=as.numeric(B$date-B$deb)/as.numeric(B$fin-B$deb)

Je crée ici une une variable qui me dit où je suis dans l’année en cours (0 pour le 1er janvier, et 1 pour pour 31 décembre). Ça me permet de contourner un petit problème technique : les années bissextiles. On peut ensuite représenter toutes les années,

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",xlim=0:1,ylim=c(1200,4000),col="white")
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")

On veut ensuite avoir la tendance moyenne. Baptiste prenait la moyenne par jour (la courbe noire). Ici, on va lisser avec des splines

library(gam)
reg = gam(NormCpte~bs(dif,40),data=B2[B2$year != "2020",])
vx = seq(0,1,length=501)
vy = predict(reg,newdata=data.frame(dif=vx))
lines(vx,vy,col="blue",lwd=3)

ce qui donne la courbe bleue. C’est l’écart à cette courbe qui donne la surmortalité. Par exemple (prenons un exemple bien connu, et largement visible sur le graphique précédant) l’année 2003, et les quinze premiers jours d’août, correspondant (plus ou moins) à la période de canicule.

i = 2003
x = B2$dif[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
y = B2$NormCpte[(B2$year == i)&((B2$month == "08"))&(B2$day %in% c(paste("0",1:9,sep=""),10:15))]
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 14294.39

Plusieurs chiffres sont mentionnés sur wikipédia, dont le rapport de l’INSERM qui annonçait 14,802 morts. On n’est pas trop loin… On peut visualiser cette surmortalité sur le graphique ci-dessous

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",col="white",xlim=0:1,ylim=c(1200,4000))
for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")
lines(vx,vy,col="blue")
lines(B2$dif[B2$year == 2003],B2$NormCpte[B2$year == 2003],col="red")
for(u in 1:length(x)) segments(x[u],y[u],x[u],yp[u],lwd=3,col="red")

Prenons un autre exemple qu’on voit bien sur le graphique de Baptiste, la grippe de 2017 (dite de 2016-2017). Si on prend les 45 premiers jours de 2017, on obtient les chiffres suivants

i = 2017
x=B2$dif[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
y=B2$NormCpte[(B2$year == i)&((B2$month == "01")|
  ((B2$month == "02")&(B2$day %in% c(paste("0",1:9,sep=""),10:15))))]
yp=predict(reg,newdata=data.frame(dif=x))
yp = predict(reg,newdata=data.frame(dif=x))
e = y-yp
sum(e)
[1] 21177.33

Le rapport de Santé Publique France parle de 21,200 décès (pour la “surmortalité toute cause”) ce qui est là aussi comparable…

Maintenant, avant de conclure, et avant de me faire troller pendant quelques jours, je tiens à préciser que c’est juste un exercice de modélisation, et oui, j’ai exclu 2020. Je ne tiens pas à faire de comparaison avec la Covid-19, à qui 21,000 morts sont attribués, en France, ce soir. Si l’ordre de grandeur est le même, qu’on en me fasse pas dire qu’on est dans une situation comparable, et que la Covid-19 est finalement une simple grippe. Je ne me prononcerais pas sur ce point car je ne suis pas virologue. Mais surtout, on a 21,000 morts malgré des mesures incroyables de confinement. Je reviendrais probablement dans les jours qui viennent sur l’impact du contrôle par confinement, mais l’objet du court billet de ce soir était de voir comment on pouvait rapidement quantifier une surmortalité.

  • petite note technique : pour les 3 graphiques (mais pas pour les codes et les calculs) j’avais initialement normalisé les données. En effet, la population française a relativement pas mal changé en 20 ans, augmentant de 10%… Aussi, 1000 morts en 2003, ça n’est pas autant (relativement à la taille de la population) qu’en 2020.

Modeling Pandemics (3)

In Statistical Inference in a Stochastic Epidemic SEIR Model with Control Intervention, a more complex model than the one we’ve seen yesterday was considered (and is called the SEIR model). Consider a population of size N, and assume that S is the number of susceptible, E the number of exposed, I the number of infectious, and R for the number recovered (or immune) individuals, \displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-\beta {\frac {I}{N}}S\\[8pt]{\frac {dE}{dt}}&=\beta {\frac {I}{N}}S-aE\\[8pt]{\frac {dI}{dt}}&=aE-b I\\[8pt]{\frac {dR}{dt}}&=b I\end{aligned}}Between S and I, the transition rate is \beta I, where \beta is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject. Between I and R, the transition rate is b (simply the rate of recovered or dead, that is, number of recovered or dead during a period of time divided by the total number of infected on that same period of time). And finally, the incubation period is a random variable with exponential distribution with parameter a, so that the average incubation period is a^{-1}.

Probably more interesting, Understanding the dynamics of ebola epidemics suggested a more complex model, with susceptible people S, exposed E, Infectious, but either in community I, or in hospitals H, some people who died F and finally those who either recover or are buried and therefore are no longer susceptible R.

Thus, the following dynamic model is considered\displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}\\[8pt]\frac {dE}{dt}&=(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}-\alpha E\\[8pt]\frac {dI}{dt}&=\alpha E+\theta\gamma_H I-(1-\theta)(1-\delta)\gamma_RI-(1-\theta)\delta\gamma_FI\\[8pt]\frac {dH}{dt}&=\theta\gamma_HI-\delta\lambda_FH-(1-\delta)\lambda_RH\\[8pt]\frac {dF}{dt}&=(1-\theta)(1-\delta)\gamma_RI+\delta\lambda_FH-\nu F\\[8pt]\frac {dR}{dt}&=(1-\theta)(1-\delta)\gamma_RI+(1-\delta)\lambda_FH+\nu F\end{aligned}}In that model, parameters are \alpha^{-1} is the (average) incubation period (7 days), \gamma_H^{-1} the onset to hospitalization (5 days), \gamma_F^{-1} the onset to death (9 days), \gamma_R^{-1} the onset to “recovery” (10 days), \lambda_F^{-1} the hospitalisation to death (4 days) while \lambda_R^{-1} is the hospitalisation to recovery (5 days), \eta^{-1} is the death to burial (2 days). Here, numbers are from Understanding the dynamics of ebola epidemics (in the context of ebola). The other parameters are \beta_I the transmission rate in community (0.588), \beta_H the transmission rate in hospital (0.794) and \beta_F the transmission rate at funeral (7.653). Thus

epsilon = 0.001 
Z = c(S = 1-epsilon, E = epsilon, I=0,H=0,F=0,R=0)
p=c(alpha=1/7*7, theta=0.81, delta=0.81, betai=0.588,
    betah=0.794, blambdaf=7.653,N=1, gammah=1/5*7,
    gammaf=1/9.6*7, gammar=1/10*7, lambdaf=1/4.6*7,
    lambdar=1/5*7, nu=1/2*7)

If \boldsymbol{Z}=(S,E,I,H,F,R), if we write \frac{\partial \boldsymbol{Z}}{\partial t} = SEIHFR(\boldsymbol{Z})where SEIHFR is

SEIHFR = function(t,Z,p){
  S=Z[1]; E=Z[2]; I=Z[3]; H=Z[4]; F=Z[5]; R=Z[6]
  alpha=p["alpha"]; theta=p["theta"]; delta=p["delta"]
  betai=p["betai"]; betah=p["betah"]; gammah=p["gammah"]
  gammaf=p["gammaf"]; gammar=p["gammar"]; lambdaf=p["lambdaf"]
  lambdar=p["lambdar"]; nu=p["nu"]; blambdaf=p["blambdaf"]
  N=S+E+I+H+F+R
  dS=-(betai*I+betah*H+blambdaf*F)*S/N
  dE=(betai*I+betah*H+blambdaf*F)*S/N-alpha*E
  dI=alpha*E-theta*gammah*I-(1-theta)*(1-delta)*gammar*I-(1-theta)*delta*gammaf*I
  dH=theta*gammah*I-delta*lambdaf*H-(1-delta)*lambdaf*H
  dF=(1-theta)*(1-delta)*gammar*I+delta*lambdaf*H-nu*F
  dR=(1-theta)*(1-delta)*gammar*I+(1-delta)*lambdar*H+nu*F
  dZ=c(dS,dE,dI,dH,dF,dR)
  list(dZ)}

We can solve it, or at least study the dynamics from some starting values

library(deSolve)
times = seq(0, 50, by = .1)
resol = ode(y=Z, times=times, func=SEIHFR, parms=p)

For instance, the proportion of people infected is the following

plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="",col="red")
lines(resol[,"time"],resol[,"H"],col="blue")

Modeling pandemics (2)

When introducing the SIR model, in our initial post, we got an ordinary differential equation, but we did not really discuss stability, and periodicity. It has to do with the Jacobian matrix of the system. But first of all, we had three equations for three function, but actually\displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}so it means that our problem is here simply in dimension 2. Hence\displaystyle {\begin{aligned}&X={\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&Y={\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\mu+\gamma)I\end{aligned}}and therefore, the Jacobian of the system is\begin{pmatrix}\displaystyle{\frac{\partial X}{\partial S}}&\displaystyle{\frac{\partial X}{\partial I}}\\[9pt]\displaystyle{\frac{\partial Y}{\partial S}}&\displaystyle{\frac{\partial Y}{\partial I}}\end{pmatrix}=\begin{pmatrix}\displaystyle{-\mu-\beta\frac{I}{N}}&\displaystyle{-\beta\frac{S}{N}}\\[9pt]\displaystyle{\beta\frac{I}{N}}&\displaystyle{\beta\frac{S}{N}-(\mu+\gamma)}\end{pmatrix}We should evaluate the Jacobian at the equilibrium, i.e. S^\star=\frac{\gamma+\mu}{\beta}=\frac{1}{R_0}andI^\star=\frac{\mu(R_0-1)}{\beta}We should then look at eigenvalues of the matrix.

Our very last example was

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")

We can compute values at the equilibrium

mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
N=1
S = (gamma + mu)/beta
I = mu * (beta/(gamma + mu) - 1)/beta

and the Jacobian matrix

J=matrix(c(-(mu + beta * I/N),-(beta * S/N),
         beta * I/N,beta * S/N - (mu + gamma)),2,2,byrow = TRUE)

Now, if we look at the eigenvalues,

eigen(J)$values
[1] -0.024975+0.6318831i -0.024975-0.6318831i

or more precisely 2\pi/b where a\pm ib are the conjuguate eigenvalues

2 * pi/(Im(eigen(J)$values[1]))
[1] 9.943588

we have a damping period of 10 time lengths (10 days, or 10 weeks), which is more or less what we’ve seen above,

The graph above was obtained using

p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[1:1e5,"time"],resol[1:1e5,"I"],type="l",xlab="time",ylab="",lwd=3,col="red")
yi=resol[,"I"]
dyi=diff(yi)
i=which((dyi[2:length(dyi)]*dyi[1:(length(dyi)-1)])<0)
t=resol[i,"time"]
arrows(t[2],.008,t[4],.008,length=.1,code=3)

If we look carefully. at the begining, the duration is (much) longer than 10 (about 13)… but it does converge towards 9.94

plot(diff(t[seq(2,40,by=2)]),type="b")
abline(h=2 * pi/(Im(eigen(J)$values[1]))

So here, theoretically, every 10 weeks (assuming that our time length is a week), we should observe an outbreak, smaller than the previous one. In practice, initially it is every 13 or 12 weeks, but the time to wait between outbreaks decreases (until it reaches 10 weeks).

Modeling pandemics (1)

The most popular model to model epidemics is the so-called SIR model – or Kermack-McKendrick. Consider a population of size N, and assume that S is the number of susceptible, I the number of infectious, and R for the number recovered (or immune) individuals, \displaystyle {\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}so that \displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}which implies that S+I+R=N. In order to be more realistic, consider some (constant) birth rate \mu, so that the model becomes\displaystyle {\begin{aligned}&{\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\gamma+\mu) I,\\[6pt]&{\frac {dR}{dt}}=\gamma I-\mu R,\end{aligned}}Note, in this model, that people get sick (infected) but they do not die, they recover. So here, we can model chickenpox, for instance, not SARS.

The dynamics of the infectious class depends on the following ratio:\displaystyle{R_{0}={\frac {\beta }{\gamma +\mu}}} which is the so-called basic reproduction number (or reproductive ratio). The effective reproductive ratio is R_0S/N, and the turnover of the epidemic happens exactly when R_0S/N=1, or when the fraction of remaining susceptibles is R_0^{-1}. As shown in Directly transmitted infectious diseases:Control by vaccination, if S/N<R_0^{-1} the disease (the number of people infected) will start to decrease.

Want to see it  ? Start with

mu = 0
beta = 2
gamma = 1/2

for the parameters. Here,  R_0=4. We also need starting values

epsilon = .001
N = 1
S = 1-epsilon
I = epsilon
R = 0

Then use the ordinary differential equation solver, in R. The idea is to say that \boldsymbol{Z}=(S,I,R) and we have the gradient \frac{\partial \boldsymbol{Z}}{\partial t} = SIR(\boldsymbol{Z})where SIR is function of the various parameters. Hence, set

p = c(mu = 0, N = 1, beta = 2, gamma = 1/2)
start_SIR = c(S = 1-epsilon, I = epsilon, R = 0)

The we must define the time, and the function that returns the gradient,

times = seq(0, 10, by = .1)
SIR = function(t,Z,p){
S=Z[1]; I=Z[2]; R=Z[3]; N=S+I+R
mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
dS=mu*(N-S)-beta*S*I/N
dI=beta*S*I/N-(mu+gamma)*I
dR=gamma*I-mu*R
dZ=c(dS,dI,dR)
return(list(dZ))}

To solve this problem use

library(deSolve)
resol = ode(y=start_SIR, times=times, func=SIR, parms=p)

We can visualize the dynamics below

par(mfrow=c(1,2))
t=resol[,"time"]
plot(t,resol[,"S"],type="l",xlab="time",ylab="")
lines(t,resol[,"I"],col="red")
lines(t,resol[,"R"],col="blue")
plot(t,t*0+1,type="l",xlab="time",ylab="",ylim=0:1)
polygon(c(t,rev(t)),c(resol[,"R"],rep(0,nrow(resol))),col="blue")
polygon(c(t,rev(t)),c(resol[,"R"]+resol[,"I"],rev(resol[,"R"])),col="red")

We can actually also visualize the effective reproductive number is R_0S/N, where

R0=p["beta"]/(p["gamma"]+p["mu"])

The effective reproductive number is on the left, and as we mentioned above, when we reach 1, we actually reach the maximum of the infected,

plot(t,resol[,"S"]*R0,type="l",xlab="time",ylab="")
abline(h=1,lty=2,col="red")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),1,pch=19)
plot(t,resol[,"S"],type="l",xlab="time",ylab="",col="grey")
lines(t,resol[,"I"],col="red",lwd=3)
lines(t,resol[,"R"],col="light blue")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),max(resol[,"I"]),pch=19)

And when adding a \mu parameter, we can obtain some interesting dynamics on the number of infected,

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")