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

Generated by summarytools 0.9.6 (R version 4.0.2)
2020-09-16

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)

You must enable Javascript to view this page properly.

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)

You must enable Javascript to view this page properly.

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)
k=ncol(X)
X=as.matrix(cbind(array(1,n),X))
round(t(X)%*%X/(n),2)

output

            array(1, n) ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
array(1, n)           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

b=solve(t(X)%*%X)%*%t(X)%*%Y
print(b)

output

                  [,1]
array(1, n) 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

            array(1, n) ALTEZZA TORACECM CRANIOCM BISACROM BITROCAN SPANCM
array(1, n)           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]
array(1, n) 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.91036

source

print(sqrt(s2))

output

[1] 3.45114

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) -5.045e-16  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] -5.04e-16 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] -1.40e-14 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) 1.393e-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) -5.045e-16  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 \(\mathbf{X}^{\mathsf{T}}\mathbf{X}\)