Depuis la rentrée, je donne un cours de modélisation appliquée, ou pour être plus précis, intitulé “modèles linéaires appliqués“. Ça a été l’occasion de me rendre compte que le mot “appliqué” pouvait prêter à confusion. Dans ma tête, on allait appliquer une théorie (la théorie des modèles linéaires) à cas concrets. Généralement, je présente les grandes lignes de la théorie : les hypothèses, les principaux résultats, et je fais des démonstrations “avec les mains” pour expliquer de manière intuitive certains passages (genre la perte des degrés de liberté dans le test de Student, le signe de la covariance des estimateurs dans le modèle simple, etc) en renvoyant à des livres de référence pour les preuves complètes. Et j’essaye, autant que faire se peut, d’illustrer sur des données, avec du code R (par exemple vendredi dernier, c’était sur la notion de “modalité de référence” quand on régresse sur une variable catégorielle, montrant que tous les modèles sont équivalent et que ce n’est vraiment qu’une histoire de conventions). J’ai aussi donné un gros devoir, avec des vraies données, 2400 observations, 34 variables (quantitatives mais aussi qualitatives). Bref, un vrai cas d’étude. Et j’ai même proposé de revenir sur leurs modélisations en fin de cours…

Et j’ai été surpris quand un groupe d’élèves est venu me voir à la pause, pour me demander s’il y aurait des “chiffres” : si au lieu d’écrire {\displaystyle {\hat {\beta }}_{1}={\frac {\sum x_{i}\sum y_{i}-n\sum x_{i}y_{i}}{\left(\sum x_{i}\right)^{2}-n\sum x_{i}^{2}}}={\frac {\sum (x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum (x_{i}-{\bar {x}})^{2}}}}ou encore {\displaystyle {\hat {\beta }}_{0}={\frac {\sum y_{i}-{\hat {\beta }}_{1}\sum x_{i}}{n}}={\bar {y}}-{\hat {\beta }}_{1}{\bar {x}}}il y aurait des calculs, à la calculatrice, avec des chiffres au lieu des x et des y. J’ai été très supris, parce qu’on est en 2018, que tout le monde parle de “big data” et revenir à la calculatrice pour calculer des sommes de 10 nombres, c’est possible. En tous cas, certains profs le font encore. J’ai dit que si ça pouvait les rassurer, je pouvais chercher des bases dans des vieux livres, mais que c’était pas ma manière de comprendre le terme “appliqué” dans le nom du cours. Je me suis souvenu du livre de Robert Pindyck et Daniel Rubinfeld, Econometric Models and Economic Forecasts, qui était le livre de référence de la SOA il y a 20 ans, qui avait des tableaux de ce genre. Mais je ne l’ai pas trouvé en ligne. Par contre, je suis tombé sur la Table 4.3. dans le livre de John Neter, Applied Linear Regression Models, paru en 1974 (repris dans la Table 3.4. de Applied Linear Statistical Models)

On a ici juste 11 observations, un modèle linéaire simple. Il est possible de faire tous les calculs à la main. Mais je me suis dit que ce n’était pas un cadeau : on ne peut pas faire grand chose sans faire un minimum de maths… Faisons juste une régression linéaire sur ces données, pour “voir”

base = data.frame(x=c(125,100,200,75,150,175,75,175,125,200,100),
                  y=c(160,112,124,28,152,156,42,124,150,104,136))
plot(base,xlab="Size of Minumum Deposit",
     ylab="Number of New Accounts",pch=19,ylim=c(30,180))
abline(lm(y~x,data=base),col="red")

La droite de régression (qui minimise la somme des carrés des erreurs) est la courbe (ou ici “droite”, on est sur un modèle linéaire) rouge

summary(lm(y~x,data=base))
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  50.7225    39.3979   1.287     0.23
x             0.4867     0.2747   1.772     0.11

On ne voit rien, littéralement ! Surtout que sur cet exemple, la question qui nous brule les lèvres, c’est “où est le maximum ?” car on voit bien qu’on a une courbe qui semble concave, croissante au début, décroissante sur la fin… On peut faire du lissage non-paramétrique pour s’en convaincre

with(base, 
     scatter.smooth(x, y,
        lpars = list(col = "red"),
        xlab="Size of Minumum Deposit",
       ylab="Number of New Accounts",
       ylim=c(30,180),pch=19))

Peut-aller plus loin, pour répondre à la question “où est le maximum ?”. En particulier, peut-on suggérer une valeur, voir trouver un intervalle de confiance ?

On peut continuer avec un modèle paramétrique – l’essentiel de ce qu’on voit dans les cours de statistique tourne autour des modèles paramétrique. Mais si on est plus pragmatique, on peut considérer un modèle quadratique, autrement dit, notre prévision sera une parabole. Ça pourrait être cohérent avec ce qu’on voit quand on trace les points, et puis souvent, les paraboles ont un maximum

with(base, 
     scatter.smooth(x, y,
        lpars = list(col = "red",lty=2),
        xlab="Size of Minumum Deposit",
        ylab="Number of New Accounts",
        ylim=c(30,180),pch=19))
reg=lm(y~x+I(x^2),data=base)
vx=seq(min(base$x),max(base$x),length=251)
vy=predict(reg,newdata=data.frame(x=vx))
lines(vx,vy,col="red")

On peut voir que le modèle pourrait sembler réaliste, non seulement visuellement, mais en plus si on regarde la sortie de la régression (même si on pourrait discuter ici l’interprétation de la significativité…)

summary(reg)
 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.255e+02  5.589e+01  -5.824 0.000394 ***
x            6.569e+00  8.744e-01   7.513 6.84e-05 ***
I(x^2)      -2.203e-02  3.143e-03  -7.011 0.000111 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Or on peut montrer que pour une parabole de la forme y=\beta_2x^2+\beta_1x+\beta_0 l’optimum (qui sera un maximum si \beta_2<0 – la courbe part en -\infty à gauche et à droite) sera obtenu en x^\star=\theta=-\beta_1/2\beta_2. Un estimateur naturel de \theta est alors \widehat{\theta}=-\widehat{\beta}_1/2\widehat{\beta}_2 où les estimateurs sont obtenus ici par minimisation de la somme des carrés des erreurs. Mais comment avoir la variance de cet estimateur ? On peut naturellement tenter la \Delta-méthode (décrite sur wikipedia), en considérant \displaystyle{\theta=g(\beta_1,\beta_2)=\frac{-\beta_1}{2\beta_2}}Comme \displaystyle{\frac{\partial \theta}{\partial \beta_1}=\frac{-1}{2\beta_2}}et \displaystyle{\frac{\partial \theta}{\partial \beta_2}=\frac{\beta_1}{2\beta_2^2}}la variance (asymptotique) vaut ici \left[ \frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right] \left[ \begin{array}{cc} \sigma_{1}^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{array} \right] \left[\frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right]^{\text{ T}}soit \frac{\sigma_1^2\beta_2^2-2\beta_1\beta_2\sigma_{12}+\beta_1^2\sigma_2^2}{4\beta_2^2}On obtient alors un estimateur de cette variance asymptotique en remplaçant les valeurs inconnues par leurs estimateurs. Ici, on aurait, numériquement

beta=coefficients(reg)
n=nrow(base)
sigma=n*vcov(reg)[2:3,2:3]
dg = c(-.5/beta[3],.5*beta[2]/beta[3]^2)
theta=-beta[2]/(2*beta[3])
theta
[1]  149.0676 
s2=t(dg) %*% sigma %*% dg
s2
         [,1]
[1,] 94.59911
sqrt(s2)
         [,1]
[1,] 9.726207

Autrement dit, si on suppose la normalité de notre estimateur (ce que nous dit d’ailleurs la \Delta-méthode, tout du moins asymptotiquement), on a l’intervalle de confiance suivant

theta+qt(c(.025,.975),n-3)*sqrt(s2)
[1] 126.6389 171.4963
vx=theta
vy=predict(reg,newdata=data.frame(x=vx))
points(vx,vy,pch=19,cex=1.5,col="red")
arrows(vx-qt(.975,n-3)*sqrt(s2),vy,
       vx+qt(.975,n-3)*sqrt(s2),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

On peut aussi tenter une autre stratégie, basée sur l’utilisation de la vraisemblance profilée. Je renvoie à mon précédent billet pour expliquer l’idée de la méthode, malheureusement pas assez enseignée dans les cours de statistique… Pour notre modèle quadratique, on maximisait la log-vraisemblance, classiquement, sous une hypothèse de résidus Gaussiens et homoscédastiques

logL = function(pm){
  b0=pm[1]
  b1=pm[2]
  b2=pm[3]
  b3=pm[4]
  -sum(log(dnorm( (base$y-(b0+b1*base$x+b2*base$x^2)) ,0,b3)))
}

Ici, la première astuce est d’introduire \theta (là où se trouve le maximum de la parabole) comme un des paramètres du modèle – par exemple à la place de \beta_2

logL = function(pm){
  b0=pm[1]
  b1=pm[2]
  theta=pm[3]
  b3=pm[4]
  -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3)))
}

