Category Archives: STT5100

Slides 1 – loi binomiale

Allez, on reprend les cours, après un peu plus de deux semaines d’interruption… Mais je reprends tranquillement, en revenant sur des rappels que j’avais fait au premier cours, et dont on aura finalement besoin très bientôt ! On revient sur la loi de Bernoulli, la loi binomiale et la loi multinomiale. Et on parlera inférence, et tests. Les slides sont en ligne (slides 1) et la vidéo aussi (slides 1)

Il y a une coquille dans les slides (7 et 8) mais je n’ai pas interrompu la vidéo et j’essayais de continuer comme si de rien n’était. Pour la proportion \widehat{p}, on a juste la loi asymptotique Gaussienne (ou une transformation de la loi Binomiale à distance finie, mais pas de loi de Student, désolé pour le copier-coller depuis des vieux slides qui parlaient d’estimation de moyenne \overline{x} pour des observations supposées Gaussiennes)

J’essaye de limiter le nombre d’erreurs, mais ces cours à distance me prennent beaucoup de temps, et je n’ai pas le courage pour refaire les vidéos, pour l’instant. Toutes mes excuses.

 

 

Reprise des cours (à venir)

Lundi, les cours de la session d’hiver reprennent… à distance. Je vais mettre en ligne une série de capsules vidéos pour finir le cours, sur les GLM. J’ai mis en ligne une première vidéo (slides 0) pour annoncer le plan. J’ai fait des slides rapidement (ça changera des cours que je faisais au tableau) et je fais des enregistrement unique, sans montage, histoire de mettre en ligne rapidement le cours en ligne. Le pdf des slides est aussi en ligne (slides 0). J’essayerais de mettre les vidéos en ligne au fur et à mesure. Je ne suis pas particulièrement fier de moi, mais quitte à perdre du temps à faire le guignol devant la caméra, autant que ça serve au plus grand nombre.

Ah oui, il y a probablement des coquilles, voire des erreurs dans les slides… les commentaires sont ouverts pour me faire part de toute suggestion !

Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic

Vendredi dernier, j’ai mis en ligne un petit quiz pour prédire qui a survécu au naufrage du Titanic, en 1912. J’avais mis en ligne une sous-base avec quelques passagers

loc_fichier = "http://freakonometrics.free.fr/titanic.RData"
download.file(loc_fichier, "titanic.RData")
load("titanic.RData")
str(base)

et plus précisément

  • Survived: Passenger survival indicator (1 if survived)
  • Pclass: Passenger class
  • Sex: Sex of the passenger
  • Age: Age of the passenger
  • SibSp: Number of siblings/spouses aboard
  • Parch: Number of parents/children aboard
  • Embarked: Port of Embarkation
  • Name: Name of passenger

La dernière variable ne sert pas à grand chose, donc on la vire,

base = base[,1:7]

On va maintenant pouvoir répondre aux questions.

Continue reading Il était prévisible que Léonard DiCaprio ne survive pas au naufrage du Titanic

Qui a survécu au naufrage du Titanic?

also known as quiz numéro 4 du cours STT5100 (de la session d’hiver). Vendredi dernier, alors que nous terminions le cours vers midi, François Legault a décrété l’état d’urgence sanitaire pour une quinzaine de jours. Le cours est donc sur la glace depuis une semaine, et (pour l’instant) encore une bonne semaine.

Les deux derniers cours ont été l’occasion de voir comment faire, et comment analyser, une régression logistique. Histoire de briser un peu le confinement que l’on connaît depuis une semaine, j’ai envoyé ce matin un quiz aux étudiants, et j’en profite pour le mettre en ligne (la correction arrivera dans le milieu de la semaine prochaine) pour tous ceux qui voudraient se pratiquer un peu sur les modèles de classification…

J’ai mis en ligne une base de données avec des informations sur les passagers du Titanic, qui est un paquebot qui a fait naufrage il y a plus d’un siècle (certains ont peut être vu le documentaire)

loc_fichier = "http://freakonometrics.free.fr/titanic.RData"
download.file(loc_fichier, "titanic.RData")
load("titanic.RData")
str(base)

On a dans la base

  • Survived: Passenger survival indicator (1 if survived)
  • Pclass: Passenger class
  • Sex: Sex of the passenger
  • Age: Age of the passenger
  • SibSp: Number of siblings/spouses aboard
  • Parch: Number of parents/children aboard
  • Embarked: Port of Embarkation
  • Name: Name of passenger

La dernière variable ne sert pas à grand chose, donc on la vire,

base = base[,1:7]

J’ai ensuite posé une série de questions,

  1. quelle proportion de passagers ont survécu ?
  2. quelle proportion de passagers de première classe ont survécu ?
  3. quelle est la statistique de test pour un test d’indépendance du chi-deux [vu en cours vendredi dernier] entre le fait de survivre (ou pas) et la classe ?
  4. ensuite il faut faire une régression logistique, et prévoir la probabilité de survie pour deux passagers (fictifs) du paquebot
newbase = data.frame(
  Pclass = as.factor(c(1,3)),
  Sex = as.factor(c("female","male")),
  Age = c(17,20),
  SibSp = c(1,0),
  Parch = c(2,0),
  Embarked = as.factor(c("S","S")),
  Name = as.factor(c("Winslet, Miss. Kate","DiCaprio, Mr. Leonardo")))

Je ne peux mettre le formulaire de réponse en ligne ici, mais les commentaires sont ouverts pour ceux qui voudraient d’entraîner un peu pendant cette période un peu longue de confinement…

De la qualité d’un classifieur

On va profiter de la quarantaine pour mettre en ligne un billet sur la courbe ROC, la receiver operating characteristic. Considérons une petite base de données avec n=10 observations, deux variables continues, x_1 et x_2, et la variable d’intérêt binaire y\in\{0,1\}. On peut représenter les points dans le plan (x_1,x_2), et on utilise une couleur différente pour y\in\{0,1\}.

x1 = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
x2 = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
y = c(1,1,1,1,1,0,0,1,0,0)
df = data.frame(x1=x1,x2=x2,y=as.factor(y))
plot(x1,x2,col=c("red","blue")[1+y],pch=19,cex=1.5)

On peut alors faire une régression logistique, de telle sorte que \mathbb{P}(Y=1|x_1,x_2)=\frac{e^{\beta_0+\beta_1x_1+\beta_2x_2}}{1+e^{\beta_0+\beta_1x_1+\beta_2x_2}}On peut visualiser l’ensemble des points (x_1,x_2) pour lesquels \mathbb{P}(Y=0|x_1,x_2)=\mathbb{P}(Y=1|x_1,x_2)(on a alors autant de chance d’être rouge et bleu) soit \beta_0+\beta_1x_1+\beta_2x_2=0 qui correspond à une droite,

reg = glm(y~x1+x2,data=df,family=binomial(link = "logit"))
b = coefficients(reg)
abline(a=-b[1]/b[3],b=-b[2]/b[3])

On peut alors représenter y_i en fonction du score, i.e. l’estimation de \mathbb{P}(Y=1|x_{1,i},x_{2,i}),

Y = df$y
S = predict(reg,type="response")
plot(S,y,xlab="probabilité prédite",ylab="y")

On va alors se donner un seuil (par exemple 50\%) : si la probabilité que Y prenne la valeur 1 excède le seuil, on prédit 1 (et sinon 0). Sur la figure ci-dessus, on a alors 4 sortes de points : ceux à gauche du seuil (et qui sont prédits 0), qui sont bien classés s’ils sont en bas, et mal classés en haut; à droite du seuil (et qui sont prédits 1), ils sont bien classés s’ils sont en haut, et mal classés en bas [dans le code ci-dessous, le symbole &gt désigne l’opérateur “supérieur“, qui malheureusement ne passe pas dans cet éditeur]

seuil = .5
Yhat = (S>seuil)*1
plot(S,y,xlab="probabilité prédite",ylab="y",pch=19,
     col=c("red","blue")[1+(y==Yhat)])
abline(v=seuil,lty=2)

Les couleurs reflètent le bon ou mauvais classement : les points rouges correspondent à des erreurs de classement. On peut retrouver tout ça dans le tableau de contingence ci-dessous,  qui correspond au tableau standard d’un test d’hypothèse

table(Yhat,Y)
    Y
Yhat 0 1
   0 3 1
   1 1 5

Ce qui va nous intéresser ici à deux grandeurs particulières : le taux de faux positifs et le taux de vrais positifs,

  FP=sum((Ps==1)*(Y==0))/sum(Y==0)    
  TP=sum((Ps==1)*(Y==1))/sum(Y==1)

