Category Archives: STT5100

Traitement des valeurs manquantes : remplacer les NA par une constante ?

Un rapide billet pour répondre à une question posée à la fin du cours de ce matin (ST5100), par Jean-Pierre Liégeois, jeune lecteur du Var (pour préserver un peu d’anonymat)

Dans mon stage, quand on avait des valeurs manquantes, on me disait de remplacer par -1, puis de rajouter une indicatrice comme quoi la variable vaut -1. Ça permet de ne supprimer ni variable, ni observations. On peut faire ça ?

Si je formalise un peu, on va simuler ici des données, disons x_1 et x_2, on génère ensuite des données suivant un modèle, de la forme y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon. Une proportion \alpha de x_1 sera transformée en NA.  Ce que suggérais Jean-Pierre, c’est de remplacer les valeurs manquantes par -1, puis d’ajuster un modèle y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon. Côté code, c’est assez simple. Par défaut, la stratégie de R est de supprimer les valeurs manquantes. Si 50% des données de x_1 sont manquantes, la moitié des lignes sont supprimées

n=1000
x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.05
indice=sample(1:n,size=round(n*alpha))
base=data.frame(y=y,x1=x1)
base$x1[indice]=NA
reg=lm(y~x1+x2,data=base)

Au lieu de générer un unique échantillon, on va en simuler 10,000, et regarder la distribution de \widehat{\beta}_1,

