Tomorrow morning, I will give a lecture online on Balance, Calibration and Related Concepts, for actuaries, working in an insurance company. We will start with calibration


We will also discuss Lorenz and Concentration curves

Tomorrow morning, I will give a lecture online on Balance, Calibration and Related Concepts, for actuaries, working in an insurance company. We will start with calibration


We will also discuss Lorenz and Concentration curves

Just some simple codes to illustrate some points we will discuss this week, for the last course on GLMs, before the final exam. We have mentioned that the Gamma distribution belongs to the exponential, so we can run a regression, and compute the associated AIC,
> set.seed(123) > test.data = rgamma(n=2000, scale=1, shape=1) > m1 = glm( test.data~1, family=Gamma(link=log)) > AIC(m1) [1] 3997.332
The Gamma distribution is also a special case of the Tweedie distribution, with power 2
> library(statmod) > library(tweedie) > m2 = glm( test.data~1, family=tweedie(link.power=0, var.power=2) ) > AIC(m2) [1] NA
Unfortunately, we cannot compute the AIC, and we need a trick (with the appropriate R function).
> AICtweedie(m2) [1] 3997.332
Of course, we can do the same with the Poisson distribution, which also belongs to the exponential family
> test.data = rpois(n=2000, lambda=1) > m3 = glm( test.data~1, family=poisson(link=log)) > m4 = glm( test.data~1, family=tweedie(link.power=0, var.power=1) ) > AIC(m3) [1] 5124.61
Here, we have a problem with the AICtweedie function
> AICtweedie(m4) [1] Inf
because we need to specify the dispersion parameter
> AICtweedie(m4, dispersion=1) [1] 5124.61
We can now check: we generate some Gamma sample, and fit various Tweedie distribution, changing simply the variance function (which is a power function)
> set.seed(123)
> test.data = rgamma(n=2000, scale=1, shape=1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+
+ }
> vt = seq(1,2.7,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")
The minimum of the AIC is close to 2, corresponding to the Gamma distribution

We can also try with a Poisson
> set.seed(123)
> test.data = rpois(n=2000, lambda=1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+
+ }
> vt = seq(1,2,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")
The minimum is now close to 1, corresponding to the Poisson distriubtion (the variance is equal to the average)

Let us now try some compound Poisson distribution,
> rcpd=function(n,lambda,shape,scale){
+ N=rpois(n,lambda)
+ X=rgamma(sum(N),shape=shape, scale=scale)
+ I=as.factor(rep(1:n,N))
+ S=tapply(X,I,sum)
+ V=as.numeric(S[as.character(1:n)])
+ V[is.na(V)]=0
+ return(V)}
Let us generate some compound Poisson random variables, with Poisson distribution with average 1, and with Gamma jumps, with average and variance 1,
> set.seed(123)
> test.data = rcpd(n=2000, 1,1,1)
> glmtw = function(t){
+ m1 = glm( test.data~1, family=tweedie(link.power=0, var.power=t) )
+ d = NULL
+ if(t == 1) d = 1
+ AICtweedie(m1, dispersion = d)
+ }
> vt = seq(1.1,1.9,length=100)
> vg = Vectorize(glmtw)(vt)
> plot(vt,vg,log="y",type="l")
The optimal value for the power function is here 1.5, based on the AIC (relationships between Tweedie parameters and the compound Poisson ones are given in the slides)

We can now play a little bit with the variance of the jumps: they still have aveage 1, but they now have a smaller variance
> set.seed(123) > test.data = rcpd(n=2000, 1,3,1/3) > vt = seq(1.05,1.95,length=100) > vg = Vectorize(glmtw)(vt) > plot(vt,vg,log="y",type="l")
The optimal power for the Tweedie is closer to one, closer to the Poison case

while if we increase the variance of the jumps
> set.seed(123) > test.data = rcpd(n=2000, 1,1/3,3) > vt = seq(1.05,1.95,length=100) > vg = Vectorize(glmtw)(vt) > plot(vt,vg,log="y",type="l")
the optimal power is higher, closer to the Gamma distribution.

Dans Confessiones, Saint Augustin écrivait
quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio
que l’on traduit
Qu’est-ce-que le temps ? Si personne ne me le demande, je le sais. Si je veux l’expliquer à qui me le demande, je ne le sais plus.
Pour aller un peu plus loin (car souvent, si on nous demande d’expliquer, on a quelques idées), dans Une étude en rouge de Sir Arthur Conan Doyle, paru en 1887, on a l’échange suivant, entre Sherlock Holmes et le docteur Watson,
– Je me demande ce que cherche ce type là-bas, demandai-je,
désignant un grand individu habillé simplement qui suivait
l’autre côté de la rue, en examinant anxieusement les numéros.
Il tenait à la main une grande enveloppe bleue et, de toute
évidence, portait un message.
– Vous parlez de ce sergent d’infanterie de marine ? dit Sherlock Holmes.
puis, comme il s’avère que la personne est effectivement sergent dans la marine (tout comme un autre personnage de l’histoire, un certain Arthur Charpentier), le docteur Holmes lui demande une explication, il veut savoir comment il est arrivé à cette conclusion
« Comment diable avez-vous pu deviner cela ? demandai-je.
– Deviner quoi ? fit-il sans aménité.
– Eh bien, qu’il était un sergent de marine en retraite ?
– Je n’ai pas de temps à perdre en bagatelles ! répondit-il
avec brusquerie avant d’ajouter dans un sourire : excusez ma rudesse ! Vous avez rompu le fil de mes pensées. Mais c’est peut-être aussi bien. Ainsi donc vous ne voyiez pas que cet homme était un sergent de marine ?
– Non, certainement pas !
– Décidément, l’explication de ma méthode me coûte plus que son application ! Si l’on vous demandait de prouver que deux et deux font quatre, vous seriez peut-être embarrassé ; et cependant, vous êtes sûr qu’il en est ainsi. Malgré la largeur de la rue, j’avais pu voir une grosse ancre bleue tatouée sur le dos de la main du gaillard. Cela sentait la mer. Il avait la démarche militaire et les favoris réglementaires ; c’était, à n’en pas douter, un marin. Il avait un certain air de commandement et d’importance. Rappelez-vous son port de tête et le balancement de sa canne ! En outre, son visage annonçait un homme d’âge moyen, sérieux, respectable. Tous ces détails m’ont amené à penser qu’il était sergent.
– C’est merveilleux ! m’écriai-je.
– Peuh ! L’enfance de l’art ! dit Holmes, mais d’un air qui me
parut trahir sa satisfaction devant ma surprise et mon admiration manifestes.
(pour être honnête, c’est Liu Cixin qui en parle dans Le problème à trois corps). Pour l’anecdote, c’est la première histoire du couple Holmes-Watson, qui introduit la méthode de travail de Sherlock Holmes. Pour ceux qui sont familier avec les nouvelles, cette approche narrative sera largement reprise par la suite: Sherlock Holmes énonce un fait, le docteur Watson est étonné et demande une explication, et Sherlock Holmes explique, point par point, comment il est arrivé à cette conclusion. C’est un peu cette approche qu’on tente de mettre en place quand on va construire un modèle prédictif : sur la base des données du Titanic, si on prédit que telle personne va mourir, et que telle autre va survivre, on veut comprendre pourquoi le modèle arrive à cette conclusion.
Continue reading Interprétabilité et explicabilité (formalisé) des modèles prédictifs
In 400 AD, in his Confessiones, Augustine wrote
quid est ergo tempus? si nemo ex me quaerat, scio; si quaerenti explicare velim, nescio
that can be translated as
What then is time? If no one asks me, I know what it is. If I wish to explain it to him who asks, I do not know.
To go a little further (because often, if we are asked to explain, we have some ideas), in A Study in Scarlet by Sir Arthur Conan Doyle, published in 1887, we have the following exchange, between Sherlock Holmes and Doctor Watson
– “I wonder what that fellow is looking for?” I asked, pointing to a stalwart, plainly-dressed individual who was walking slowly down the other side of the street, looking anxiously at the numbers. He had a large blue envelope in his hand, and was evidently the bearer of a message.
– “You mean the retired sergeant of Marines,” said Sherlock Holmes.
then, as it turns out that the person is indeed a sergeant in the navy (as is another character in the story, someone named Arthur Charpentier), Dr. Holmes asks him for an explanation, he wants to know how he arrived at this conclusion
– “How in the world did you deduce that?” I asked.
– “Deduce what?” said he, petulantly.
– “Why, that he was a retired sergeant of Marines.”
– “I have no time for trifles,” he answered, brusquely; then with a smile, “Excuse my rudeness. You broke the thread of my thoughts; but perhaps it is as well. So you actually were not able to see that that man was a sergeant of Marines?”
– “No, indeed.”
– “It was easier to know it than to explain why I knew it. If you were asked to prove that two and two made four, you might find some difficulty, and yet you are quite sure of the fact. Even across the street I could see a great blue anchor tattooed on the back of the fellow’s hand. That smacked of the sea. He had a military carriage, however, and regulation side whiskers. There we have the marine. He was a man with some amount of self-importance and a certain air of command. You must have observed the way in which he held his head and swung his cane. A steady, respectable, middle-aged man, too, on the face of him – all facts which led me to believe that he had been a sergeant.”
(to be honest, it is Liu Cixin who talks about it in The Three-Body Problem). For the record, this is the first story of the Holmes-Watson couple, which introduces Sherlock Holmes’ working method. For those who are familiar with the short stories, this narrative approach will be widely used thereafter: Sherlock Holmes states a fact, Dr. Watson is astonished and asks for an explanation, and Sherlock Holmes explains, point by point, how he arrived at this conclusion. This is a bit like the approach we try to implement when we build a predictive model: on the basis of the Titanic data, if we predict that such and such a person will die, and that such and such a person will survive, we want to understand why the model arrives at this conclusion.
Continue reading Interpretability and explainability of predictive models
Dernière ligne droite dans le cours STT5100 de modèles linéaires appliqués. Les supports de cours sont en ligne sur https://github.com/freakonometrics/STT5100. Cette session étant en distanciel, le cours est asynchrone, et je poste régulièrement des capsules vidéos. Les capsules en lien avec les modèles linéaires généralisés (GLM) sont maintenant en ligne,
Pour la fin de la semaine, j’avais donné un petit quiz, basé sur la base suivante, qui donne le nombre de cyclistes à une intersection, à New York City,
download.file("http://freakonometrics.free.fr/NYCVelo.RData","velo.RData") load("velo.RData") str(base) 'data.frame': 214 obs. of 7 variables: $ Date : chr "1-Apr-17" "2-Apr-17" "3-Apr-17" "4-Apr-17" ... $ HIGH_T : num 46 62.1 63 51.1 63 48.9 48 55.9 66 73.9 ... $ LOW_T : num 37 41 50 46 46 41 43 39.9 45 55 ... $ PRECIP : num 0 0 0.03 1.18 0 0.73 0.01 0 0 0 ... $ BB_COUNT: int 606 2021 2470 723 2807 461 1222 1674 2375 3324 ... $ DAY : chr "Sam" "Dim" "Lun" "Mar" ... $ DIFF_T : num 9 21.1 13 5.1 17 7.9 5 16 21 18.9 ...... |
En utiisant une régression de Poisson, il fallait prédire combien de cyclistes passeront un dimanche, s’il fait une température maximale de 85F, minimale de 70F, et s’il ne pleut pas. Et il fallait voir ce que donnerait la prévision pour un lundi.
newbase = data.frame(DAY=as.factor(c("Lun","Dim")), HIGH_T=c(85,85),LOW_T=c(70,70), PRECIP=c(0,0)) |
Faisons un modèle avec toutes les variables explicatives.
reg0 = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+DAY,data=base,family=poisson) |
Rajoutons aussi une indicatrice pour indiquer les jours où il ne pleut pas,
reg = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+I(PRECIP==0)+DAY,data=base,family=poisson) summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 6.8844970 0.0110463 623.241 <2e-16 *** HIGH_T 0.0210950 0.0003133 67.328 <2e-16 *** LOW_T -0.0114006 0.0003351 -34.024 <2e-16 *** PRECIP -0.6570450 0.0071899 -91.384 <2e-16 *** I(PRECIP == 0)TRUE 0.1303908 0.0033283 39.176 <2e-16 *** DAYJeu 0.1683475 0.0049690 33.880 <2e-16 *** DAYLun 0.1144129 0.0050480 22.665 <2e-16 *** DAYMar 0.1655886 0.0049936 33.160 <2e-16 *** DAYMer 0.1688035 0.0049190 34.317 <2e-16 *** DAYSam 0.0466447 0.0051838 8.998 <2e-16 *** DAYVen 0.1003536 0.0050919 19.709 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 70021 on 213 degrees of freedom Residual deviance: 26493 on 203 degrees of freedom AIC: 28580 Number of Fisher Scoring iterations: 4 |
Tout semble vraiment significatif. Et même avoir du sens (si on regarde les signes, par exemple). Et si on a peur de rater un effet non-linéaire, on peut mettre des splines sur toutes les variables continues
library(gam) reggam = gam(BB_COUNT~bs(HIGH_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson) plot(reggam, se=TRUE) |

pour la température maximale, ou pour la température minimale

et la courbe suivante pour la précipitation, avec un lissage linéaire entre l’observation maximale (3) et celle juste avant (vers 1.8)

On peut aussi régresser sur la température minimale, et l’écart entre le maximum et le minimum (dans un modèle linéaire les modèles sont équivalents, mais avec de transformation nonlinéaire, la différence pourrait donner un modèle plus simple)
library(gam) reggam2 = gam(BB_COUNT~bs(DIFF_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson) plot(reggam2, se=TRUE) |

On peut maintenant comparer les quatre modèles, et leurs prévision. Par exemple, pour le modèle linéaire (avec la variable indicatrice indiquant qu’il ne pleut pas),
P = predict(reg,newdata=newbase,type="response",se.fit=TRUE) |
pour le lundi, on obtient l’intervalle de confiance à 95% pour \widehat{\lambda}
P$fit[1]+c(-2,2)*P$se.fit[1] [1] 3349.842 3401.395 |
et pour le dimanche l’intervalle de confiance à 95% est
P$fit[2]+c(-2,2)*P$se.fit[2] [1] 2987.497 3033.861 |
On peut visualiser les quatre intervalles de confiance (avec, de haut en bas, le second modèle gam, le premier, le modèle linéaire avec l’indicatrice d’absence de pluie, puis le premier modèle linéaire)

alors que pour le dimanche, on a

Autrement dit, en changeant de modèle, on change beaucoup les intervalles de confiance sur la prévision (et les intervalles sont parfois complètement disjoints). Ce qui n’est pas forcément une bonne nouvelle.