1 Test per la divergenza dalla linearità

1.1 Inizializzazione:

vengono caricati i packages necessari per realizzare questo documento

Warning in smooth.spline(y ~ x, cv = TRUE): cross-validation with non-unique 'x' values seems doubtful

1.2 Selezione dei dati

Utilizziamo il dataset children1 di 26039 osservazioni.

Rappresentiamo la relazione fra peso alla nascita e numero di settimane di gestazione senza effettuare alcuna selezione. Evidentemente l’insieme utilizzato contiene molti errori di rilevazione.

Inoltre per questo esempio, da prendere esclusivamente come spunto per introdurre l’argomento della divergenza dalla linearità, 1 è opportuno selezionare solo le osservazioni complete, in un range di valori fra 33 e 42 settimane di gestazione, e con ulteriori criteri di selezione nel tentativo di lavorare con un insieme dei dati approssimativamente omogeneo2.

ch1=complete(children1)

eserc1=ch1[(ch1$Complicanzioni==0)&(ch1$Ospedale==1)&(ch1$modo.di.arrivo==0)&(ch1$X.decesso==0)&(ch1$X.gravidanze==0)&(ch1$X...Parto==0)&(ch1$Fumatrici==0)&(ch1$gestazione>32)&(ch1$gestazione<43),]
attach(eserc1)
table(gestazione)
gestazione
  33   34   35   36   37   38   39   40   41   42 
   6   31   41   61  230  596 1146 1599  692  102 
dim(eserc1)
[1] 4504   26
with(eserc1,boxplot(peso~gestazione))

1.3 Analisi della varianza con un criterio di classificazione quantitativo

Nella figura è rappresentata quindi la distribuzione dei pesi del campione di 4504 neonati, suddivisi secondo la durata della gestazione, misurata in settimane. Appare evidente che il peso medio dipende dal numero di settimane della gestazione, ed un test di analisi della varianza ad una via risulterà certamente significativo. Tuttavia in questo caso, dato che il criterio di classificazione (ossia il numero di settimane della gestazione) è ovviamente quantitativo, ovviamente possiamo fare di più:

Ci possiamo chiedere se una retta descrive adeguatamente la dipendenza in media dei pesi dal numero di settimane di gestazione, o se invece non vi sia una significativa divergenza dalla linearità. In effetti in questo caso la rappresentazione sembra suggerire che una relazione polinomiale descriva meglio l’andamento delle medie, e più avanti affronteremo anche questa possibilità.

Ovviamente possiamo condurre questa analisi in termini di test di significatività perchè siamo in presenza di per ogni modalità della variabile esplicativa: se avessimo osservazioni singole per ciascuna modalità potremmo giudicare della linearità di una relazione con altri strumenti (grafici o mediante analisi dei residui a posteriori o mediante confronto con regressioni non parametriche).

1.4 Boxplot delle due variabili

Lavoriamo quindi con il dataset ridotto eserc1 di 4504 osservazioni.

Miglioriamo il boxplot usando la function MLA.boxplot del package didattico. Oltre il classico boxplot, sono rappresentati i singoli punti, le medie, la spezzata di regressione e la retta di regressione.

output1=MLA.boxplot(peso,gestazione)

1.5 Due modelli lineari… per le stesse variabili

Ora sfruttiamo il fatto che la variabile gestazione può essere vista come un criterio di classificazione quantitativo (oppure come una variabile quantitativa con modalità ripetute…), e adattiamo due modelli lineari in cui la variabile di risposta è sempre la stessa (peso) e la variabile esplicativa è sempre la variabile gestazione, prima vista come fattore di classificazione (modello m0) e poi come variabile quantitativa (modello m1)

m0=lm(peso~as.factor(gestazione))
print(m0)

Call:
lm(formula = peso ~ as.factor(gestazione))

Coefficients:
            (Intercept)  as.factor(gestazione)34  as.factor(gestazione)35  as.factor(gestazione)36  
                 2266.7                     17.2                    248.7                    482.5  
