1 Stima parametri modello lineare

1.1 Inizializzazione del codice necessario

Vengono caricati i packages necessari per realizzare questo documento

1.2 Due variabili esplicative

Estraiamo alcune variabili dal data.set antropometric e usiamo per questi esempi elementari sui modelli lineari la variabile peso come variabile dipendente.

source 
data(antropometric)

x1=antropometric$ALTEZZA
x2=antropometric$TORACECM
y=antropometric$PESOKG
view(dfSummary(as.data.frame(cbind(x1,x2,y))),method = "render")

Data Frame Summary

Dimensions: 1427 x 3
Duplicates: 146
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 x1 [numeric]
Mean (sd) : 151.9 (10.1)
min ≤ med ≤ max:
127 ≤ 151 ≤ 183
IQR (CV) : 15 (0.1)
 54 distinct values 1427 (100.0%) 0 (0.0%)
2 x2 [numeric]
Mean (sd) : 75.6 (7.8)
min ≤ med ≤ max:
57 ≤ 74 ≤ 104
IQR (CV) : 10 (0.1)
 44 distinct values 1427 (100.0%) 0 (0.0%)
3 y [numeric]
Mean (sd) : 45 (10.7)
min ≤ med ≤ max:
21 ≤ 43 ≤ 100
IQR (CV) : 14 (0.2)
 65 distinct values 1427 (100.0%) 0 (0.0%)

Generated by summarytools 1.0.0 (R version 4.1.2)
2022-02-25

source 
## oppure summary(as.data.frame(cbind(x1,x2,y)))
plot(x1,y)
plot

source 
plot(x2,y)
plot

source 
# DIPENDENZA DI Y DA X1 E X2

MLA.regression3d(cbind(x1,x2,y),ell=FALSE,conf.exp = FALSE)

Potremmo considerare anche una terza variabile nominale, le cui modalità indico nel grafico con colori diversi.

source 
rgl.close()
x3=as.numeric(antropometric$HABITUS)
MLA.regression3d(cbind(x1,x2,y),level=0.9, ext=0.1,ell=FALSE,conf.exp = FALSE,col.obs=x3)

L’analisi dettagliata dei modelli lineari con variabili esplicative sia quantitative che qualitative verrà fatta in seguito. E’ evidente che potremo porci il problema di costruire delle relazioni diverse per ciascuna modalità o una unica. Ma ce ne occuperemo più in là

1.3 Stima dei parametri del modello

Abbiamo visto che lo stimatore dei minimi quadrati è \[ \mathbf{b}=(X^T X)^{-1}X^T \mathbf{y} \] Applichiamo la formula diretta a un piccolo esempio sullo stesso dataset:

source 
X=as.data.frame(scale(antropometric[,c(7,9:13)]))
Y=as.matrix(antropometric[,8])
n=nrow(X)
ones=array(1,n)
X=as.matrix(cbind(ones,X))
round(t(X)%*%X/(n),2)
output 
         ones ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
ones        1    0.00     0.00     0.00     0.00     0.00   0.00
ALTEZZA     0    1.00     0.59     0.48     0.75     0.74   0.89
TORACECM    0    0.59     1.00     0.48     0.69     0.77   0.58
CRANIOCM    0    0.48     0.48     1.00     0.50     0.49   0.47
BISACROM    0    0.75     0.69     0.50     1.00     0.76   0.78
BITROCAN    0    0.74     0.77     0.49     0.76     1.00   0.71
SPANCM      0    0.89     0.58     0.47     0.78     0.71   1.00
source 
k=ncol(X)

b=solve(t(X)%*%X)%*%t(X)%*%Y
print(b)
output 
               [,1]
ones     44.9810792
ALTEZZA   2.1561275
TORACECM  6.7819084
CRANIOCM  0.5103593
BISACROM  0.1232399
BITROCAN  1.9288991
SPANCM   -0.1548494
source 
ys=X%*%as.matrix(b)
rgl.close()

Calcoliamo anche l’inversa dal momento che \(V(\mathbf{b})=(X^T X)^{-1} \sigma^2\):

source 
round(solve(t(X)%*%X/n),2)
output 
         ones ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