Si on cherche juste le maximum de la vraisemblance, on obtient

 optim(par=c(-323,6.5,140,5),logL)
$par
[1] -325.521244    6.569328  149.070597   13.693488
 
$value
[1] 44.39629

ce qui est cohérent avec notre calcul précédant (et c’est normal, c’est la propriété du cours de statistique qui dit qu’on peut reparamétrer le modèle, ça ne change rien). La seconde astuce est celle de la vraisemblance profilée : on dit que dans notre paramètre multivarié, un est plus important que les autres. Et les autres sont juste maximisés, à \theta donné. Techniquement, la log-vraisemblance profilée du paramètre \theta est

logL=function(theta){
  f=function(pm){
  b0=pm[1]
  b1=pm[2]
  b3=pm[3]
  -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3)))
  }
  optim(par=c(-323,6.5,5),f)$value
}

On peut d’ailleurs la tracer

v1 = seq(130,190)
v2 = Vectorize(logL)(v1)
plot(v1,-v2,type="l",xlim=range(base$x),xlab="",ylab="log vraisemblance")

Le maximum est atteint, ici en

opt
$minimum
[1] 149.068

Ce qui est cohérent avec nos calculs. On peut ensuite un résultat sur le test du rapport de vraisemblance, qui nous dit que la différence entre la log-vraisemblance au maximum et en la vraie valeur suit la moitié du quantile d’une loi du chi-deux à 1 degré de liberté

