Tag Archives: STT5100

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…

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…(?)