ones        1    0.00     0.00     0.00     0.00     0.00   0.00
ALTEZZA     0    5.58     0.18    -0.17    -0.26    -1.17  -3.96
TORACECM    0    0.18     2.71    -0.31    -0.74    -1.59   0.12
CRANIOCM    0   -0.17    -0.31     1.44    -0.22    -0.08  -0.12
BISACROM    0   -0.26    -0.74    -0.22     3.58    -0.87  -1.40
BITROCAN    0   -1.17    -1.59    -0.08    -0.87     3.81  -0.02
SPANCM      0   -3.96     0.12    -0.12    -1.40    -0.02   5.62

1.4 Modelli lineari con R: lm()

Fortunatamente non dobbiamo impostare calcoli matriciali esplicitamente con R, ma utilizziamo la funzione lm()

source 
X=as.data.frame(X)
l1 =lm(Y~ALTEZZA+TORACECM+CRANIOCM+BISACROM+BITROCAN+SPANCM,data=X)

summary(l1)
output 

Call:
lm(formula = Y ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + 
    SPANCM, data = X)

Residuals:
     Min       1Q   Median       3Q      Max 
-28.6479  -1.8852  -0.0398   1.7615  24.9824 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 44.98108    0.09139 492.183  < 2e-16 ***
ALTEZZA      2.15613    0.21585   9.989  < 2e-16 ***
TORACECM     6.78191    0.15045  45.078  < 2e-16 ***
CRANIOCM     0.51036    0.10956   4.658 3.49e-06 ***
BISACROM     0.12324    0.17303   0.712    0.476    
BITROCAN     1.92890    0.17837  10.814  < 2e-16 ***
SPANCM      -0.15485    0.21673  -0.714    0.475    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.452 on 1420 degrees of freedom
Multiple R-squared:  0.8968,    Adjusted R-squared:  0.8964 
F-statistic:  2056 on 6 and 1420 DF,  p-value: < 2.2e-16
source 
print(b)
output 
               [,1]
ones     44.9810792
ALTEZZA   2.1561275
TORACECM  6.7819084
CRANIOCM  0.5103593
BISACROM  0.1232399
BITROCAN  1.9288991
SPANCM   -0.1548494

1.5 Stima di \(\sigma^2\)

\[\hat{\sigma}^2=\frac{\sum (y_i-\hat{y}_i)^2}{n-k}\]
source 
s2=sum((Y-ys)^2)/(n-k)
print(s2)
output 
[1] 11.91875
source 
print(sqrt(s2))
output 
[1] 3.452355

1.6 L’istruzione lm()

Ovviamente in generale non è necessario ritrasformare i dati, definire una matrice X, un vettore Y, etc. (anche se consiglio vivamente lo studente di Statistica di farlo, almeno le prime volte che usa modelli lineari con più variabili esplicative, R è molto adatto al calcolo matriciale), ma possiamo direttamente usare l’istruzione lm(var.dipendente~var1+var2+...+vark,data=data.frame) utilizzando per le variabili e il data.frame i nomi originali, cosa che rende poi più leggibili i risultati (ho messo nell’istruzione i dati in forma standardizzata, con scale(). Se lasciate i dati originari le stime dei parametri saranno ovviamente diverse)

Ricordo che traslare i dati (ossia p.e. lavorare con gli scarti dalle medie) ha un effetto solo sul termine noto, mentre standardizzare o comunque cambiare scala ha un effetto in generale su tutti i coefficienti

source 
l1 =lm(PESOKG~ALTEZZA+TORACECM+CRANIOCM+BISACROM+BITROCAN+SPANCM,data=as.data.frame(scale(antropometric[,7:13])))

summary(l1)
output 

Call:
lm(formula = PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + 
    BITROCAN + SPANCM, data = as.data.frame(scale(antropometric[, 
    7:13])))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.67154 -0.17580 -0.00371  0.16427  2.32972 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.675e-15  8.523e-03   0.000    1.000    