On a obtenu ce tableau à un seuil donné (ici 50\%) mais on peut regarder ce qui se passe lors que le seuil change, comme sur l’animation ci-dessous, où on trace, à droite, le taux de vrais positifs (sur l’axe y) en fonction du taux de faux positifs (sur l’axe x

L’ensemble des points donne la courbe ROC.

roc.curve=function(s,print=FALSE){
  Ps=(S>s)*1 
  FP=sum((Ps==1)*(Y==0))/sum(Y==0)    
  TP=sum((Ps==1)*(Y==1))/sum(Y==1)   
  if(print==TRUE){  
    print(table(Observed=Y,Predicted=Ps))   
  }
  vect=c(FP,TP)
  names(vect)=c("FPR","TPR") 
  return(vect)}
u = seq(0,1,length=251)
V = Vectorize(roc.curve)(u)
plot(t(V),type="s",xlab="Faux Positifs",ylab="Vrais Positifs")
segments(0,0,1,1,col="light blue")

On peut vérifier que le point qu’on avait obtenu avec un seuil de 50\% est bien sur la courbe

table(Yhat,Y)
    Y
Yhat 0 1
   0 3 1
   1 1 5
(FP = sum((Yhat)*(Y==0))/sum(Y==0))
[1] 0.25
(TP = sum((Yhat==1)*(Y==1))/sum(Y==1))
[1] 0.8333333
abline(v=FP,lty=2,col="blue")
abline(h=TP,lty=2,col="blue")
points(FP,TP,pch=19,cex=1.5)

Bien entendu, il y a (beaucoup) de packages R qui permettent d’avoir cette courbe,

library(ROCR)
pred = prediction(S,Y)
plot(performance(pred,"tpr","fpr"))

Une grandeur intéressante est appelée aire sous la courbe (ou AUC) qu’on peut calculer ici à la main (on a une simple fonction en escalier)

p1 = roc.curve(1/3)
p2 = roc.curve(.7)
p2[1]*p2[2]+(p1[1]-p2[1])*p1[2]+(1-p1[1]) 
[1] 0.875


mais qu’on peut avoir automatiquement

auc.perf = performance(pred, measure = "auc")
auc.perf@y.values[[1]]
[1] 0.875

Allez, tentons un autre classifieur : toujours une régression logistique, mais sur un facteur obtenu en coupant la seconde variable en deux, \boldsymbol{1}_{[s,\infty)}(x_2)

reg = glm(y~I(x2>.525),data=df,family=binomial(link = "logit"))
abline(h=.525)

La droite horizontale n’est plus la droite qui donne autant de chance d’être rouge que bleu, mais qui coupe la variable x_2. Ici, on prédit juste deux valeurs : 40\% de chance d’être bleu en bas, et 80\% de chance d’être bleu, en haut. Si on représente observations y_i en fonction des probabilités prédites, on obtient

Y = df$y
S = predict(reg,type="response")
plot(S,y,xlab="probabilité prédite",ylab="y",xlim=0:1)

Avec un seuil à 50\%, on obtient le tableau de contingence suivant (avec 3 erreurs, contre 2 auparavant)

seuil = .5
Yhat = (S>seuil)*1
table(Yhat,Y)
    Y
Yhat 0 1
   0 3 2
   1 1 4

Si on trace la courbe ROC, on obtient

roc.curve=function(s,print=FALSE){
  Ps=(S>s)*1 
  FP=sum((Ps==1)*(Y==0))/sum(Y==0)    
  TP=sum((Ps==1)*(Y==1))/sum(Y==1)   
  if(print==TRUE){  
    print(table(Observed=Y,Predicted=Ps))   
  }
  vect=c(FP,TP)
  names(vect)=c("FPR","TPR") 
  return(vect)}
u = seq(0,1,length=251)
V = Vectorize(roc.curve)(u)
plot(t(V),type="l",xlab="Faux Positifs",ylab="Vrais Positifs")
segments(0,0,1,1,col="light blue")

Cette fois, la courbe n’est plus constante par morceaux, mais linéaire par morceaux, et continue… L’interprétation est un peu plus subtile: cette fois, on a deux régions de l’espace, et dans chaque région, on ne sait pas trop comment distinguer (la probabilité est plate partout ici, contraitement à la régression précédante). Autrement dit, dans cette région, on a une proba constante, par exemple 40\% (en bas) : quand on doit prévoir pour un individu dans cette région, on lui attribue les valeurs \{0,1\} respectivement avec les probabilités \{40\%,60\%\}. Quand on a une probabilité constante, on parle de classifieur aléatoire… La diagonale bleue sur la figure ci-dessus est justement un classifieur aléatoire… C’est ce qu’on obtient si on prédit au hasard

pred = prediction(S,Y)
plot(performance(pred,"tpr","fpr"))

Le point est obtenu avec un seuil de 50\% (ou en fait, n’importe quelle valeur entre 40\% et 80\%). On peut là encore calculer l’aire sous la courbe, cette fois à l’aide de trapèzes (ou de triangles)

p1 = roc.curve(.5)
p2[1]*p2[2]/2+(1-p1[1])*p1[2]+(1-p1[1])*(1-p1[2])/2
[1] 0.7083333
auc.perf = performance(pred, measure = "auc")
auc.perf@y.values[[1]]
[1] 0.7083333

Function basis and regression

In the first part of the course on linear models, we’ve seen how to construct a linear model when the vector of covariates \boldsymbol{x} is given, so that \mathbb{E}(Y|\boldsymbol{X}=\boldsymbol{x}) is either simply \boldsymbol{x}^\top\boldsymbol{\beta} (for standard linear models) or a functional of \boldsymbol{x}^\top\boldsymbol{\beta} (in GLMs). But more generally, we can consider transformations of the covariates, so that a linear model can be used. In a very general setting, consider \sum_{j=1}^m\beta_j h_j(\boldsymbol{x})with h_j:\mathbb{R}^p\rightarrow\mathbb{R}. The standard linear model is obtained when m=p and h_j(\boldsymbol{x})=x_j , but of course, much more general models can be obtained, for instance with h_k(\boldsymbol{x})=x_j^2 or h_k(\boldsymbol{x})=x_{j}x_{j'}, that could be used to achieve high-order Taylor expansions. In that case, we will obtain the polynomial regression, that we will discuss first. We might also think of piecewise constant functions, h_k(\boldsymbol{x})=\boldsymbol{1}(x_j\in [a,b]) , that could be related to regression trees (but that is not in the scope in the STT5100 course). And if we go on step futher, we might think of piecewise linear or piecewise polynomial function, possibly with additional continuity constraints, that will lead us to spline basis.

  • Polynomial regression

For pedagogical purpose, when I talk about polynomial regression, I always have in mind (in the univariate case) y=\beta_0+\beta_1x+\beta_2x^2+\cdots+\beta_kx^k+\varepsilonbut if we use

lm(y~poly(x,k))

in R, the output is not the \beta_j‘s.

As discussed in Kennedy & Gentle (1980) Statistical Computing,

Recall that orthogonal polynomials are defined with respect to the classical inner-product (on the finite interval (a,b)){\displaystyle \langle f,g\rangle =\int _{a}^{b}f(x)g(x)~\mathrm {d} x} And a sequence of orthogonal polynomials is (P_n) where P_n is a polynomial of degree n, for all n, and such that P_m\perp P_n for all m\neq n. Note that those polyomials are orthogonal with respect to the inner product defined above, i.e. given some finite interval (a,b). But if (a,b) changes, the polynomials will be different.

A popular family of orthogonal polynomial, on finite interval (-1,+1) is the family of Legendre polynomials, satisfying{\displaystyle \int _{-1}^{1}P_{m}(x)P_{n}(x)~\mathrm {d} x=0}as soon as m\neq n. Those polynomials satisfy Bonnet’s recursion formula{\displaystyle (n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x)} or Rodrigues’ formula {\displaystyle P_{n}(x)={\frac {1}{2^{n}n!}}{\frac {d^{n}}{dx^{n}}}(x^{2}-1)^{n}}The first values are here{\displaystyle P_{0}(x)=1} {\displaystyle P_{1}(x)=x}{\displaystyle P_{2}(x)={\frac {3x^{2}-1}{2}}}{\displaystyle P_{3}(x)={\frac {5x^{3}-3x}{2}}} {\displaystyle P_{4}(x)={\frac {35x^{4}-30x^{2}+3}{8}}}

Interestingly, we can get those polynomial functions using

library(orthopolynom)
(leg4coef = legendre.polynomials(n=4))
[[1]]
1 
 
[[2]]
x 
 
[[3]]
-0.5 + 1.5*x^2 
 
[[4]]
-1.5*x + 2.5*x^3 
 
[[5]]
0.375 - 3.75*x^2 + 4.375*x^4

Of course, there are many families of orthogonal polynomials (Jacobi polynomials, Laguerre polynomials, Hermite polynomials, etc). Now, in R, there is the standard poly function, that we use in polynomial regression.

x = seq(-1,1,length=101)
y = poly(x,4)
y
                   1            2             3            4
  [1,] -1.706475e-01  0.215984813 -2.480753e-01  0.270362873
  [2,] -1.672345e-01  0.203025724 -2.183063e-01  0.216290298
...
[100,]  1.672345e-01  0.203025724  2.183063e-01  0.216290298
[101,]  1.706475e-01  0.215984813  2.480753e-01  0.270362873
attr(,"coefs")
attr(,"coefs")$alpha
[1] 3.157229e-17 2.655145e-16 9.799244e-17 5.368224e-16
 
attr(,"coefs")$norm2
[1]   1.0000000 101.0000000  34.3400000   9.3377328   2.4472330   0.6330176
 
attr(,"degree")
[1] 1 2 3 4
attr(,"class")
[1] "poly"   "matrix"

But these are not Legendre polynomials… As explained in 李哲源‘s post on stackoverflow, the idea is to start with P_{-1}(x)=0, P_{0}(x)=1 and P_{1}(x)=x, and then define \ell_n=\langle P_n,P_n\rangle  as well as \alpha_n=\langle P_nP_1,P_1\rangle/\ell_n=\langle P_n^2,P_1\rangle/\ell_i= and \beta_n=\ell_n/\ell_{n-1}. Finally, define recursively{\displaystyle P_{n}(x)=(x-\alpha_{n-1})P_{n-1}(x)-\beta_{i-1}P_{i-2}(x)}and its normalized version, \tilde{P}_{n}=P_n/\sqrt{\ell_n}. That is what poly computes.

So, for pedagogical purpose, I said that I like to use y=\boldsymbol{x}^\top\boldsymbol{\beta}+\varepsilon where\boldsymbol{x}=(1,x,x^2,\cdots,xˆ{k-1},x^k)And actually, when using poly, we use the QR decomposition of that matrix. As discussed in in 李哲源‘s post, we can almost reproduce the poly function using

my_poly - function (x, degree = 1) {
    xbar = mean(x)
    x = x - xbar
    QR = qr(outer(x, 0:degree, "^"))
    X = qr.qy(QR, diag(diag(QR$qr), length(x), degree + 1))[, -1, drop = FALSE]
    X2 = X * X
    norm2 = colSums(X * X)   
    alpha = drop(crossprod(X2, x)) / norm2
    beta = norm2 / (c(length(x), norm2[-degree]))
    colnames(X) = 1:degree
    scale = sqrt(norm2)
    X = X * rep(1 / scale, each = length(x))
    X}

Nevertheless, the two models are equivalent. More precisely,

plot(cars)
reg1 = lm(dist~speed+I(speed^2)+I(speed^3),data=cars)
reg2 = lm(dist~poly(speed,3),data=cars)
u = seq(3,26,by=.1)
v1 = predict(reg1,newdata=data.frame(speed=u))
v2 = predict(reg2,newdata=data.frame(speed=u))
lines(u,v1,col="blue")
lines(u,v2,col="red",lty=2)

We have exactly the same prediction here

v1[u==15]
     121 
38.43919 
v2[u==15]
     121 
38.43919

And probably also quite interesting : the coefficients do not have the same interpretation (since we do not have the same basis), but the p-value for the highest degree is exactly the same here ! Here the two models reject, with the same confidence, the polynomial of degree three,

summary(reg1)
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -19.50505   28.40530  -0.687    0.496
speed         6.80111    6.80113   1.000    0.323
I(speed^2)   -0.34966    0.49988  -0.699    0.488
I(speed^3)    0.01025    0.01130   0.907    0.369
 
Residual standard error: 15.2 on 46 degrees of freedom
Multiple R-squared:  0.6732,	Adjusted R-squared:  0.6519 
F-statistic: 31.58 on 3 and 46 DF,  p-value: 3.074e-11
 
summary(reg2)
 
Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)        42.98       2.15  19.988  < 2e-16 ***
poly(speed, 3)1   145.55      15.21   9.573  1.6e-12 ***
poly(speed, 3)2    23.00      15.21   1.512    0.137    
poly(speed, 3)3    13.80      15.21   0.907    0.369    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.2 on 46 degrees of freedom
Multiple R-squared:  0.6732,	Adjusted R-squared:  0.6519 
F-statistic: 31.58 on 3 and 46 DF,  p-value: 3.074e-11
  • B-splines regression (and GAMs)

Splines are also important in regression models, especially when we start talking about Generalized Additive Models. See Perperoglou, Sauerbrei, Abrahamowicz & Schmid (2019) for a review. In the univariate case, I introduce (linear) splines through positive parts, in the sense thaty=\beta_0+\beta_1x+\beta_2(x-s_1)_++\cdots+\beta_k(x-s_{k-1})_++\varepsilonwhere (x-s)_+ equals 0 if x<s and x-s if x>s. Those functions are nice since they are continuous, so the model is continuous (the weighted sum of continuous functions is continuous). And we can go one step further, with y=\beta_0+\beta_1x+\beta_2x^2+\beta_3(x-s_1)^2_++\cdots+\beta_k(x-s_{k-2})^2_++\varepsilonwith quadratic splines, or y=\beta_0+\beta_1x+\beta_2x^2+\beta_3x^3+\beta_4(x-s_1)^3_++\cdots+\beta_k(x-s_{k-3})^3_++\varepsilonfor cubic splines. Interestingly, quadratic splines are not only continuous, but their first derivative is also continuous (and the second one for cubic splines). So the knot discontinuity is s_1,s_2,\cdots is now invisible…

I like those models since they are easy to interprete. For example, the simple model \beta_1 x+\beta_2(x-s)_+ is the following piecewise linear function, continuous, with a “rupture” at knot s.

Observe also the following interpretation: for small values of x, there is a linear increase, with slope \beta_1, and for lager values of x, there is a linear decrease, with slope \beta_1+\beta_2. Hence, \beta_2 is interpreted as a change of the slope.

Unfortunately, it is now what R is using when using the bs function in R, which are the standard B-splines. Just to visualize (I will skip the maths here), with R, we have

library(splines)
clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02")
x = seq(5,25,by=.25)
B = bs(x,knots=c(10,20),Boundary.knots=c(5,55),degre=1)
matplot(x,B,type="l",lty=1,lwd=2,col=clr6)
B=bs(x,knots=c(10,20),Boundary.knots=c(5,55),degre=2)
matplot(x,B,type="l",col=clr6,lty=1,lwd=2)

while the functions I mentioned were (more or less) the following

pos = function(x,s) (x-s)*(x&gt;s)
par(mfrow=c(1,2))
clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02")
x = seq(5,25,by=.25)
B = cbind(pos(x,5),pos(x,10),pos(x,20))
matplot(x,B,type="l",lty=1,lwd=2,col=clr6)
pos2 = function(x,s) (x-s)^2*(x&gt;s)
B = cbind(pos(x,5)*20,pos2(x,5),pos2(x,10),pos2(x,20))
matplot(x,B,type="l",col=clr6,lty=1,lwd=2)

And as for the polynomial regression, the two models are equivalent. For example

plot(cars)
reg1 = lm(dist~speed+pos(speed,10)+pos(speed,20),data=cars)
reg2 = lm(dist~bs(speed,degree=1,knots=c(10,20)),data=cars)
v1 = predict(reg1,newdata=data.frame(speed=u))
v2 = predict(reg2,newdata=data.frame(speed=u))
lines(u,v1,col="blue")
lines(u,v2,col="red",lty=2)

or more specifically

v1[u==15]
     121 
39.35747 
v2[u==15]
     121 
39.35747

So one more time, the two models are equivalent, but I still find the approach with the positive part more intuitive, and easy to understand. As well as the interpretation of coefficients,

summary(reg1)
 
Coefficients:
               Estimate Std. Error t value Pr(&gt;|t|)  
(Intercept)     -7.6305    16.2941  -0.468   0.6418  
speed            3.0630     1.8238   1.679   0.0998 .
pos(speed, 10)   0.2087     2.2453   0.093   0.9263  
pos(speed, 20)   4.2812     2.2843   1.874   0.0673 .
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11
 
summary(reg2)
 
Coefficients:
                                          Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)                                  4.621      9.344   0.495   0.6233    