as.factor(gestazione)37  as.factor(gestazione)38  as.factor(gestazione)39  as.factor(gestazione)40  
                  691.2                    847.5                    989.7                   1104.2  
as.factor(gestazione)41  as.factor(gestazione)42  
                 1202.4                   1212.1  
m1=lm(peso~gestazione)
print(m1)

Call:
lm(formula = peso ~ gestazione)

Coefficients:
(Intercept)   gestazione  
    -2042.8        135.2

Intanto abbiamo una prima informazione sintetica dalla stima del modello di regressione lineare m1, ossia che il peso aumenta media mente di -2042.8477545 grammi la settimana, nel range di settimane considerate, ovviamente, poichè non potremmo mai usare questo valore al di fuori del range osservato della x.

n=length(peso)
k=length(table(gestazione))
dev.tot=sum((peso-mean(peso))^2)
dev.spieg.m0=sum((m0$fitted.values-mean(peso))^2)
dev.spieg.m1=sum((m1$fitted.values-mean(peso))^2)

1.6 Un test di divergenza dalla linearità

Vediamo quindi che il modello m0 spiega una frazione \(\eta^2_{yx} = 0.196\) della devianza totale della variabile peso, mentre il modello m1 spiega una frazione \(r^2_{yx} =0.186\).

Abbiamo dunque:

devianza totale 821220184.87
devianza spieg. m0 160923737.23
devianza spieg. m1 152693406.95
devianza res. m0 660296447.64

1.7 Il test con l’istruzione ‘anova’

Possiamo quindi semplicemente saggiare l’ipotesi di linearità fra le medie, contro l’ipotesi di divergenza dalla linearità, mediante l’istruzione:

anova(m1,m0)
Analysis of Variance Table

Model 1: peso ~ gestazione
Model 2: peso ~ as.factor(gestazione)
  Res.Df       RSS Df Sum of Sq     F    Pr(>F)    
1   4502 668526778                                 
2   4494 660296448  8   8230330 7.002 3.217e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.8 Formule esplicite

Che test abbiamo impiegato (e perchè)?

Si può costruire il test per la verifica dell’ipotesi di linearità della regressione: \[ F=\frac{\frac{\sum_{j=1}^{k}n_{j} (M_{j}-\hat{M}_{j})^{2}}{k-2}} {\frac{\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j})^{2}}{n-k}} =\frac{\frac{8.2303303\times 10^{6}}{8}}{\frac{6.6029645\times 10^{8}}{4494}}=7.0019883 \] E’ facile vedere che questo test si può esprimere anche in questo modo, evidenziando il ruolo delle frazioni di devianza \(e_{yz}^{2}-r_{yz}^{2}\)

Il valore calcolato del test F chiaramente è uguale a quello trovato con anova(m1,m0) \[ F=\frac{\frac{\eta_{yz}^{2}-r_{yz}^{2}}{k-2}}{\frac{1-\eta_{yz}^{2}}{n-k}} =\frac{\frac{0.1959569-0.1859348}{8}}{\frac{1-0.1959569}{4494}}= 7.0019883 \] Avendo indicato con \(y, z\) rispettivamente la variabile dipendente e quella esplicativa.

Il valore ottenuto è significativo, pertanto rifiutiamo l’ipotesi di linearità.

1.9 Divergenza da una curva di grado superiore

In effetti se ritorniamo al grafico delle medie , notiamo un andamento non lineare con un possibile punto di flesso.

Proviamo allora ad adattare due nuovi modelli:

  • m2 con una componente cubica (il minimo per potere avere un punto di flesso, come sembrano mostrare i dati).

  • ed un ulteriore modello m3 in cui adattiamo un polinomio di terzo grado completo.

m2=lm(peso~gestazione+I(gestazione^3))
m3=lm(peso~gestazione+I(gestazione^2)+I(gestazione^3))
print(m2)

Call:
lm(formula = peso ~ gestazione + I(gestazione^3))

Coefficients:
    (Intercept)       gestazione  I(gestazione^3)  
     -1.499e+04        6.382e+02       -1.119e-01  
print(m3)