points(opt$minimum,-opt$objective,pch=19,cex=1.5,col="red")
ref=-opt$objective-.5*qchisq(.95,1)
abline(h=ref,lty=2,col="red")

On peut alors alors tracer l’intervalle de confiance sur le graphique initial

with(base, 
     scatter.smooth(x, y,
                    lpars = list(col = "red",lty=2),
                    xlab="Size of Minumum Deposit",
                    ylab="Number of New Accounts",
                    ylim=c(30,180),pch=19))
reg=lm(y~x+I(x^2),data=base)
vx=seq(min(base$x),max(base$x),length=251)
vy=predict(reg,newdata=data.frame(x=vx))
lines(vx,vy,col="red")
vx=opt$minimum
vy=predict(reg,newdata=data.frame(x=vx))
points(vx,vy,pch=19,cex=1.5,col="red")
arrows(min(v1[id]),vy,max(v1[id]),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

Voilà au moins deux méthodes, basées sur deux résultats du cours de statistique. Après, comme la majorité des techniques de statistique, ce sont des résultats asymptotiques, et rien ne garantie leur validité avec seulement 11 observations (“rien” est un peu fort… on pourrait faire des développements sur ce point).

Une autre solution est d’utiliser des simulations. On suppose que les observations, c’est un modèle, et du bruit. On peut prendre comme modèle le modèle non-paramétrique (lissage local) et supposer le bruit Gaussien. Pour générer d’autres échantillons, on va garder nos observations en x, par contre, pour y, on va utiliser \widehat{y}+\varepsilon où \varepsilon sera tiré suivant une loi normale

reg=loess.smooth(x = base$x, y= base$y, evaluation = 501)
pred=function(x) reg$y[which.min(abs(x-reg$x))]
ypred=Vectorize(pred)(base$x)
epsilon = base$y-ypred
simu = function(graph=FALSE){
newbase = data.frame(x= base$x,
                     y= ypred + rnorm(n,mean(epsilon),sd(epsilon)))
reg=loess.smooth(x = newbase$x, y= newbase$y, evaluation = 501)
lines(reg$x,reg$y,col=rgb(1,0,0,.2))
reg$x[which.max(reg$y)]}
for(i in 1:20) simu(TRUE)
lines(loess.smooth(x = base$x, y= base$y, evaluation = 501),lwd=1.2,col="red")

On récupère ainsi plein de modèles, estimés sur ces nouvelles données. Compte tenu de l’asymétrie des données, on va oublier le modèle quadratique, et utiliser un modèle non-paramétrique, là encore. Et on calcule numériquement le maximum. Et on répète 10,000 fois.

V=Vectorize(simu)(1:1e5)
hist(V,probability = TRUE, col="light blue", border="white",ylim=c(0,.07))
lines(density(V),col="blue")
quantile(V,c(.025,.975))

On a ici la distribution empirique du maximum, observée sur nos 10,000 échantillons simulés. On peut même obtenir un intervalle de confiance, en prenant les quantiles empiriques

vy=pred(mean(V))
arrows(quantile(V,.025),vy,quantile(V,.975),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

Bref, j’ai l’impression que pour faire des applications, il vaut mieux bien connaître la théorie… non ?

La semaine prochaine commencera le premier cours sur les modèles linéaires appliqués, STT5100. Pas de notes de cours ou de slides/powerpoint, a priori je ferais le cours au tableau, en interagissant autant que possible avec R. Je mettrais sur github les codes et les bases de données qui seront utilises pour illustrer. Il y a aussi des examens des années passées (avec les corrections). Le cours aura lieu les vendredi matin, salle SH-3140. Le blog (via le tag STT5100) permettra d’apporter certains compléments, ici ou la, a priori sur des points techniques.