bs(speed, degree = 1, knots = c(10, 20))1   18.378     10.943   1.679   0.0998 .  
bs(speed, degree = 1, knots = c(10, 20))2   51.094     10.040   5.089 6.51e-06 ***
bs(speed, degree = 1, knots = c(10, 20))3   88.859     12.047   7.376 2.49e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11

Here we can see directly that the first knot was not interesting (the slope did not change significantly) while the second one was…

De la pratique de la régression

Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne.

Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données,

download.file("http://freakonometrics.free.fr/data3.RData","data3.RData")
load("data3.RData")
str(df)
'data.frame':	147 obs. of  10 variables:
 $ Y : num  11.72 15.91 14.19 11.15 8.31 ...
 $ X1: num  1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...
 $ X2: num  3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ...
 $ X3: num  1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...
 $ X4: num  -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ...
 $ X5: num  0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...
 $ X6: num  -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ...
 $ X7: num  1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...
 $ X8: num  8.31 8.52 9.78 7.34 7.26 ...
 $ X9: num  6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ...

Je vais mettre ici les questions que je posais, et donner des pistes de réflexions, non pas sur les réponses attendues (je n’attends rien de cet exercice à part une réflexion), mais sur la discussion que peut amener chacune des questions,

  • Faîtes un modèle linéaire pour expliquer Y en fonction des neuf variables explicatives. Combien de variables explicatives garderiez-vous ?

Commençons par faire une régression sur toutes les variables

summary(lm(Y~., data=df))
 
Coefficients:
             Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  8.046133   2.193448   3.668 0.000349 ***
X1           0.342293   0.915894   0.374 0.709186    
X2          -0.040479   0.103409  -0.391 0.696073    
X3           1.683875   0.897278   1.877 0.062693 .  
X4          -0.009254   0.062382  -0.148 0.882295    
X5          -1.085367   0.113840  -9.534  &lt; 2e-16 ***
X6           0.983207   0.111830   8.792 5.49e-15 ***
X7          -0.015646   0.087483  -0.179 0.858327    
X8           0.012165   0.094756   0.128 0.898033    
X9           0.172210   0.180605   0.954 0.342005    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.019 on 137 degrees of freedom
Multiple R-squared:  0.9228,	Adjusted R-squared:  0.9178 
F-statistic:   182 on 9 and 137 DF,  p-value: &lt; 2.2e-16

On a, dans la sortie, une dizaine de tests de Student qui sont évoqués, correspondant au test de l’hypothèse H_0:\beta_j=0 (contre l’hypothèse alternative (bilatérale) H_1:\beta_j\neq 0), dans un modèle de la forme y_i=\beta_0+\beta_1x_{1,i}+\dots+\beta_9x_{9,i}+\varepsilon_iC’est ce qu’on appelle le test de significativité de la variable x_j (oui, comme on l’a vu en cours, on peut dire le test parce que les autres tests classiques – Fisher, ou Wald – sont équivalent – sauf qu’au lieu de regarder t on regarde la statistique t^2 – ce qui présente l’avantage de voir la statistique de test comme une forme de distance à l’hypothèse H_0 : si c’est trop grand, on rejette…). Avec un seuil d’acceptation de l’ordre de 5%, on nous dit que 6 variables ne sont pas significatives. Mais gardons bien en mémoire que le test de significativité de x_j est fait ici en supposant que toutes les autres variables restent dans le modèle. Autre chose: avec un seuil d’acceptation de l’ordre de 10%, une des variables (la troisième) peut être vue comme significative.

Faisons un test multiple, pour savoir si on peut supprimer 6 des 9 variables explicatives possibles (faisons le à la main, inutile d’aller chercher un package pour le faire)

reg1 = lm(formula = Y ~ ., data = df)
reg0 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg0)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.11031    0.16788   54.27   &lt;2e-16 ***
X3           1.96784    0.07604   25.88   &lt;2e-16 ***
X5          -0.99296    0.08751  -11.35   &lt;2e-16 ***
X6           1.08761    0.05934   18.33   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: &lt; 2.2e-16 anova(reg0,reg1) Analysis of Variance Table Model 1: Y ~ X3 + X5 + X6 Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 Res.Df RSS Df Sum of Sq F Pr(&gt;F)
1    143 146.16                           
2    137 142.12  6    4.0338 0.6481 0.6916

Le test de Fisher nous dit qu’on peut accepter l’hypothèse que les 6 coefficients sont nuls – ici H_0:\beta_1=\beta_2=\beta_4=\beta_7=\beta_8=\beta_9=0(qui est un test multiple, contrairement au test de Student précédant). Mais il nous dit aussi, qu’individuellement, les trois variables restants semblent significatives. Donc j’aurais tendance à garder 3 variables explicatives (je ne parle pas de la constante : on garde toujours la constante – qui n’explique pas grand chose, sauf la valeur moyenne de y).

  • faites une prévision pour un individu dont on sait que X3=1, X5=1 et X6=8