ALTEZZA      2.011e-01  2.013e-02   9.989  < 2e-16 ***
TORACECM     6.324e-01  1.403e-02  45.078  < 2e-16 ***
CRANIOCM     4.759e-02  1.022e-02   4.658 3.49e-06 ***
BISACROM     1.149e-02  1.614e-02   0.712    0.476    
BITROCAN     1.799e-01  1.663e-02  10.814  < 2e-16 ***
SPANCM      -1.444e-02  2.021e-02  -0.714    0.475    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3219 on 1420 degrees of freedom
Multiple R-squared:  0.8968,    Adjusted R-squared:  0.8964 
F-statistic:  2056 on 6 and 1420 DF,  p-value: < 2.2e-16

consiglio poi di leggere le istruzioni del comando lm() che sarà una delle istruzioni che userete di più durante i vostri studi, e un esempio con str() di lista di oggetti restituiti dalla funzione lm() (che è un oggetto di classe lm ). Ne useremo solo alcuni

source 
help(lm)
str(l1)
output 
List of 12
 $ coefficients : Named num [1:7] -1.68e-15 2.01e-01 6.32e-01 4.76e-02 1.15e-02 ...
  ..- attr(*, "names")= chr [1:7] "(Intercept)" "ALTEZZA" "TORACECM" "CRANIOCM" ...
 $ residuals    : Named num [1:1427] 0.3042 0.152 -0.0237 0.1818 -0.0908 ...
  ..- attr(*, "names")= chr [1:1427] "1" "2" "3" "4" ...
 $ effects      : Named num [1:1427] -4.66e-15 2.73e+01 2.27e+01 -1.77 1.07 ...
  ..- attr(*, "names")= chr [1:1427] "(Intercept)" "ALTEZZA" "TORACECM" "CRANIOCM" ...
 $ rank         : int 7
 $ fitted.values: Named num [1:1427] -1.142 -1.176 -0.907 -0.46 -0.187 ...
  ..- attr(*, "names")= chr [1:1427] "1" "2" "3" "4" ...
 $ assign       : int [1:7] 0 1 2 3 4 5 6
 $ qr           :List of 5
  ..$ qr   : num [1:1427, 1:7] -37.7757 0.0265 0.0265 0.0265 0.0265 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:1427] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:7] "(Intercept)" "ALTEZZA" "TORACECM" "CRANIOCM" ...
  .. ..- attr(*, "assign")= int [1:7] 0 1 2 3 4 5 6
  ..$ qraux: num [1:7] 1.03 1.02 1.01 1 1.02 ...
  ..$ pivot: int [1:7] 1 2 3 4 5 6 7
  ..$ tol  : num 1e-07
  ..$ rank : int 7
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 1420
 $ xlevels      : Named list()
 $ call         : language lm(formula = PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + SPANCM, data = as.data.frame(scale(an| __truncated__
 $ terms        :Classes 'terms', 'formula'  language PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + SPANCM
  .. ..- attr(*, "variables")= language list(PESOKG, ALTEZZA, TORACECM, CRANIOCM, BISACROM, BITROCAN, SPANCM)
  .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:7] "PESOKG" "ALTEZZA" "TORACECM" "CRANIOCM" ...
  .. .. .. ..$ : chr [1:6] "ALTEZZA" "TORACECM" "CRANIOCM" "BISACROM" ...
  .. ..- attr(*, "term.labels")= chr [1:6] "ALTEZZA" "TORACECM" "CRANIOCM" "BISACROM" ...
  .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(PESOKG, ALTEZZA, TORACECM, CRANIOCM, BISACROM, BITROCAN, SPANCM)
  .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
  .. .. ..- attr(*, "names")= chr [1:7] "PESOKG" "ALTEZZA" "TORACECM" "CRANIOCM" ...
 $ model        :'data.frame':  1427 obs. of  7 variables:
  ..$ PESOKG  : num [1:1427] -0.838 -1.024 -0.931 -0.278 -0.278 ...
  ..$ ALTEZZA : num [1:1427] -0.886 -0.787 -0.986 -1.481 -0.787 ...
  ..$ TORACECM: num [1:1427] -1.105 -1.2331 -0.849 -0.2088 -0.0808 ...
  ..$ CRANIOCM: num [1:1427] -1.086 -0.47 -0.47 -0.47 0.762 ...
  ..$ BISACROM: num [1:1427] -1.521 -0.512 -1.521 -0.848 -0.848 ...
  ..$ BITROCAN: num [1:1427] -1.199 -1.199 -0.84 -0.122 -0.122 ...
  ..$ SPANCM  : num [1:1427] -1.394 -0.411 -1.305 -1.662 -1.215 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + SPANCM
  .. .. ..- attr(*, "variables")= language list(PESOKG, ALTEZZA, TORACECM, CRANIOCM, BISACROM, BITROCAN, SPANCM)
  .. .. ..- attr(*, "factors")= int [1:7, 1:6] 0 1 0 0 0 0 0 0 0 1 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:7] "PESOKG" "ALTEZZA" "TORACECM" "CRANIOCM" ...
  .. .. .. .. ..$ : chr [1:6] "ALTEZZA" "TORACECM" "CRANIOCM" "BISACROM" ...
  .. .. ..- attr(*, "term.labels")= chr [1:6] "ALTEZZA" "TORACECM" "CRANIOCM" "BISACROM" ...
  .. .. ..- attr(*, "order")= int [1:6] 1 1 1 1 1 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(PESOKG, ALTEZZA, TORACECM, CRANIOCM, BISACROM, BITROCAN, SPANCM)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:7] "numeric" "numeric" "numeric" "numeric" ...
  .. .. .. ..- attr(*, "names")= chr [1:7] "PESOKG" "ALTEZZA" "TORACECM" "CRANIOCM" ...
 - attr(*, "class")= chr "lm"

