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