(en réalité, la vraie question que j’ai posée contenait une typo ce qui la rendait vicieuse parce que ce n’est pas le modèle qu’on vient de construire… mais pour commencer, regardons cette question)

On vient de calibrer ce modèle, donc il suffit de faire une prévision,

predict(reg0, newdata=data.frame(X3=1, X5=1, X6=8))
       1 
18.78603

Mais ce n’était pas la vraie question…

  • faites une prévision pour un individu dont on sait que X1=3, X5=1 et X6=8

Je pense que pour répondre à cette question, il convient d’oublier tout ce qu’on vient de voir. On nous donne trois informations, et on va voir si on peut les exploiter. Autrement dit, on va commencer par regarder la régression sur ces trois variables,

reg156 = lm(formula = Y ~ X1+X5+X6, data = df)
summary(reg156)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.06688    0.17222   52.65   &lt;2e-16 ***
X1           1.98342    0.07800   25.43   &lt;2e-16 ***
X5          -1.01438    0.08947  -11.34   &lt;2e-16 ***
X6           1.07824    0.06039   17.86   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.026 on 143 degrees of freedom
Multiple R-squared:  0.9183,	Adjusted R-squared:  0.9166 
F-statistic: 535.8 on 3 and 143 DF,  p-value: &lt; 2.2e-16

Le modèle est ici bon, les trois variables étant significatives. Cela dit, je dis qu’il est “bon” mais on verra vendredi comment discuter davantage ce point… En tous cas, on peut tenter de faire une prévision, et on obtient

predict(reg156, newdata=data.frame(X1=3, X5=1, X6=8))
       1 
22.62867
  • faites une prévision pour un individu dont on sait que X1=3, X3=2, X5=1 et X6=8

Comme auparavant, on nous donne 4 informations… on va regarder le modèle avec les 4 variables

reg1356 = lm(formula = Y ~ X1+X3+X5+X6, data = df)
summary(reg1356)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.10458    0.17131  53.148   &lt;2e-16 ***
X1           0.16362    0.89028   0.184    0.854    
X3           1.80663    0.88052   2.052    0.042 *  
X5          -0.99558    0.08895 -11.192   &lt;2e-16 ***
X6           1.08649    0.05986  18.152   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.014 on 142 degrees of freedom
Multiple R-squared:  0.9207,	Adjusted R-squared:  0.9184 
F-statistic: 411.9 on 4 and 142 DF,  p-value: &lt; 2.2e-16

Cette fois, une des variables n’est pas significative (la première) donc on devrait l’enlever: on fait alors la régression juste sur les trois autres variables

reg356 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg356)
 
Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)  9.11031    0.16788   54.27   &lt;2e-16 ***
X3           1.96784    0.07604   25.88   &lt;2e-16 ***
X5          -0.99296    0.08751  -11.35   &lt;2e-16 ***
X6           1.08761    0.05934   18.33   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: &lt; 2.2e-16

Cette fois, le modèle est “bon” et on peut alors faire la prévision

predict(reg356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.75388

A titre d’information, si on avait fait la prévision sans enlever la variable non significative, on aurait obtenu la valeur suivante

predict(reg1356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.90501

qui est globalement assez proche. On aura l’occasion d’en reparler en cours, quand on se demandera s’il est plus grave d’avoir une valeur non-significative dans la régression, ou d’oublier une variable importante… Mais encore une fois, le but de ces petits exercices est d’appliquer ce qu’on a vu en cours, et d’introduire des questions auxquelles j’apporterai des réponses au prochain cours !

Combiner les modalités d’une variable factorielle

Un billet rapide pour reprendre un point qu’on a vu ce matin en cours STT5100 pour illustrer le test de Fisher. On va utiliser les données de prix d’appartements en Pologne (données pas mal utilisées dans mon ébauche de notes de cours)

library(DALEX)
data(apartments)
with(data = apartments, boxplot(m2.price ~ district))

On a envie de faire ici des regroupements de modalités (c’est d’ailleurs suggéré par la régression simple, 5 variables explicatives étant ici non significatives). Pour mieux voir, on peut réordonner les modalités

A = with(data = apartments, aggregate(m2.price,by=list(district),FUN=mean))
A = A[order(A$x),]
L = as.character(A$Group.1)
apartments$district = factor(apartments$district, level=L)
with(data = apartments, boxplot(m2.price ~ district))

On va prendre ici le district le moins cher comme référence,

reg=lm(m2.price ~ district, data=apartments)
&gt; summary(reg)
 
Coefficients:
                    Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept)          2968.36      58.02  51.160   &lt;2e-16 ***
districtBielany        17.38      84.16   0.207    0.836    
districtPraga          26.45      85.12   0.311    0.756    
districtUrsynow        42.01      82.65   0.508    0.611    
districtBemowo         80.10      83.71   0.957    0.339    
districtUrsus         102.01      82.25   1.240    0.215    
districtZoliborz      829.59      83.94   9.884   &lt;2e-16 ***
districtMokotow       887.10      81.86  10.837   &lt;2e-16 ***
districtOchota        987.93      84.16  11.738   &lt;2e-16 ***
districtSrodmiescie  2214.39      83.28  26.591   &lt;2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 597.4 on 990 degrees of freedom
Multiple R-squared:  0.5698,	Adjusted R-squared:  0.5659 
F-statistic: 145.7 on 9 and 990 DF,  p-value: &lt; 2.2e-16

On peut tester si les 5 premières modalités sont nulles, ce qui est un test multiple, et on va utiliser le test de Ficher :

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(&gt;F)
1    995 354051715                           
2    990 353269202  5    782513 0.4386 0.8217

La statistique de Fisher est faible, et avec une p-value de 82%. On peut tenter le diable, et rajouter encore une modalité

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0",
                        "districtZoliborz = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F    Pr(&gt;F)    
1    996 405455409                                  
2    990 353269202  6  52186207 24.374 &lt; 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Mais on a peut être été trop gourmand cette fois. On va regrouper les 6 premières modalités (et appeler A le regroupement de districts). Si on regarde les prix moyens, par districts, on obtient

levels(apartments$district) = c(rep("A",6),levels(apartments$district)[7:11])
with(data = apartments, boxplot(m2.price ~ district))

apartments$district = relevel(apartments$district,"Zoliborz")

On recommence, en mettant le district le moins cher comme référence, et on veut tester si les deux suivants ont des coefficients nuls dans la régression linéaire.

reg=lm(m2.price ~ district, data=apartments)
linearHypothesis(reg, c("districtMokotow = 0",
                        "districtOchota = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(&gt;F)
1    997 355292524                           
2    995 354051715  2   1240809 1.7435 0.1754

Avec une p-value de 17%, on peut accepter de regrouper les trois modalités ensemble. On a alors trois groupes de districts, dont les noms sont A, B et C. On obtient les boîtes à moustaches suivantes

levels(apartments$district) = c("B","A",rep("B",2),"C")
apartments$district = relevel(apartments$district,"A")
with(data = apartments, boxplot(m2.price ~ district))

Je laisse les plus courageux vérifier, mais on a trois districts effectivement différents, et si le but est de prévoir le prix des logements, inutiles d’utiliser un découpage avec 10 modalités, un découpage avec 3 suffit !

De l’abus de notation dans les modèles de régression

De manière un peu rituelle, je commence toujours mon cours de régression en revenant sur un point important de la statistique : les abus de notation !  Car tout le monde utilise les mêmes lettres (surtout les grecques) pour désigner des objets de nature différente. Dans la majorité des livres, on pourra nous dire sur la même page que \widehat{\theta}=2.35 et que \text{Var}(\widehat{\theta})=1.07, autrement dit \widehat{\theta} peut désigner en même temps un nombre (dans le premier cas) et une variable aléatoire (dans le second). C’est pour le moins déroutant ! En fait, la raison est assez simple. La statistique commence toujours par un échantillon \{y_1,y_2,\cdots,y_n\}, des données, des chiffres. Si on reste là, on fait du descriptif. L’étape classique est ensuite de supposer que les observations y_i sont des réalisations de variables aléatoires Y_i, qu’on supposera bien souvent indépendantes et identiquement distribuées. Et \widehat{\theta} sera alors une statistique, c’est à dire une fonction de mes observations. Je peux alors définir \widehat{\theta}=t(y_1,\cdots,y_n) comme étant la statistique observée sur mon échantillon, mais je peux aussi considérer \widehat{\theta}=t(Y_1,\cdots,Y_n), qui est alors une variable aléatoire, mais avec la même notation. Si on voulait aider à comprendre, on utiliserait \widehat{\Theta}, mais bon, les choses sont ce qu’elles sont… Et en économétrie, ça devient rapidement un cauchemar quand on commence à parler des résidus… Autre particularité en statistique, c’est que si on distingue l’espérance et la moyenne (empirique), on a un seul mot pour parler de la variance, que ce soit pour une variable aléatoire, ou un vecteur de \mathbb{R}^n. On aura ainsi \mathbb{E}[Y]=\int y dF(y)et\overline{y}=\widehat{\mathbb{E}}[\boldsymbol{y}]=\frac{1}{n}\sum_{i=1} y_ialors que\text{Var}[Y]=\int [y-\mathbb{E}[Y]]^2 dF(y)et\widehat{\text{Var}}[\boldsymbol{y}]=\frac{1}{n}\sum_{i=1} (y_i-\overline{y})^2

Considérons un problème de régression maintenant, avec un modèle de la forme y_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}+\varepsilon_i. Ici, \varepsilon_i est un nombre réel, inconnu. Dans une écriture matricielle, on a \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}, où cette fois \boldsymbol{\varepsilon} est un vecteur de \mathbb{R}^n (et oui, je suis désolé mais ici \boldsymbol{X} désigne la matrice des covariates, et non pas un vecteur aléatoire… je ferais un billet un jour pour parler du fait que parfois on dit que les \boldsymbol{x} sont donnés et des fois – comme on conditionne suivant \boldsymbol{X}, autrement dit, on les voit comme aléatoires). On peut parfois faire une hypothèse quant à la distribution des résidus. Autrement dit, les \varepsilon_i sont vues comme des réalisations de variables aléatoires \varepsilon_i, ainsi que \boldsymbol{\varepsilon}. On notera ainsi \boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}). Ah oui, autre point juste pour perdre les élèves : \text{Var}(\boldsymbol{\varepsilon})=\boldsymbol{\Sigma} alors que \text{Var}(\varepsilon_i)=\sigma^2… Bon, ici comme on suppose les observations indépendentes, et identiquement distribuées, on supposera que \text{Var}(\boldsymbol{\varepsilon})=\boldsymbol{\Sigma}=\sigma^2\mathbb{I}.