Abbiamo visto il test F per la saggiare l’ipotesi \(H_0: \mathbf{\beta}=\mathbf{0}\) contro un alternativa generica (in effetti l’ipotesi nulla riguarda tutti i parametri tranne il termine noto, ma abbiamo visto come sia possibile rendere lo stimatore del termine noto indipendente da quelli relativi alle variabili esplicative; è come dire che l’ipotesi nulla è semplicemente \(H_0: E(\mathbf{Y}_i)=\mu, \forall i\))

source 
vtot=sum((Y-mean(Y))^2)
vspieg=sum((ys-mean(ys))^2)
vres=sum((Y-ys)^2)

s2=vres/1420
s2
output 
[1] 11.91875
source 
F=(vspieg/6)/(vres/1420)
print(F)
output 
[1] 2056.31

Il valore di F ovviamente coincide con quello fornito da lm().

1.7 Saggiare un ipotesi generica

Ci poniamo ora un altro problema: saggiare invece un ipotesi che non riguarda tutti i parametri, ma solo alcuni ad esempio: \(H_0: \beta_1=\beta_2=0\) (sebbene non sia particolarmente utile per questi dati, ma si tratta solo di un esempio numerico, senza una reale valenza applicativa).

E’ come dire che l’ipotesi nulla prevede un modello in cui mancano \(X_1, X_2\), nel nostro piccolo esempio ALTEZZA, TORACECM e che prevede 4 variabili esplicative invece di sei.

Che fare? Sappiamo già dalla teoria sviluppata, che la verosimiglianza dipende da \(R(\mathbf{\beta)}|H_0\), sotto l’ipotesi nulla e da \(R(\mathbf{\beta)}|H_1\) sotto l’ipotesi alternativa (che qui comprende l’intero spazio dei parametri \(\Omega\))

Occorre quindi adattare il modello lineare specificato da \(H_0\)
source 
l0 =lm(PESOKG~CRANIOCM+BISACROM+BITROCAN+SPANCM,
       data=as.data.frame(scale(antropometric[,7:13])))

summary(l0)
output 

Call:
lm(formula = PESOKG ~ CRANIOCM + BISACROM + BITROCAN + SPANCM, 
    data = as.data.frame(scale(antropometric[, 7:13])))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.37123 -0.28379 -0.02104  0.25310  2.86101 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 9.015e-16  1.340e-02   0.000  1.00000    