Call:
lm(formula = peso ~ gestazione + I(gestazione^2) + I(gestazione^3))

Coefficients:
    (Intercept)       gestazione  I(gestazione^2)  I(gestazione^3)  
     20374.0419       -2149.5473          73.0917          -0.7494

3

Costruiamo le due polinomiali e le rappresentiamo graficamente sopra in boxplot precedente (in blue la polinomiale completa, in verde quella incompleta, senza il termine di secondo grado)

Vediamo che le due polinomiali sono molto simili ed effettuiamo il test per vedere se m2 o m3 sono sufficient1 per spiegare la variabilità delle medie, ancora impiegando la funzione anova (che quando si mettono a confronto due modelli nidificati fornisce il test F relativo ad una \(H_0\) qualsiasi):

anova(m2,m0) # divergenza da una relazione cubica
Analysis of Variance Table

Model 1: peso ~ gestazione + I(gestazione^3)
Model 2: peso ~ as.factor(gestazione)
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1   4501 660901739                           
2   4494 660296448  7    605291 0.5885 0.7659
anova(m3,m0) # divergenza da una relazione cubica completa
Analysis of Variance Table

Model 1: peso ~ gestazione + I(gestazione^2) + I(gestazione^3)
Model 2: peso ~ as.factor(gestazione)
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1   4500 660787809                           
2   4494 660296448  6    491361 0.5574 0.7645

1.10 Modello finale (cubico)

Entrambi i test risultano non significativi, per cui già il modello m2 (con tre soli parametri) sarà sufficiente per spiegare la variabilità delle medie!

\[ \hat{E}(Y_j)=\beta_0+\beta_1 x_j +\beta_3 x_j^3= -1.499484\times 10^{4}+638.1699588 x_j +-0.1118815 x_j^3 \]

1.11 Scomposizione della devianza per un criterio di classificazione quantitativo

Riassumo qui un po’ della teoria impiegata per affrontare l’esempio precedente (almeno per la parte di divergenza della linearità: il principio poi si può applicare al confronto con una polinomiale o altra curva parametrica qualsiasi)

1.11.1 Scomposizione della devianza empirica in 3 componenti

Se un criterio di classificazione è numerico, con \(k\) livelli \(z_{j}\), si può scomporre come segue la devianza fra i gruppi:

\[ \sum_{j=1}^{k}n_{j} (M_{j}-M)^{2}, \] per vedere quanta parte di essa è spiegata dalla regressione lineare delle medie di \(\mathbf{y}\) \(M_{j}\) sui valori \(z_{j}.\)

Se \(\hat{M}_{j}\) è il valore stimato dalla regressione lineare:

\(\hat{M}_{j}=b_{yz}z_{j}\),

si può dimostrare che si può operare algebricamente la scomposizione:

Fra i gruppi Divergenza dalla linearità Regressione lineare di su \(z\)
\(\sum_{j=1}^{k}n_{j} (M_{j}-M)^{2}=\) \(\sum_{j=1}^{k}n_{j} (M_{j}-\hat{M}_{j})^{2}+\) \(\sum_{j=1}^{k}n_{j} (\hat{M}_{j}-M)^{2}\)

la somma dei doppi prodotti è nulla per le equazioni normali.

I due termini a destra si distribuiscono (sotto \(H_{0}\)) secondo due v.c. \(\chi^{2}\) indipendenti, con \(k-2\) e 1 grado di libertà. 4

In effetti possiamo considerare lo scarto di ogni osservazione dalla media aritmetica generale \((y_{ij}-M)\); aggiungendo e sottraendo \(M_{j}\) e \(\hat{M}_{j}\) e riordinando i termini otteniamo tre componenti ciascuna con un suo preciso significato:

\[ \underbrace{(y_{ij}-M)}_{\mbox{scarto totale}} = \underbrace{(y_{ij}-M_{j})}_{\mbox{ scarto dalla media del gruppo }}+ \underbrace{(M_{j}-\hat{M}_{j})}_{\mbox{ Divergenza dalla linearità }} + \underbrace{(\hat{M}_{j}-M)}_ {\mbox{ Regressione lineare }} \]
In definitiva nel caso di un fattore quantitativo \(Z\) possiamo scomporre la devianza totale in tre parti (avendo indicato con \(e_{yz}\) il rapporto di correlazione che misura la dipendenza in media di \(y\) da \(z\)):