Encore une fois, \boldsymbol{\varepsilon} est (par définition) non observable. Par contre, on peut estimer ces résidus : à partir d’un estimateur \widehat{\boldsymbol{\beta}} de \boldsymbol{\beta}, on peut définir \widehat{\boldsymbol{\varepsilon}}=\boldsymbol{y}-\widehat{\boldsymbol{y}}=\boldsymbol{y}-\boldsymbol{x}^\top\widehat{\boldsymbol{\beta}}Histoire de clarifier, je vais plutôt noter  \widehat{\boldsymbol{e}} ces résidus estimés, en utilisant l’estimateur par moindres carrés de \boldsymbol{\beta}. On peut noter que \widehat{\boldsymbol{e}}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{y} où classiquement \boldsymbol{H}=\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top est la matrice de projection sur l’espace engendré par toutes les combinaisons linéaires des variables explicatives. Mais là encore, on peut voir le vecteur (numérique) \widehat{\boldsymbol{e}} comme la réalisation d’une variable aléatoire \widehat{\boldsymbol{E}}. En particulier, \widehat{\boldsymbol{E}}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{Y}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{\varepsilon}\boldsymbol{\varepsilon} est notre vecteur aléatoire, centré, de matrice de variance-covariance \text{Var}(\boldsymbol{\varepsilon})=\sigma^2\mathbb{I}. On peut alors en déduire que\mathbb{E}[\widehat{\boldsymbol{E}}]=(\mathbb{I}-\boldsymbol{H})\mathbb{E}[\boldsymbol{\varepsilon}]=\boldsymbol{0}et\text{Var}[\widehat{\boldsymbol{E}}]=(\mathbb{I}-\boldsymbol{H})\text{Var}[\boldsymbol{\varepsilon}](\mathbb{I}-\boldsymbol{H})^\top=\sigma^2(\mathbb{I}-\boldsymbol{H})(car (\mathbb{I}-\boldsymbol{H}) est idenpotent). Cette dernière relation est particulièrement importante, car on notera que \text{Var}(\widehat{\boldsymbol{E}})\neq\sigma^2\mathbb{I}. En particulier, si on prend un résidu estimé au hasard \text{Var}(\widehat{E}_i)=\sigma^2(1-\boldsymbol{H}_{i,i}) (on avait parlé longuement de \boldsymbol{H}_{i,i} dans un billet récent, on leverage, en particulier on avait vu que \boldsymbol{H}_{i,i}\in[0,1] (on avait discuté la borne inférieur, qui peut être améliorée, en fait \boldsymbol{H}_{i,i}\in(0,1]) de telle sorte que \text{Var}(\widehat{E}_i)\leq\sigma^2. Si on poursuit un peu, on peut regarder la somme des carrés estimés, et noter que\mathbb{E}\big[\sum_{i=1}^n \widehat{E}_i^2\big]=\mathbb{E}[\text{trace}( \widehat{\boldsymbol{E}}\widehat{\boldsymbol{E}}^\top)] =\text{trace}(\mathbb{E}[\text{trace}( \widehat{\boldsymbol{E}}\widehat{\boldsymbol{E}}^\top])i.e.\mathbb{E}\big[\sum_{i=1}^n \widehat{E}_i^2\big]=\sigma^2\text{trace}(\mathbb{I}-\boldsymbol{H})or \text{trace}(\mathbb{I}-\boldsymbol{H})=n-p donc\widehat{\sigma}^2=\frac{1}{n-p}\sum_{i=1}^n \widehat{E}_i^2est un estimateur sans biais de \sigma^2. Et classiquement, on considèrera les résidus Studentisés\widehat{R}_i=\frac{\widehat{E}_i}{\widehat{\sigma}\sqrt{1-\boldsymbol{H}_{i,i}}}Si je voulais résumer un peu, on pourrait dire que\text{Var}(\boldsymbol{E})=\sigma^2\mathbb{I}\widehat{\text{Var}}(\boldsymbol{E})=\widehat{\sigma}^2\mathbb{I}\text{Var}(\widehat{\boldsymbol{E}})=\sigma^2(\mathbb{I}-\boldsymbol{H})\widehat{\text{Var}}(\widehat{\boldsymbol{E}})=\widehat{\sigma}^2(\mathbb{I}-\boldsymbol{H})En espérant que ça clarifie un peu…(?)

On leverage

Last week, in our STT5100 (applied linear models) class, I’ve introduce the hat matrix, and the notion of leverage. In a classical regression model, \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta} (in a matrix form), the ordinary least square estimator of parameter \boldsymbol{\beta} is \widehat{\boldsymbol{\beta}}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y}The prediction can then be written\widehat{\boldsymbol{y}}=\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\underbrace{\color{blue}{\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}}_{\color{blue}{\boldsymbol{H}}}\boldsymbol{y}where \color{blue}{\boldsymbol{H}} is called the hat matrix.

The matrix is idempotent, i.e. \boldsymbol{H}^2={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\textcolor{grey}{\boldsymbol{X}^\top{\boldsymbol{X}}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}}\boldsymbol{X}^\top}={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}=\boldsymbol{H}so it can be interpreted as a projection matrix. Furthermore, since\boldsymbol{H}\boldsymbol{X}=\boldsymbol{X} (just do the maths), the projection is on a subspace that contains all the linear combinations of columns of \boldsymbol{X}. One can also observe that \mathbb{I}-\boldsymbol{H} is also a projection matrix. And we can write\boldsymbol{y}=\underbrace{\boldsymbol{H}\boldsymbol{y}}_{\widehat{\boldsymbol{y}}}+\underbrace{(\mathbb{I}-\boldsymbol{H})\boldsymbol{y}}_{\widehat{\boldsymbol{\varepsilon}}}where \widehat{\boldsymbol{y}} is the orthogonal projection of \boldsymbol{y} on the (linear) space of linear combinations of columns of \boldsymbol{X}, and \widehat{\boldsymbol{y}}\perp\widehat{\boldsymbol{\varepsilon}}, which gives the classical interpretation of residuals, being unpredictible (at least with a linear model using variables \boldsymbol{X}).

Let’s move a bit faster now (we’ve seen many other properties last week), and consider elements on the diagonal of matrix \boldsymbol{H}. Recall that we have

so entry \boldsymbol{H}_{i,i} is a measure of the influence of entry \boldsymbol{y}_i on its prediction latex]\widehat{\boldsymbol{y}}_i[/latex].

We have seen that\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}(\boldsymbol{H})=\text{trace}(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top)which can be written\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}\boldsymbol{X}^\top(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1})=\text{trace}(\mathbb{I})=pwhere classically p=k+1, where k is the number of explanatory variables. Further, since \boldsymbol{H} is idempotent, we can write (from \boldsymbol{H}=\boldsymbol{H}^2) that\boldsymbol{H}_{i,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}\boldsymbol{H}_{j,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}^2One the one hand, since the second term is positive \boldsymbol{H}_{i,i}\geq\boldsymbol{H}_{i,i}^2, i.e. 1\geq\boldsymbol{H}_{i,i}. And since both terms are positive, then \boldsymbol{H}_{i,i}\in[0,1]. And there was a question in the course on the sharpeness of the bounds.

Using Anscombe’s dataset, we’ve seen that it was possible to get a leverage of 1. Using something rather similar

df = data.frame(x = c(rep(1,10),6), y = c(1:10,8))
plot(df)

we obtain

model = lm(y~x,data=df)
abline(model,col="red",lwd=2)
H = lm.influence(model)$hat
plot(1:11,H,type="h")

The very last observation, the one one the right, is here extremely influencial : if we remove it, the model is completely different ! And here, we reach the upper bound, \boldsymbol{H}_{11,11}=1. Observe that all other points are equally influencial, and because on the constraint on the trace of the matrix, \boldsymbol{H}_{i,i}=1/10 when i\in\{1,2,\cdots,10\}.

Now, what about the lower bound ? In order to have some sort of “non-influencial” observations, consider the two following case.

  • the case where one observation (below the first one) is such that \widehat{\boldsymbol{y}}_{i}=\boldsymbol{y}_{i} (perfect prediction)
  • the case where one observation (below the tenth one) is such that \boldsymbol{x}_{i}=\overline{\boldsymbol{x}} and \boldsymbol{y}_{i}=\overline{\boldsymbol{y}} (from the first order condition – or normal equation), the fitted regression line always go through point (\overline{\boldsymbol{x}},\overline{\boldsymbol{y}})

Let us move two observations from our dataset,

mean(c(4,rep(1,8),6))
[1] 1.8
df = data.frame(x = c(4,rep(1,8),6,1.8),
y = c(predict(model,newdata=data.frame(x=4)),
2:9,8,
predict(model,newdata=data.frame(x=1.8))))

We now have

If we compute the leverages, we obtain

model = lm(y~x,data=df)
H = lm.influence(model)$hat
plot(1:11,H,type="h")

so, for the first observation, its leverage actually increased (the blue part), and for the tenth one, we have the lowest influence, but it is not zero. Is it possible to reach zero ?