CRANIOCM    1.238e-01  1.584e-02   7.817 1.05e-14 ***
BISACROM    1.913e-01  2.460e-02   7.777 1.41e-14 ***
BITROCAN    5.804e-01  2.192e-02  26.479  < 2e-16 ***
SPANCM      7.057e-02  2.236e-02   3.157  0.00163 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5061 on 1422 degrees of freedom
Multiple R-squared:  0.7446,    Adjusted R-squared:  0.7439 
F-statistic:  1036 on 4 and 1422 DF,  p-value: < 2.2e-16

Confrontiamo ora il modello l1 col modello l0 (cioè confrontiamo i modelli stimati sotto H1 e sotto H0)

source 
my=0
spieg0=sum((l0$fitted.values-my)^2)
res0=sum(l0$residuals^2)
spieg1=sum((l1$fitted.values-my)^2)
res1=sum(l1$residuals^2)

print(c(spieg0,res0))
output 
[1] 1061.7619  364.2381
source 
print(c(spieg1,res1))
output 
[1] 1278.8172  147.1828

Come saggiare H0 adesso? Come dimostreremo, il test del rapporto delle verosimiglianze è basato sul confronto fra res0 e res1 (è ovvio che si avrà sempre res0 > res1); se questo incremento di devianza, dovuto ai vincoli imposti da H0, è grande avremo evidenza contro H0!

In termini di devianza spiegata, ovviamente l1 spiega più di l0…

Come giudicare dunque l’incremento di devianza residua res0-res1= 364.2381041 - 147.1828031?

** Ancora con un test F!!! **
source 
F=((res0-res1)/2)/(res1/1420)
print(F)
output 
[1] 1047.06

Potevamo ottenere lo stesso risultato con

source 
anova(l0,l1)
output 
Analysis of Variance Table

Model 1: PESOKG ~ CRANIOCM + BISACROM + BITROCAN + SPANCM
Model 2: PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + 
    SPANCM
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1   1422 364.24                                  
2   1420 147.18  2    217.06 1047.1 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adesso vediamo nella teoria il collegamento esplicito col rapporto delle verosimiglianze, il rapporto F che costruiamo, e la verifica di ipotesi generali nel modello lineare.

1.8 Il coefficiente di determinazione multipla

Il coeffciente di determinazione multipla indica la parte di devianza (o varianza) di y spiegata dalle X e varia ovviamente fra zero e uno. Vediamo una relazione interessante con la matrice di varianza e covarianza unifica di y e delle X (prima in r e poi spiegazione)

source 
xfull=as.data.frame(scale(antropometric[,c(8,7,9:13)]))
source 
print(names(xfull))
output 
[1] "PESOKG"   "ALTEZZA"  "TORACECM" "CRANIOCM" "BISACROM" "BITROCAN" "SPANCM"
source 
print(round(cov(xfull),2))
output 
         PESOKG ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
PESOKG     1.00    0.72     0.91     0.54     0.75     0.84   0.69
ALTEZZA    0.72    1.00     0.59     0.48     0.75     0.74   0.89
TORACECM   0.91    0.59     1.00     0.48     0.69     0.77   0.58
CRANIOCM   0.54    0.48     0.48     1.00     0.50     0.49   0.47
BISACROM   0.75    0.75     0.69     0.50     1.00     0.76   0.78
BITROCAN   0.84    0.74     0.77     0.49     0.76     1.00   0.71
SPANCM     0.69    0.89     0.58     0.47     0.78     0.71   1.00
source 
print(round(cor(xfull),2))
output 
         PESOKG ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
PESOKG     1.00    0.72     0.91     0.54     0.75     0.84   0.69
ALTEZZA    0.72    1.00     0.59     0.48     0.75     0.74   0.89
TORACECM   0.91    0.59     1.00     0.48     0.69     0.77   0.58
CRANIOCM   0.54    0.48     0.48     1.00     0.50     0.49   0.47
BISACROM   0.75    0.75     0.69     0.50     1.00     0.76   0.78
BITROCAN   0.84    0.74     0.77     0.49     0.76     1.00   0.71
SPANCM     0.69    0.89     0.58     0.47     0.78     0.71   1.00
source 
cyy=solve(cor(xfull))[1,1]
print(cyy)
output 
[1] 9.688632
source 
print(1-1/cyy)
output 
[1] 0.8967863