Tavola di scomposizione della devianza empirica di \(\mathbf{y}\) per un criterio di classificazione semplice con \(k\) livelli quantitativi \(z_{j}\)

Dev TIPO g.d.l. Prop. di dev. tot. Elem.
\(\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M)^{2}\) Totale \(n-1\) 1 \(\mathbf{y}_{ij}-M\)
\(\sum_{j=1}^{k}\sum_{i=1}^{n_j}(y_{ij}-M_{j})^{2}\) Entro i gruppi \(n-k\) \(1-e_{yz}^{2}\) \(\mathbf{y}_{ij}-M_{j}\)
\(\sum_{j=1}^{k}n_{j} (M_{j}-\hat{M}_{j})^{2}\) Divergenza dalla linearita’ \(k-2\) \(e_{yz}^{2}-r_{yz}^{2}\) \(M_{j}-\hat{M}_{j}\)
\(\sum_{j=1}^{k}n_{j} (\hat{M}_{j}-M)^{2}\) Regressione lineare di \(\mathbf{y}\) su \(z\) 1 \(r_{yz}^{2}\) \(\hat{M}_{j}-M\)

1.11.2 Differenza fra i test di omogeneità

Il test con \(r_{yz}^{2}\) (ossia quello sulla regressione lineare semplice) saggia l’ipotesi che le \(k\) medie non varino, contro l’alternativa che varino in modo lineare rispetto a \(z\). Saggia quindi l’ipotesi che \(\beta_{yz}=0\), contro l’alternativa che \(\beta_{yz} \neq 0\), essendo comunque i valori attesi di \(\mathbf{y}\) funzioni lineari di \(Z\).

Il test con \(e_{yz}^{2}\) (ossia quello sull’analisi della varianza a una via) saggia l’ipotesi che le \(k\) medie non varino, contro l’alternativa che varino in modo qualsiasi (anche non linearmente rispetto a \(z\)): saggia l’ipotesi che \(\mu_{j} = \mu_{r}, \ \forall i \neq r\), contro l’alternativa che per almeno due gruppi si abbia: \(\mu_{j} \neq \mu_{r}\), essendo i \(k\) valori attesi di \(\mathbf{y}\) funzioni qualsiasi di \(Z\))

In linea generale questi due test dovrebbero differire per quanto riguarda il potere, dal momento che si riferiscono ad alternative differenti.

Il test su \(e_{yz}^{2}-r_{yz}^{2}\) saggia l’ipotesi che le \(k\) medie varino solo per effetto di una relazione lineare rispetto a \(z\), contro l’alternativa che varino in modo non lineare.

Per esempio, supponendo l’esistenza di una relazione polinomiale di grado \(k-1\) dei valori attesi di \(\mathbf{y}\) rispetto a \(Z\).

\[ \mathrm{E}\left[Y_{i}\right]= \sum_{j=0}^{k-1}\beta_{j}z_{i}^{j}: \]

saggia l’ipotesi che i \(k-2\) coefficienti dei termini di grado \(2^o\) e superiore siano nulli,

\[ H_{0}: \beta_{2}= \beta_{3}= \dots = \beta_{k-1}=0; \qquad \forall \ \beta_{0},\beta_{1} \]

ossia l’ipotesi che la relazione sia lineare, contro l’alternativa che almeno un coefficiente sia diverso da zero, ossia che la relazione sia curvilinea.

  1. e senza pretese di trarre conclusioni generali sulle leggi della crescita intrauterina↩︎

  2. Questi dati, come altri nel corso, non costituisono a rigore un campione casuale.↩︎

  3. la notazione I(x^r) si usa per trasformazioni semplici di variabili in lm()↩︎

  4. Si può applicare il teorema di Cochran (perché la somma dei gradi di libertà coincide col totale)↩︎