Here, observe that for the tenth observation, \boldsymbol{H}_{i,i}=1/n. And actually, that’s the best we can do… We can prove that, in the case of a simple regression (as above)\boldsymbol{H}_{i,i}=\frac{1}{n}+\frac{(x_i-\overline{x})^2}{n\text{Var}(x)}which is minimum when x_i=\overline{x}, and then \boldsymbol{H}_{i,i}=1/n, otherwise \boldsymbol{H}_{i,i}>1/n. And this property is also valid in a multiple regression (as soon as an intercept is included in the regression – which should always be the case). To prove that result, let \tilde{\boldsymbol{X}} denote the matrix of centered variables \boldsymbol{X}, then we can prove that \boldsymbol{H}_{i,i}=\frac{1}{n}+\big[\tilde{\boldsymbol{X}}(\tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}})^{-1}\tilde{\boldsymbol{X}}^\top\big]_{i,i}(which is basically a matrix version of the previous equation).

I can maybe add another comment on Anscombe’s data. We’ve seen that on the right that we did reach 1. But I did not prove it. One way to prove it is actually to focus on the remaining n-1 points, on the left. Those have all the same x values. We can prove that if \boldsymbol{X}_{i_1}=\boldsymbol{X}_{i_2}, then \boldsymbol{H}_{i_1,i_2}=\boldsymbol{X}_{i_1}^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_{i_2}=\boldsymbol{H}_{i_1,i_1}hence, using the relationship obtained since the hat matrix is idempotent\boldsymbol{H}_{i_1,i_1}=2\boldsymbol{H}_{i_1,i_1}^2+\sum_{j\notin\{i_1,i_2\}}\boldsymbol{H}_{i_1,j}^2thus, we now have\boldsymbol{H}_{i_1,i_1}\big(1-2\boldsymbol{H}_{i_1,i_1}\big)>0i.e. \boldsymbol{H}_{i_1,i_1}\in[0,1/2], where the upper bound becomes 1/(n-1) “duplicates”. So for n-1 \boldsymbol{H}_{i,i}‘s, we have values below 1/(n-1), the last one should be below 1 and the sum has to be k=2 . So we have the value of the n \boldsymbol{H}_{i,i}‘s.

 

What it the interpretation of the diagonal for a ROC curve

Last Friday, we discussed the use of ROC curves to describe the goodness of a classifier. I did say that I will post a brief paragraph on the interpretation of the diagonal. If you look around some say that it describes the “strategy of randomly guessing a class“, that it is obtained with “a diagnostic test that is no better than chance level“, even obtained by “making a prediction by tossing of an unbiased coin“.

Let us get back to ROC curves to illustrate those points. Consider a very simple dataset with 10 observations (that is not linearly separable)

x1 = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
x2 = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
y = c(1,1,1,1,1,0,0,1,0,0)
df = data.frame(x1=x1,x2=x2,y=as.factor(y))

here we can check that, indeed, it is not separable

plot(x1,x2,col=c("red","blue")[1+y],pch=19)

Consider a logistic regression (the course is on linear models)

reg = glm(y~x1+x2,data=df,family=binomial(link = "logit"))

but any model here can be used… We can use our own function

Y=df$y
S=predict(reg)
roc.curve=function(s,print=FALSE){
  Ps=(S&gt;=s)*1
   
  FP=sum((Ps==1)*(Y==0))/sum(Y==0)
     
  TP=sum((Ps==1)*(Y==1))/sum(Y==1)if(print==TRUE){print(table(Observed=Y,Predicted=Ps))}
   
vect=c(FP,TP)names(vect)=c("FPR","TPR")return(vect)}

or any R package actually

library(ROCR)

perf=performance(prediction(S,Y),"tpr","fpr")

We can plot the two simultaneously here

plot(performance(prediction(S,Y),"tpr","fpr"))
V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

So our code works just fine, here. Let us consider various strategies that should lead us to the diagonal.

The first one is : everyone has the same probability (say 50%)

S=rep(.5,10)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])

Indeed, we have the diagonal. But to be honest, we have only two points here : (0,0) and (1,1). Claiming that we have a straight line is not very satisfying… Actually, note that we have this situation whatever the probability we choose

S=rep(.2,10)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])

We can try another strategy, like “making a prediction by tossing of an unbiased coin“. This is what we obtain

set.seed(1)

S=sample(0:1,size=10,replace=TRUE)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

We can also try some sort of “random classifier”, where we choose the score randomly, say uniform on the unit interval

set.seed(1)

S=runif(10)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(0,1,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

Let us try to go further on that one. For convenience, let us consider another function to plot the ROC curve

V=Vectorize(roc.curve)(seq(0,1,length=251))

roc_curve=Vectorize(function(x) max(V[2,which(V[1,]&lt;=x)]))

We have the same line as previously

x=seq(0,1,by=.025)

y=roc_curve(x)lines(x,y,type="s",col="red")

But now, consider many scoring strategies, all randomly chosen

MY=matrix(NA,500,length(y))for(i in 1:500){
  
S=runif(10)
  
V=Vectorize(roc.curve)(seq(0,1,length=251))
  
MY[i,]=roc_curve(x)}plot(performance(prediction(S,df$y),"tpr","fpr"),col="white")for(i in 1:500){lines(x,MY[i,],col=rgb(0,0,1,.3),type="s")}lines(c(0,x),c(0,apply(MY,2,mean)),col="red",type="s",lwd=3)segments(0,0,1,1,col="light blue")

The red line is the average of all random classifiers. It is not a straight line, be we observe oscillations around the diagonal.

Consider a dataset with more observations


myocarde = read.table("http://freakonometrics.free.fr/myocarde.csv",head=TRUE, sep=";")

myocarde$PRONO = (myocarde$PRONO=="SURVIE")*1

reg = glm(PRONO~.,data=myocarde,family=binomial(link = "logit"))

Y=myocarde$PRONO

S=predict(reg)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

Here is a “random classifier” where we draw scores randomly on the unit interval

S=runif(nrow(myocarde)plot(performance(prediction(S,Y),"tpr","fpr"))

V=Vectorize(roc.curve)(seq(-5,5,length=251))points(V[1,],V[2,])segments(0,0,1,1,col="light blue")

And if we do that 500 times, we obtain, on average

MY=matrix(NA,500,length(y))for(i in 1:500){
  
S=runif(length(Y))
  
V=Vectorize(roc.curve)(seq(0,1,length=251))
  
MY[i,]=roc_curve(x)}plot(performance(prediction(S,Y),"tpr","fpr"),col="white")for(i in 1:500){lines(x,MY[i,],col=rgb(0,0,1,.3),type="s")}lines(c(0,x),c(0,apply(MY,2,mean)),col="red",type="s",lwd=3)segments(0,0,1,1,col="light blue")

So, it looks like me might say that the diagonal is what we have, on average, when drawing randomly scores on the unit interval…

I did mention that an interesting visual tool could be related to the use of the Kolmogorov Smirnov statistic on classifiers. We can plot the two empirical cumulative distribution functions of the scores, given the response Y

score=data.frame(yobs=Y,
                 ypred=predict(reg,type="response"))

f0=c(0,sort(score$ypred[score$yobs==0]),1)

f1=c(0,sort(score$ypred[score$yobs==1]),1)plot(f0,(0:(length(f0)-1))/(length(f0)-1),col="red",type="s",lwd=2,xlim=0:1)lines(f1,(0:(length(f1)-1))/(length(f1)-1),col="blue",type="s",lwd=2)

we can also look at the distribution of the score, with the histogram (or density estimates)

S=score$ypred

hist(S[Y==0],col=rgb(1,0,0,.2),
     probability=TRUE,breaks=(0:10)/10,border="white")hist(S[Y==1],col=rgb(0,0,1,.2),
     probability=TRUE,breaks=(0:10)/10,border="white",add=TRUE)lines(density(S[Y==0]),col="red",lwd=2,xlim=c(0,1))lines(density(S[Y==1]),col="blue",lwd=2)

The underlying idea is the following : we do have a “perfect classifier” (top left corner)

is the supports of the scores do not overlap

otherwise, we should have errors. That the case below

we in 10% of the cases, we might have misclassification

or even more missclassification, with overlapping supports

Now, we have the diagonal

when the two conditional distributions of the scores are identical

Of course, that only valid when n is very large, otherwise, it is only what we observe on average….

Exotic link functions for GLMs

In my previous post on GLMs, I discussed power link functions. But there are much more links that can be used :

  • The square root link (for the Poisson model)

Consider some random variable Y with mean \mu and variance \sigma^2. Using Taylor’s expansion,g(Y)\sim g(\mu)+(Y-\mu)g'(\mu)+\frac{1}{2}(Y-\mu)^2g''(\mu)we can write\mathbb{E}[g(Y)]\sim g(\mu)+\frac{\sigma^2}{2}g''(\mu) \text{Var}[g(Y)]\sim [g'(\mu)]^2\sigma^2

Assume that Y\sim\mathcal{P}(\lambda), a consider a square root transformation, g(y)=\sqrt{y}, then the second equality becomes \text{Var}[\sqrt{Y}]\sim \left[\frac{1}{2\sqrt{\mathbb{E}[Y]}}\right]^2\text{Var}[Y]=\frac{1}{4}

So, somehow, with a square-root transformation, we have variance stability, which might be interpreted as some homoscedasticity.

  • The complementary log-log function for the Bernoulli model

Assume that the true variable of interest is a Poisson one, N|\mathbf{X}=\mathbf{x}\sim\mathcal{P}(\lambda_{\mathbf{x}}) where \lambda_{\mathbf{x}}=\exp[\mathbf{x}^T\mathbf{\beta}]Thus,\mathbb{P}[N=0|\mathbf{X}=\mathbf{x}]=\exp[-\lambda_{\mathbf{x}}]=\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]while\mathbb{P}[N>0|\mathbf{X}=\mathbf{x}]=1-\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]=H(\mathbf{x}^T\mathbf{\beta})where H(s)=1-\exp[-\exp(s)]. Let Y=\mathbf{1}(N>0). The previous model seems like a Bernoulli regression with H as link function,\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=H(\mathbf{x}^T\mathbf{\beta})

So, assume now that instead of observing N we observe Y=\boldsymbol{1}(N>0). In that case, running a Bernoulli regression with a complementary log-log link function would be the same (?) as running first a Poisson regression on the original data, and then use it on our binary variable, zero vs. non-zero. Let us generate some data, and see what’s going on. Let us compare e^{\lambda_{\mathbf{x}}} and p_{\mathbf{x}} obtained from a standard logistic regression