Possiamo vedere che quest’ultimo valore è proprio il coefficiente di determinazione multipla (modello l1 con 6 variabili esplicative):

source 
inv.xfull=solve(cor(xfull))
print(round(inv.xfull,2))
output 
         PESOKG ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
PESOKG     9.69   -1.95    -6.13    -0.46    -0.11    -1.74   0.14
ALTEZZA   -1.95    5.97     1.41    -0.08    -0.24    -0.82  -3.99
TORACECM  -6.13    1.41     6.58    -0.02    -0.67    -0.48   0.03
CRANIOCM  -0.46   -0.08    -0.02     1.46    -0.21     0.00  -0.12
BISACROM  -0.11   -0.24    -0.67    -0.21     3.58    -0.85  -1.40
BITROCAN  -1.74   -0.82    -0.48     0.00    -0.85     4.12  -0.04
SPANCM     0.14   -3.99     0.03    -0.12    -1.40    -0.04   5.62
source 
print(inv.xfull[1,1])
output 
[1] 9.688632
source 
print(1-1/cyy)
output 
[1] 0.8967863
source 
print(spieg1/(spieg1+res1))
output 
[1] 0.8967863

Si dimostra utilizzando le proprietà delle inverse di matrici orlate e sfruttando il fatto che: \[1-R^2_{y.1,2,...,k}=\frac{y'(I-H)y}{y'y}\] avendo ipotizzato y a media nulla (e quindi le y sono scartate dalla propria media)

e si dimostra che \(\frac{1}{c_{yy}}=\sigma_y(1-R^2_{y.1,2,...,k})\) (il termine \(\sigma_y\) si può omettere se si sta lavorando con variabili standardizzate e quindi con matrici di correlazione )

questo risultato vale in generale per qualsiasi matrice di correlazione: \(\frac{1}{c_{jj}}=\sigma_y(1-R^2_{j})\)

1.9 Test su singoli coefficienti di regressione

(questi esempi hanno servono per ora solo come esercizi su dati già noti: non avrebbe senso effettuare diversi test con diverse H0 sugli stessi dati, perlomeno non la teoria che stiamo studiando finora)

Riprendiamo il modello adattato l1, e confrontiamolo ora con un modello in cui manca una sola variabile, per esempio la variabile numero 2, ossia la variabile TORACECM (come se volessimo saggiare l’ipotesi \(H_0: \beta_2=0\), contro l’alternativa \(H_0: \beta_2 \ne 0\))

source 
l1 =lm(PESOKG~ALTEZZA+TORACECM+CRANIOCM+BISACROM+BITROCAN+SPANCM,data=as.data.frame(scale(antropometric[,7:13])))

l1.senzaX2 =lm(PESOKG~ALTEZZA+CRANIOCM+BISACROM+BITROCAN+SPANCM,data=as.data.frame(scale(antropometric[,7:13])))

anova(l1.senzaX2,l1)
output 
Analysis of Variance Table

Model 1: PESOKG ~ ALTEZZA + CRANIOCM + BISACROM + BITROCAN + SPANCM
Model 2: PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + BITROCAN + 
    SPANCM
  Res.Df    RSS Df Sum of Sq    F    Pr(>F)    
1   1421 357.80                                
2   1420 147.18  1    210.62 2032 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
source 
summary(l1)
output 

Call:
lm(formula = PESOKG ~ ALTEZZA + TORACECM + CRANIOCM + BISACROM + 
    BITROCAN + SPANCM, data = as.data.frame(scale(antropometric[, 
    7:13])))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.67154 -0.17580 -0.00371  0.16427  2.32972 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.675e-15  8.523e-03   0.000    1.000    
ALTEZZA      2.011e-01  2.013e-02   9.989  < 2e-16 ***
TORACECM     6.324e-01  1.403e-02  45.078  < 2e-16 ***
CRANIOCM     4.759e-02  1.022e-02   4.658 3.49e-06 ***
BISACROM     1.149e-02  1.614e-02   0.712    0.476    
BITROCAN     1.799e-01  1.663e-02  10.814  < 2e-16 ***
SPANCM      -1.444e-02  2.021e-02  -0.714    0.475    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3219 on 1420 degrees of freedom
Multiple R-squared:  0.8968,    Adjusted R-squared:  0.8964 
F-statistic:  2056 on 6 and 1420 DF,  p-value: < 2.2e-16
source 
45.078^2
output 
[1] 2032.026

