# Model selection, AIC and Tweedie regression

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.

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

Dans Confessiones, Saint Augustin écrivait

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

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.

# Interpretability and explainability of predictive models

In 400 AD, in his Confessiones, Augustine wrote

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

that can be translated as

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

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

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

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

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

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

# GLM, STT5100

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,

1. introduction générale video pdf (10:15)
2. lois de Bernoulli, binomiale, multnomiale video pdf (29:23)
3. régression logistique (Bernoulli) video pdf (23:04)
4. régression multinomiale video pdf (21:45)
5. régression logistique sur variables catégorielles video pdf (30:24)
6. régression logistique sur variables continues video pdf (21:35)
7. analyse discriminante et courbe ROC video pdf (56:53)
8. modèles de comptage et loi de Poisson video pdf (19:28)
9. régression de Poisson video pdf (25:38)
10. régression de Poisson et interprétations video pdf (40:15)
11. régression de Poisson et méthode des marges video pdf (25:29)
12. régression de Poisson et application en assurance video pdf (36:23)
13. famille exponentielle video pdf (30:01)
14. famille exponentielle et GLM video pdf (41:02)
15. loi et lien video video pdf (36:05)
16. déviance et résidus video pdf (15:10)
17. modèle Tweedie et poids video pdf (29:33)
18. surdispersion video pdf (21:01)
19. tests et GLM video pdf (19:41)
20. GLM en petite dimension video pdf (23:34)
21. méthode stepwise video pdf (17:57)
22. Poisson vs. Binomiale, application en démographie video pdf (22:25)
23. Exemple (1) video + pdf

# STT5100, quiz (régression de Poisson #1)

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(&gt;|z|) (Intercept) 6.8844970 0.0110463 623.241 &lt;2e-16 *** HIGH_T 0.0210950 0.0003133 67.328 &lt;2e-16 *** LOW_T -0.0114006 0.0003351 -34.024 &lt;2e-16 *** PRECIP -0.6570450 0.0071899 -91.384 &lt;2e-16 *** I(PRECIP == 0)TRUE 0.1303908 0.0033283 39.176 &lt;2e-16 *** DAYJeu 0.1683475 0.0049690 33.880 &lt;2e-16 *** DAYLun 0.1144129 0.0050480 22.665 &lt;2e-16 *** DAYMar 0.1655886 0.0049936 33.160 &lt;2e-16 *** DAYMer 0.1688035 0.0049190 34.317 &lt;2e-16 *** DAYSam 0.0466447 0.0051838 8.998 &lt;2e-16 *** DAYVen 0.1003536 0.0050919 19.709 &lt;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.

# Slides 21 – Poisson vs. Binomiale

Poisson ou Binomiale ? Application sur un problème de démographie. Les slides sont en ligne (slides 21) et la vidéo aussi (slides 21)

# Slides 20 – GLM et sélection de variables (stepwise)

Rapide discussion sur l’utilisation des méthodes de simulations en petite dimension dans les GLM. Les slides sont en ligne (slides 20) et la vidéo aussi (slides 20)

# Slides 13 – famille exponentielle et GLM

Après avoir présenté la famille exponentielle, on va voir comme l’utiliser dans un contexte de régression. Les slides sont en ligne (slides 13) et la vidéo aussi (slides 13)

# Slides 11 – Modèle de comptage, application

Pour finir avec la régression de Poisson, une application sur des données d’assurance automobile. Les slides sont en ligne (slides 11) et la vidéo aussi (slides 11)

# Slides 9 – régression de Poisson et interprétations

On va revenir un peu sur les interprétations de la régression de Poisson. Les slides sont en ligne (slides 9) et la vidéo aussi (slides 9)

# Slides 7 – modèle de Poisson

On va maintenant parler des modèles de comptage, en revenant sur le modèle de Poisson. Les slides sont en ligne (slides 7) et la vidéo aussi (slides 7)

# Slides 6 – analyse discriminante et courbe ROC

Dernière série de slides sur la régression logistique, avec une présentation de l’analyse discriminante, et la présentation de la courbe ROC. Les slides sont en ligne (slides 6) et la vidéo aussi (slides 6)

# Slides 5 – régression logistique sur variable(s) continue(s)

Après quelques slides sur la régression d’une variable binaire sur des variables catégorielles, on va revenir un peu sur la régression sur une variable continue. Les slides sont en ligne (slides 5) et la vidéo aussi (slides 5)

# Slides 4 – régression logistique sur variable(s) catégorielle(s)

Après quelques slides sur la loi multinomiale, on va revenir un peu sur l’interprétation, en commençant par la régression sur une, ou plusieurs, variable catégorielle. Les slides sont en ligne (slides 4) et la vidéo aussi (slides 4)

# Slides 3 – régression multinomiale

Après quelques slides sur la loi logistique, on va revenir un peu sur la régression multinomiale. Les slides sont en ligne (slides 3) et la vidéo aussi (slides 3)