n=563
set.seed(1)
base=data.frame(X1=rnorm(n),X2=rnorm(n))
lambda=base$X1+base$X2
base$Y=rpois(n,exp(lambda))
regPois = glm(Y~.,data=base,family=poisson(link="log"))
lambda = predict(regPois,type="response")
regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit"))
prob = predict(regBinom, type="response")
plot(prob,exp(-lambda),xlim=0:1,ylim=0:1)
abline(a=0,b=1,lty=2,col="red")

What if p_{\mathbf{x}} was obtained from a Bernoulli regression, with a cloglog link function ?

regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog"))
prob = predict(regBinom, type="response")
plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1)
abline(a=0,b=1,lty=2,col="red")

It looks like the fit is very good here ! Now, what if we have real data, like the dataset from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy (with 563 observations, and nine variables)

base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE)
str(base)
x=base$SEX
base$SEX="M"
base$SEX[x=="0"]="F"
x=base$CHILDREN
base$CHILDREN="YES"
base$CHILDREN[x==0]="NO"
regPois = glm(Y~.,data=base,family=poisson(link="log"))
lambda = predict(regPois,type="response")
regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit"))
prob = predict(regBinom, type="response")
plot(prob,exp(-lambda),xlim=0:1,ylim=0:1)
abline(a=0,b=1,lty=2,col="red")

In that case the two models are very different. But actually, so is the second one

regBinom = glm((Y&gt;0)~.,data=base,family=binomial(link="cloglog"))
prob = predict(regBinom, type="response")
plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1)
abline(a=0,b=1,lty=2,col="red")

How can we interpret that ? Could it be because the Poisson model is not good ? Actually, if we run a zero-inflated model here,

library(pscl)
regZIP = zeroinfl(Y ~ . | ., data = base)
summary(regZIP)
 
Count model coefficients (poisson with log link):
             Estimate Std. Error z value Pr(&gt;|z|)    
(Intercept) -0.002274   0.048413  -0.047    0.963    
X1           1.019814   0.026186  38.945   &lt;2e-16 ***
X2           1.004814   0.024172  41.570   &lt;2e-16 *** 
Zero-inflation model coefficients (binomial with logit link): 
            Estimate Std. Error z value Pr(&gt;|z|)  
(Intercept) -4.90190    2.07846  -2.358   0.0184 *
X1          -2.00227    0.86897  -2.304   0.0212 *
X2          -0.01545    0.96121  -0.016   0.9872  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hence, we reject here the Poisson distribution assumption, because of the inflation of zeros… It looks like the cloglog link can be used to check if the Poisson distribution is a good model, or not…

GLMs: link vs. distribution

Usually, when I give a course on GLMs, I try to insist on the fact that the link function is probably more important than the distribution. In order to illustrate, consider the following dataset, with 5 observations

x = c(1,2,3,4,5)
y = c(1,2,4,2,6)
base = data.frame(x,y)

Then consider several model, with various distributions, and either an identity link (and in that case \mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbf{x}^T\mathbf{\beta}) or a log link function (so that \mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=e^{\mathbf{x}^T\mathbf{\beta}})

regNId = glm(y~x,family=gaussian(link="identity"),data=base)
regNlog = glm(y~x,family=gaussian(link="log"),data=base)
regPId = glm(y~x,family=poisson(link="identity"),data=base)
regPlog = glm(y~x,family=poisson(link="log"),data=base)
regGId = glm(y~x,family=Gamma(link="identity"),data=base)
regGlog = glm(y~x,family=Gamma(link="log"),data=base)
regIGId = glm(y~x,family=inverse.gaussian(link="identity"),data=base)
regIGlog = glm(y~x,family=inverse.gaussian(link="log"),data=base

One can also consider some Tweedie distribution, to be even more general

library(statmod)
regTwId = glm(y~x,family=tweedie(var.power=1.5,link.power=1),data=base)
regTwlog = glm(y~x,family=tweedie(var.power=1.5,link.power=0),data=base)

Consider the prediction obtained in the first case, with the linear link function

library(RColorBrewer)
darkcols = brewer.pal(8, "Dark2")
plot(x,y,pch=19)
abline(regNId,col=darkcols[1])
abline(regPId,col=darkcols[2])
abline(regGId,col=darkcols[3])
abline(regIGId,col=darkcols[4])
abline(regTwId,lty=2)

The predictions are very very close, aren’t they ? In the case of the exponential prediction, we obtain

plot(x,y,pch=19)
u=seq(.8,5.2,by=.01)
lines(u,predict(regNlog,newdata=data.frame(x=u),type="response"),col=darkcols[1])
lines(u,predict(regPlog,newdata=data.frame(x=u),type="response"),col=darkcols[2])
lines(u,predict(regGlog,newdata=data.frame(x=u),type="response"),col=darkcols[3])
lines(u,predict(regIGlog,newdata=data.frame(x=u),type="response"),col=darkcols[4])
lines(u,predict(regTwlog,newdata=data.frame(x=u),type="response"),lty=2)

We can actually look closer. For instance, in the linear case, consider the slope obtained with a Tweedie model (that will include all the parametric familes mentioned here, actually)

pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base))$coefficients[2,1:2]
Vgamma = seq(-.5,3.5,by=.05)
Vpente = Vectorize(pente)(Vgamma)
plot(Vgamma,Vpente[1,],type="l",lwd=3,ylim=c(.965,1.03),xlab="power",ylab="slope")

The slope here is always very very close to one ! Even more if we add a confidence interval

plot(Vgamma,Vpente[1,])
lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2)
lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2)

Heuristically, for the Gamma regression, or the Inverse Gaussian one, because the variance is a power of the prediction, if the prediction is small (here on the left), the variance should be small. So, on the left of the graph, the error should be small with a higher power for the variance function. And that’s indeed what we observe here

erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=1),data=base),newdata=data.frame(x=1),type="response")-y[x==1] 
Verreur = Vectorize(erreur)(Vgamma)
plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(-.1,.04),xlab="power",ylab="error")
abline(h=0,lty=2)

Of course, we can do the same with the exponential models

pente=function(gamma) summary(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base))$coefficients[2,1:2]
Vpente = Vectorize(pente)(Vgamma)
plot(Vgamma,Vpente[1,],type="l",lwd=3)

or, if we add the confidence bands, we obtain

plot(Vgamma,Vpente[1,],ylim=c(0,.8),type="l",lwd=3,xlab="power",ylab="slope")
lines(Vgamma,Vpente[1,]+1.96*Vpente[2,],lty=2)
lines(Vgamma,Vpente[1,]-1.96*Vpente[2,],lty=2)

So here also, the “slope” is rather similar… And if we look at the error we make on the left part of the graph, we obtain

erreur=function(gamma) predict(glm(y~x,family=tweedie(var.power=gamma,link.power=0),data=base),newdata=data.frame(x=1),type="response")-y[x==1] 
Verreur = Vectorize(erreur)(Vgamma)
plot(Vgamma,Verreur,type="l",lwd=3,ylim=c(.001,.32),xlab="power",ylab="error")

So my point is that the distribution is usually not the most important point on GLMs, even if chapters of books on GLMs are distribution based… But as mentioned in an another post, if you consider a nonlinear transformation, like we have with GAMs, the story is more complicated…

Bailey (1963) and Poisson regression on two factors

Consider the following dataset, from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy, with 563 observations, and nine variables : eight covariates, and the variable of interest, the number of extramarital affairs, over a year,

base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE)
str(base)
'data.frame':	563 obs. of  9 variables:
 $ SEX         : int  1 0 0 1 1 0 0 1 0 1 ...
 $ AGE         : num  37 27 32 57 22 32 22 57 32 22 ...
 $ YEARMARRIAGE: num  10 4 15 15 0.75 1.5 0.75 15 15 1.5 ...
 $ CHILDREN    : int  0 0 1 1 0 0 0 1 1 0 ...
 $ RELIGIOUS   : int  3 4 1 5 2 2 2 2 4 4 ...
 $ EDUCATION   : int  18 14 12 18 17 17 12 14 16 14 ...
 $ OCCUPATION  : int  7 6 1 6 6 5 1 4 1 4 ...
 $ SATISFACTION: int  4 4 4 5 3 5 3 4 2 5 ...
 $ Y           : int  0 0 0 0 0 0 0 0 0 0 ...

Let us focus on two categorical covariates, related to the importance of religion, and the occupation

df=data.frame(y=base$Y,
              religion=as.factor(base$RELIGIOUS),
              occupation=as.factor(base$OCCUPATION),
              expo = 1)
(E=xtabs(expo~religion+occupation,data=df))
        occupation
religion  1  2  3  4  5  6  7
       1  4  1  8  4 16  9  0
       2 23  3 11 17 56 36  6
       3 29  1 10 12 39 25  2
       4 38  7 12 21 59 44  2
       5 13  1  3 10 19 19  3
(N=xtabs(y~religion+occupation,data=df))
        occupation
religion  1  2  3  4  5  6  7
       1  4  1 13  3 13  7  0
       2  1  1 13 10 25 43 10
       3 15  0 12 11 34 35  1
       4 24  1  3 15 11  9 10
       5  6  0  0  6 11  7  0

The two tables above are the exposure (number of observations) and the number of extramarital affairs, here as contingency tables. Without any covariate, one can assume that N\sim\mathcal{P}(\lambda\cdot E), where \lambda would be

sum(N)/sum(E)
[1] 0.6305506

The idea with the margin method is to assume that N_{i,j}=E_{i,j}\cdot\lambda_{i,j} where \lambda_{i,j}=A_i\cdot B_j. Bailey (1963) added two series of constraints : per row, \sum_j N_{i,j}=\sum_j E_{i,j}\cdot A_i\cdot B_j for any i and similarly, for any j \sum_i N_{i,j}=\sum_i E_{i,j}\cdot A_i\cdot B_jFrom the first series of constraints, write A_i=\frac{\sum_j N_{i,j}}{\sum_j E_{i,j}\cdot B_j} and use the second series to write B_j=\frac{\sum_i N_{i,j}}{\sum_i E_{i,j}\cdot A_i}Because we need A_i‘s to compute B_j‘s, and conversely, it is natural to consider some iterative procedure to solve it. Observe that we do not have unicity…