m=10000
B=rep(NA,m)
for(s in 1:m){
  x1=runif(n)
  x2=runif(n)
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.5
  indice=sample(1:n,size=round(n*alpha))
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=NA
  reg=lm(y~x1+x2,data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Bien sur, avec un taux de valeurs manquantes plus faibles – disons \alpha=5\% – on perd moins d’observations, et donc l’estimateur a une variance plus faible.

Tentons maintenant la stratégie consistant à remplacer les valeurs manquantes par des valeurs numériques fixes, et de rajouter une indicatrice,

B=rep(NA,m)
for(s in 1:m){
  x1=runif(n)
  x2=runif(n)
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.5
  indice=sample(1:n,size=round(n*alpha))
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=-1
  reg=lm(y~x1+x2+I(x1==(-1)),data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Ce qui ne change pas grand chose, on en conviendra…  y compris si le taux de valeurs manquantes passe à 5%,

On peut se demander ce qui se passe si le shift n’est plus de 1 mais de 10 (a priori, c’est arbitraire, la variable x_1  pouvant être plus ou moins dispersée… -1 pour une variable entre 0 et 1, ou entre 0 et 1000, ça n’est pas tout à fait pareil). Mais non, par exemple avec toujours 5% de valeurs manquantes, on a

Si on regarde notre échantillon, en particulier le nuage de points  (x_1,y), on observe

ici, les valeurs manquantes sont choisies au hasard, de manière totalement indépendante,

x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.3333333
indice=sample(1:n,size=round(n*alpha))
clr=rep("black",n)
clr[indice]="red"
plot(x1,y,col=clr)

(ici avec 1/3 de valeurs manquantes, en rouge). Mais on pourrait supposer que les valeurs manquantes sont les plus grandes valeurs de x_1, par exemple,

x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.3333333
indice=sample(1:n,size=round(n*alpha),prob = x1^3)
clr=rep("black",n)
clr[indice]="red"
plot(x1,y,col=clr)

On peut se demander ce que ça donnerait sur l’estimateur \widehat{\beta}_1

Ça ne change pas grand chose, mais on a plus de variance, si on regarde bien. Dernier essai : que se passe-t-il si les variables x_1 et x_2 sont maintenant corrélées,

B=rep(NA,m)
library(mnormt)
r=.8
S = matrix(c(1,r,r,1),2,2)
for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=-1
  reg=lm(y~x1+x2+I(x1==(-1)),data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Cette fois, on a un estimateur biaisé (de l’ordre de 10% sur cet exemple numérique). Manifestement, cette technique n’est pas très concluante…

Je pourrais ajouter que cette méthode ne revient pas à la première, même si la distribution des estimateurs est proche

set.seed(1)
x=rmnorm(n,varcov = S)
x1=pnorm(x[,1])
x2=pnorm(x[,2])
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.2
indice=sample(1:n,size=round(n*alpha),prob = x1^3)
base=data.frame(y=y,x1=x1,x2=x2)
base$x1[indice]=-1
reg1=lm(y~x1+x2+I(x1==(-1)),data=base)
coefficients(reg1)
      (Intercept)                x1                x2 I(x1 == (-1))TRUE 
        1.0988005         1.7454385        -0.5149477         3.1000668 
base$x1[indice]=NA
reg2=lm(y~x1+x2,data=base)
coefficients(reg2)
(Intercept)          x1          x2 
  1.1123953   1.8612882  -0.6548206

Comme je le disais (lors de la discussion qui a suivi le cours) une méthode plus prometteuse est l’imputation. L’idée est de prédire une valeur pour les x_1 manquants. On pourrait être tenté de prendre  \overline{x}_1, la moyenne des x_1 observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de x_1, ici, donc on doit pouvoir faire mieux ! On sait aussi que x_1  et x_2 sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si x_2 est grand, on sait que le x_1 (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, x_1=\alpha_0+\alpha_2x_2+\eta_i, calibré sur les valeurs non-manquantes. Puis on utilise \widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2 pour les valeurs manquantes. C’est simpliste, mais pourquoi pas ? On estime alors le modèle sur cette nouvelle base.

for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
    base$x1[indice]=NA
    reg0=lm(x1~x2,data=base[-indice,])
    base$x1[indice]=predict(reg0,newdata=base[indice,])
  reg=lm(y~x1+x2,data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

sur l’exemple numérique, on obtient

base$x1[indice]=NA
reg0=lm(x1~x2,data=base[-indice,])
base$x1[indice]=predict(reg0,newdata=base[indice,])
reg3=lm(y~x1+x2,data=base)
coefficients(reg3)
(Intercept)          x1          x2 
  1.1593298   1.8612882  -0.6320339

Cette méthode a au moins pu corriger du biais…

Après, si on regarde attentivement, on a exactement la même valeur qu’avec la première méthode qui consiste à supprimer les lignes avec des valeurs manquantes !

summary(reg3)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.15933    0.06649  17.435  < 2e-16 ***
x1           1.86129    0.21967   8.473  < 2e-16 ***
x2          -0.63203    0.20148  -3.137  0.00176 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.051 on 997 degrees of freedom
Multiple R-squared:  0.1094,	Adjusted R-squared:  0.1076 
F-statistic: 61.23 on 2 and 997 DF,  p-value: < 2.2e-16 
 
summary(reg2) 
 
Coefficients: Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.11240    0.06878  16.173  < 2e-16 ***
x1           1.86129    0.21666   8.591  < 2e-16 ***
x2          -0.65482    0.20820  -3.145  0.00172 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.037 on 797 degrees of freedom
  (200 observations deleted due to missingness)
Multiple R-squared:  0.1223,	Adjusted R-squared:   0.12 
F-statistic:  55.5 on 2 and 797 DF,  p-value: < 2.2e-16

Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d’imputation. Par exemple prendre la moyenne des k valeurs de x_1 (observées) pour les k individus ayant des x_2 les plus proches du x_2 de l’individu ayant x_1 manquant :

kpp=function(i,basena,k=5){
  x2=basena$x2[i]
  sb=basena[!is.na(basena$x1),]
  idx=rank(abs(sb$x2-x2))
  mean(sb[which(idx<=k),"x1"])
}

Sur notre base simulée on obtient

base$x1[indice]=NA
base0=base
for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
reg4=lm(y~x1+x2,data=base)
coefficients(reg4)
(Intercept)          x1          x2 
   1.197944    1.804220   -0.806766

Si on regarde ce que ça donne sur nos 10,000 simulations, on a (c’est un peu long, car j’ai codé ça très rapidement, et pas du tout de manière optimale)

for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=NA
  base0=base
    for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
  reg=lm(y~x1+x2,data=base0)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Le biais semble ici plus faible que sans imputation… autrement dit, les méthodes d’imputation me semblent plus robustes que cette stratégie visant à remplacer des NA par une valeur arbitraire, et de rajouter une indicatrice dans la régression.

Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
Davis=data.frame(Y=Davis$weight * 2.204622,
                 X1=Davis$sex)

On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') 
 Y=Davis$Y

car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)

 solve(t(X)%*%X)
            [,1]       [,2]
[1,] 0.008928571 0.00000000
[2,] 0.000000000 0.01136364

et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

 solve(t(X)%*%X) %*% (t(X)%*%Y)
         [,1]
[1,] 125.4272
[2,] 167.3258

ce qui correspond effectivement à la sortie de R,

 reg=lm(Y~0+X1,data=Davis)
 summary(reg)
 
Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1F  125.427      1.960   64.00   <2e-16 ***
X1M  167.326      2.211   75.68   <2e-16 ***

Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes

x=Y[X[,1]==1]
y=Y[X[,2]==1]
nx=length(x)
ny=length(y)

On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M   est le poids moyen des hommes. Et en effet

 mean(y)
[1] 167.3258

C’est finalement très naturel, ou intuitif.

On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici

 sqrt(var(y)/ny)
[1] 2.794391
 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)
[1] 2.794391

car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2

 s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2))
 sqrt(s2/ny)
[1] 2.210863

et en effet, on retombe sur la valeur donnée dans le tableau de régresion

 sqrt(s2/nx)
[1] 1.959721

(pareil pour l’autre coefficient).

On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).

 X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y)
          [,1]
[1,] 125.42724
[2,]  41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)
[1] 41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.427      1.960   64.00   <2e-16 ***
X1M           41.899      2.954   14.18   <2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -6.4475, df = 286, p-value = 4.826e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.62603 -16.30035
sample estimates:
mean of x mean of y 
 143.8626  167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Call:
lm(formula = Y ~ X1, data = Davis)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-57.458 -22.184  -5.512  17.809 118.912 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 143.4382     2.7820   51.56   <2e-16 ***
X1M           0.9645     4.1940    0.23    0.818    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 29.44 on 198 degrees of freedom
Multiple R-squared:  0.000267,	Adjusted R-squared:  -0.004782 
F-statistic: 0.05289 on 1 and 198 DF,  p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y
 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -0.22998, df = 198, p-value = 0.8183
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.235209  7.306165
sample estimates:
mean of x mean of y 
 143.4382  144.4027

et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.

Modélisation “appliquée”

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\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 ?

STT5100, Modèles linéaires appliqués

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.