Il test F per la verifica di un ipotesi nulla sul singolo coefficiente di regressione \(\beta_2\) coincide con il quadrato del test t costruito sul singolo coefficiente nel modello completo

Il test (vedere dispense della teoria) sarebbe dato da: \(t_{\hat{\beta_2}}=\frac{\hat{\beta_2}-0}{s.e.(\hat{\beta_2})}\)

L’errore standard \(s.e.(\hat{\beta_2})\) è calcolato dall’opportuno elemento della matrice di varianza e covarianza \(V(\hat{\mathbf{\beta}})=\sigma^2 (X'X)^{-1}\)

Riprendere la parte teorica

source 
s2=sum((l1$residuals^2))/1420
Sx=cov(as.data.frame(scale(antropometric[,c(7,9:13)])))
Cx=solve(Sx)
kable(round(Sx,2))
ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
ALTEZZA 1.00 0.59 0.48 0.75 0.74 0.89
TORACECM 0.59 1.00 0.48 0.69 0.77 0.58
CRANIOCM 0.48 0.48 1.00 0.50 0.49 0.47
BISACROM 0.75 0.69 0.50 1.00 0.76 0.78
BITROCAN 0.74 0.77 0.49 0.76 1.00 0.71
SPANCM 0.89 0.58 0.47 0.78 0.71 1.00
source 
kable(round(Cx,2),caption="Inversa di Cor(X)")
Inversa di Cor(X)
ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
ALTEZZA 5.57 0.18 -0.17 -0.26 -1.17 -3.96
TORACECM 0.18 2.71 -0.31 -0.74 -1.59 0.12
CRANIOCM -0.17 -0.31 1.44 -0.22 -0.08 -0.12
BISACROM -0.26 -0.74 -0.22 3.58 -0.87 -1.40
BITROCAN -1.17 -1.59 -0.08 -0.87 3.81 -0.02
SPANCM -3.96 0.12 -0.12 -1.40 -0.02 5.62
source 
print(Cx[2,2]*s2/n)
output 
[1] 0.0001967018
source 
print(c(s2,n))
output 
[1]    0.1036499 1427.0000000
source 
print(l1$coefficients[3]/sqrt(Cx[2,2]*s2/n))
output 
TORACECM 
45.09379
source 
kable(round(Cx,2),caption="Inversa di Cor(X)")
Inversa di Cor(X)
ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
ALTEZZA 5.57 0.18 -0.17 -0.26 -1.17 -3.96
TORACECM 0.18 2.71 -0.31 -0.74 -1.59 0.12
CRANIOCM -0.17 -0.31 1.44 -0.22 -0.08 -0.12
BISACROM -0.26 -0.74 -0.22 3.58 -0.87 -1.40
BITROCAN -1.17 -1.59 -0.08 -0.87 3.81 -0.02
SPANCM -3.96 0.12 -0.12 -1.40 -0.02 5.62
source 
print(diag(Cx))
output 
 ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN   SPANCM 
5.574349 2.708093 1.436256 3.582109 3.806666 5.619999
source 
print(diag(1/Cx))
output 
  ALTEZZA  TORACECM  CRANIOCM  BISACROM  BITROCAN    SPANCM 
0.1793931 0.3692635 0.6962548 0.2791652 0.2626970 0.1779360
source 
print(sum(diag(Cx)))
output 
[1] 22.72747
source 
lambda=eigen(Sx)$values
print(lambda)
output 
[1] 4.2688245 0.6577968 0.5383502 0.2348599 0.1982743 0.1018944
source 
print(1/lambda)
output 
[1] 0.2342565 1.5202264 1.8575269 4.2578578 5.0435184 9.8140862
source 
print(sum(1/lambda))
output 
[1] 22.72747

1.10 Multicollinearità

torna alla panoramica Ancora sulla struttura della matrice \(\xx\)