Consider here some starting values for A_i‘s and B_j‘s

A=rep(1,length(levels(df$religion)))
B=rep(1,length(levels(df$occupation)))*sum(N)/sum(E)
A
[1] 1 1 1 1 1
B
[1] 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506 0.6305506

The predicted number of extramarital affairs would be \hat N_{i,j}=E_{i,j}\cdot\hat A_i\cdot \hat B_j

E * A%*%t(B)
        occupation
religion          1          2          3          4          5          6          7
       1  2.5222025  0.6305506  5.0444050  2.5222025 10.0888099  5.6749556  0.0000000
       2 14.5026643  1.8916519  6.9360568 10.7193606 35.3108348 22.6998224  3.7833037
       3 18.2859680  0.6305506  6.3055062  7.5666075 24.5914742 15.7637655  1.2611012
       4 23.9609236  4.4138544  7.5666075 13.2415631 37.2024867 27.7442274  1.2611012
       5  8.1971581  0.6305506  1.8916519  6.3055062 11.9804618 11.9804618  1.8916519
sum(B*E[1,])
[1] 26.48313
sum(B*E[2,])
[1] 95.84369
apply(t(B*t(E)),1,sum)
        1         2         3         4         5 
 26.48313  95.84369  74.40497 115.39076  42.87744 
sum(A*E[,1])
[1] 107
sum(A*E[,2])
[1] 13
apply(A*E,2,sum)
  1   2   3   4   5   6   7 
107  13  44  64 189 133  13

From expressions above, observe that one can very easily write expressions of A_i‘s and B_j‘s as functions of B_j‘s and A_i‘s respectively

A=apply(N,1,sum)/apply(t(B*t(E)),1,sum)
B=apply(N,2,sum)/apply(A*E,2,sum)

Let it iterate one thousand times

for(i in 1:1000){
  A=apply(N,1,sum)/apply(t(B*t(E)),1,sum)
  B=apply(N,2,sum)/apply(A*E,2,sum)
}

We obtain here

A
        1         2         3         4         5 
1.5404346 1.0447195 1.4825650 0.6553159 0.6634763 
B
        1         2         3         4         5         6         7 
0.4685515 0.2629769 0.8454435 0.7245310 0.4889697 0.7770553 1.6753750 
E * A%*%t(B)
        occupation
religion          1          2          3          4          5          6          7
       1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
       2 11.2586111  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
       3 20.1450811  0.3898804 12.5342484 12.8899708 28.2722423 28.8008726  4.9677044
       4 11.6678702  1.2063307  6.6483904  9.9707299 18.9053460 22.4055332  2.1957997
       5  4.0413463  0.1744790  1.6827951  4.8070914  6.1639760  9.7955975  3.3347148

That is our prediction, per category, of the number of affairs. Observe that here, sums per row are equal to observed numbers,

apply(N,1,sum)
  1   2   3   4   5 
 41 103 108  73  30 
apply(E * A%*%t(B),1,sum)
  1   2   3   4   5 
 41 103 108  73  30

as well as sums per colums

apply(N,2,sum)
  1   2   3   4   5   6   7 
 50   3  41  45  94 101  21 
apply(E * A%*%t(B),2,sum)
  1   2   3   4   5   6   7 
 50   3  41  45  94 101  21

Now, why should I mention that here, in the section on the Poisson regression in our course ? Because actually, this is exactly what we get if we run a Poisson regression on those two covariates

reg=glm(y~religion+occupation,data=df,family=poisson)
summary(reg)
Coefficients:
            Estimate Std. Error z value Pr(&gt;|z|)    
(Intercept) -0.32604    0.21325  -1.529 0.126285    
religion2   -0.38832    0.18791  -2.066 0.038783 *  
religion3   -0.03829    0.18585  -0.206 0.836771    
religion4   -0.85470    0.19757  -4.326 1.52e-05 ***
religion5   -0.84233    0.24416  -3.450 0.000561 ***
occupation2 -0.57758    0.59549  -0.970 0.332083    
occupation3  0.59022    0.21349   2.765 0.005699 ** 
occupation4  0.43588    0.20603   2.116 0.034381 *  
occupation5  0.04265    0.17590   0.242 0.808399    
occupation6  0.50587    0.17360   2.914 0.003569 ** 
occupation7  1.27415    0.26298   4.845 1.27e-06 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

First of all, observe that the total sum of predictions equals the total sum of observations

yp = predict(reg,type="response")
sum(yp)
[1] 355
sum(df$y)
[1] 355

But actually, the predicted number of affairs, for our 35 classes, is exactly what we got using Bailey’s technique

xtabs(yp~df$religion+df$occupation)
           df$occupation
df$religion          1          2          3          4          5          6          7
          1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
          2 11.2586112  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
          3 20.1450813  0.3898804 12.5342484 12.8899708 28.2722424 28.8008726  4.9677044
          4 11.6678703  1.2063307  6.6483904  9.9707300 18.9053460 22.4055332  2.1957997
          5  4.0413464  0.1744790  1.6827951  4.8070914  6.1639761  9.7955975  3.3347148
E * A%*%t(B)
        occupation
religion          1          2          3          4          5          6          7
       1  2.8870914  0.4050987 10.4188024  4.4643702 12.0516123 10.7730250  0.0000000
       2 11.2586111  0.8242113  9.7157637 12.8678376 28.6068235 29.2249717 10.5017811
       3 20.1450811  0.3898804 12.5342484 12.8899708 28.2722423 28.8008726  4.9677044
       4 11.6678702  1.2063307  6.6483904  9.9707299 18.9053460 22.4055332  2.1957997
       5  4.0413463  0.1744790  1.6827951  4.8070914  6.1639760  9.7955975  3.3347148

To be more specific, up to a multiplicate constant, series of coefficients are equal here, e.g. for A_i‘s

a=exp(coefficients(reg)[1]+c(0,coefficients(reg)[2:5]))
a/a[1]
          religion2 religion3 religion4 religion5 
1.0000000 0.6781979 0.9624329 0.4254098 0.4307072 
A/A[1]
        1         2         3         4         5 
1.0000000 0.6781979 0.9624329 0.4254098 0.4307072

but also for B_j‘s

b=exp(coefficients(reg)[1]+c(0,coefficients(reg)[6:11]))
b/b[1]
            occupation2 occupation3 occupation4 occupation5 occupation6 occupation7 
  1.0000000   0.5612551   1.8043769   1.5463210   1.0435773   1.6584203   3.5756477 
B/B[1]
        1         2         3         4         5         6         7 
1.0000000 0.5612551 1.8043770 1.5463210 1.0435773 1.6584203 3.5756478

This will have major implications in non-life insurance models (for claims reserving).

Variable explicative dans un intervalle

Suite a une question posée ce matin en cours, je vais faire un rapide billet pour expliquer comment extraire les valeurs inférieures et supérieures quand on a des intervalles, sous R. Commençons par générer des données,

n=200
set.seed(123)
X=rnorm(n)
Y=2+X+rnorm(n,sd = .3)

Supposons maintenant que l’on n’observe plus la vraie variable x mais juste une classe (on va créer huit classes, avec chacune un huitième des observations)

Q=quantile(x = X,(0:8)/8)
Q[1]=Q[1]-.00001
Xcut=cut(X,breaks = Q)
B=data.frame(Y=Y,X=Xcut)

Par exemple, pour la premiere valeur, on a

as.character(Xcut[1])
[1] "(-0.626,-0.348]"

Pour extraire des informations sur ces bornes, on peut utiliser le petit code suivant, qui renvoie la borne inférieure, la borne supérieur, et le milieu de l’intervalle

extraire = function(x){
  ax=as.character(x)
  lower1 = as.numeric( sub("\\((.+),.*", "\\1", ax) )
  lower2 = as.numeric( sub("\\[(.+),.*", "\\1", ax) )
  upper1 = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", ax) )
  upper2 = as.numeric( sub("[^,]*,([^]]*)\\)", "\\1", ax) )
  lower = c(lower1,lower2)
  lower=lower[!is.na(lower)]
  upper = c(upper1,upper2)
  upper=upper[!is.na(upper)]
  mid   = (lower+upper)/2
  return(c(lower=lower,mid=mid,upper=upper))
}

On peut vérifier sur notre première observation

extraire(Xcut[1])
 lower    mid  upper 
-0.626 -0.487 -0.348

Juste pour voir, on peut créer trois variables supplémentaires dans notre base (avec ces trois informations)

B2=Vectorize(function(i) extraire(Xcut[i]))(1:n)
B$lower=B2[1,]
B$mid  =B2[2,]
B$upper=B2[3,]

et on peut comparer 4 régressions (i) on régresse sur nos 8 classes, i.e. nos 8 facteurs (ii) on régresse sur la borne inférieure de l’intervalle, (iii) sur la valeur “moyenne” de l’intervalle (iv) sur la borne supérieure

regF=lm(Y~X,data=B)
regL=lm(Y~lower,data=B)
regM=lm(Y~mid,data=B)
regU=lm(Y~upper,data=B)

On peut comparer les prévisions avec nos quatre modeles

plot(B$Y,predict(regF),ylim=c(0,4))
points(B$Y,predict(regM),col="red")
points(B$Y,predict(regU),col="blue")
points(B$Y,predict(regL),col="purple")
abline(a=0,b=1,lty=2)

Pour aller plus loin, on peut aussi comparer les AIC de nos modèles,

AIC(regF)
[1] 204.5653
AIC(regM)
[1] 201.1201
AIC(regL)
[1] 266.5246
AIC(regU)
[1] 255.0687

Si l’utilisation des bornes inférieures et supérieures n’est pas concluante, ici, on notera qu’utiliser la valeur moyenne de l’intervalle donne des résultats un peu meilleurs que d’utiliser 8 facteurs.