Parte IIIvenus.unive.it/stone/statcomp/lucidi/seconda.pdf · 2008-02-15 · Figura 3: I dati...

Regressione lineare semplice

Parte III

Regressione lineare


Scopo

Siamo interessati alle relazioni che esistono tra due variabili:

Y , variabile quantitativa, detta variabile risposta;

X , variabile esplicativa (o variabile indipendente o predittore oregressore o, piu raramente, covariata). Per il momentoassumeremo che anche questa sia una variabile quantitativa,ma questa assunzione puo essere rimossa.

In particolare, vogliamo capire se X eserciti qualche influenzasu Y e come tale influenza possa manifestarsi.


Il modello statistico

Abbiamo bisogno di definire un modello statistico che leghi Y a X

Modello statistico

Rappresentazione semplificata, ma auspicabilmente soddisfacente,della realta osservabile



Ipotizziamo che Y possa essere modellata come una variabilecasuale: le osservazioni yi che rileviamo sulle unita statistichesono in qualche modo determinate dal caso.

La variabile X e invece considerata deterministica, noncasuale.

Sull’unita statistica i-esima la modalita xi e misurata senzaerrore



Come si manifesta l’influenza di X su Y ?

E(Y |X = x) = η(x ,β), β ∈ Rm, m ≥ 1

η(x ,β) si dice funzione di regressione: essa rappresenta unmodello di comportamento della media di Y come funzione diX

β e un vettore di parametri da cui dipende η(x ,β)



Nel modello di regressione lineare semplice abbiamo che

η(x ,β) = β0 + β1x , β = [β0, β1]′, β0, β1 ∈ R

Si assume inoltre che

Var(Y |X = x) = σ2 > 0 costante rispetto a x

L’ultima assunzione importante riguarda l’incorrelazione trale variabili casuali che generano le osservazionicampionarie di Y . Supponendo di disporre di un campione didimensione n,

Corr(Yi ,Yj) = 0 ∀i 6= j , i , j = 1, . . . , n



Le n variabili casuali campionarie Yi avranno medie diverse,dipendenti dal valore assunto da xi , e varianze uguali.

Saranno inoltre incorrelate.

Si veda la Figura 2: in ascissa si riportano i valori di X , in ordinataquelli di Y ; i cerchi hanno come coordinate le osservazioni xi e yi .Sul medesimo grafico si riportano la retta di regressione e ledistribuzioni condizionate Yi |X = xi


Figura 2: Distribuzioni di probabilita Yi condizionate a X = xi

−3 −2 −1 0 1 2 3

−10

−5

05

1015

20

x

y

●

●

●

●

●

●

retta di regressioneosservazioni campionarie

ηη((x,, ββ))


Il grafico in figura 2 ci dice una cosa molto importante.Compatibilmente con il modello che abbiamo costruito, lamodalita di Y che osserviamo in corrispondenza di una fissatamodalita di X e casuale.

Se avessimo estratto un campione diverso, avremmo osservatoqualcosa di diverso.

Tutte le conclusioni che possiamo trarre da cio che osserviamosono quindi determinate dal caso.

Dobbiamo essere consapevoli di agire in condizioni diincertezza e le nostre decisioni devono tener conto di questo.



Le assunzioni che abbiamo introdotto nei lucidi precedenti ciportano a definire sinteticamente un modello di questo tipo:

Yi = η(xi ,β) + εi

= β0 + β1xi + εi , i = 1, . . . , n

εi sono variabili casuali tali che:

E(εi ) = 0, Var(εi ) = σ2, ∀i

eCorr(εi , εj) = 0 ∀i 6= j .

Le determinazioni di εi sono dette disturbi aleatori


Alcuni dati simulati

Il file simulati.txt contiene un campione di dimensione 5generato dal modello:

Yi = 5 + 3xi + εi , εi ∼ N(0, 100)

Salvato il file in una directory sulla quale operiamo con R,importiamo i dati nel dataframe con il comando:

> simulati = read.table("simulati.txt",header=T)

Possiamo visualizzare i dati sullo schermo digitando:

> simulati

La figura 3 mostra i dati e la retta di regressione. Pone inoltrein evidenza il ruolo di β0, β1 e εi .


●

●

●

●

●

−1 0 1 2

−10

−5

05

1015

Yi == 5 ++ 3xi ++ εεi

x

y

ββ0

ββ1

ηη((x,, ββ)) == 5 ++ 3x

εε1

εε2εε5

εε4 εε3

Figura 3: I dati simulati e la retta di regressione


Stima dei parametri

Nella realta non conosceremo mai i parametri del modello: β0,β1 e σ2.

Essi devono essere stimati sulla base dei dati campionari.


Stima dei parametri

Una regola generale: fissato il campione, determiniamo lestime di β0 e di β1 minimizzando rispetto ai parametri lafunzione:

S(β0, β1) =n∑

i=1

(yi −β0−β1xi )2 metodo dei minimi quadrati

Dai corsi di base sappiamo che tali stime sono:

b0 = y − b1x e b1 =

∑ni=1(xi − x)yi∑ni=1(xi − x)2

Per i dati simulati, otteniamo le stime: b0 = 3.01 e b1 = 3.03

In Figura 4 sono rappresentati i dati, la funzione di regressionevera e quella stimata. Che deduzioni possiamo trarre?


●

●

●

●

●

−1 0 1 2

−10

−5

05

1015

x

y

● punti osservatiretta di regressioneretta di regressione stimata

Figura 4: Dati simulati: funzione di regressione vera e funzione diregressione stimata


Con il metodo dei minimi quadrati cerchiamo il valori di β1 e β2

che rendono minima

S(β0, β1) =n∑

i=1

(yi − β0 − β1xi )2

e abbiamo visto che quei valori sono b0 e b1, ovvero abbiamodefinito il vettore di stime b = [b0, b1]

′

Definiamo ora:

yi = η(xi ,b) = b0 + b1xi stima della funzione di regressionesull’unita statistica i-esima, corrispondente alla modalita xi .

ei = yi − yi residuo di regressione sull’unita statistica i-esima.


I residui della regressione

I residui rappresentano quella parte di variabilita campionariadi Y che non e spiegata dal modello di regressione (torneremopiu avanti su questo punto).

Sono molto importanti, come vedremo in seguito.

Non vanno confusi con i disturbi aleatori, εi : i residui sonosempre calcolabili dopo aver ottenuto le stime dei coefficientidel modello; i disturbi non sono mai determinabili perche percalcolarli dovremmo conoscere β0 e β1.


Un primo utilizzo dei residui: stima di σ2

La stima di σ2 e data da

s2 =1

n − 2

n∑i=1

e2i

Nel modello che abbiamo stimato con i dati simulati si ottienes2 = 59.34


Il grafico in figura 5 riassume (quasi) tutto quello che abbiamodetto sulla regressione lineare semplice.


●

●

●

●

●

−1 0 1 2

−10

−5

05

1015

x

y b0

b1

ηη((x,, b)) == 3.01 ++ 3.03x

e1

e2e5

e4 e3

Figura 5: Una sintesi grafica della stima di un modello di regressionelineare semplice


Il modello in notazione matriciale

Premessa

Il modello di regressione lineare semplice che abbiamo introdotto inprecedenza come

Yi = β0 + β1xi + εi , i = 1, . . . , n

puo essere riscritto come segue:

Yi = [1, xi ]β + εi , i = 1, . . . , n

= x′iβ + εi , i = 1, . . . , n

conx′i = [1, xi ] e β = [β0, β1]

′



vettore delle n variabili casuali campionarie: Y = [Y1, . . . ,Yn]′

vettore delle n osservazioni campionarie sulla variabile Y :y = [y1, . . . , yn]

′

matrice di regressione:

X =

1 x1...

...1 xi...

...1 xn

=

x′1...x′i...x′n



Vettore dei disturbi aleatori: ε = [ε1, . . . , εi , . . . , εn]′ tale che

E(ε) = 0 e Σε = σ2In


Il modello in forma matriciale

In forma matriciale, quindi il modello puo essere scritto come:

Y = Xβ + ε,

ovveroE(Y|X) = Xβ e ΣY|X = σ2In


Le stime in forma matriciale

Per stimare β abbiamo minimizzato la funzione

S(β) =n∑

i=1

(yi − β0 − xiβ)2

la quale puo essere scritta in forma vettoriale come:

S(β) = (y − Xβ)′(y − Xβ)



Sotto condizioni assai generali (che vedremo fra poco), il vettoreche minimizza S(β) e

b = (X′X)−1X′y ∈ R2.

Gli elementi che costituiscono i vettore b coincidono esattamentecon le stime che abbiamo definito in precedenza, cioe:

b = [b0, b1]′.



Attenzione

b e calcolabile sotto due condizioni:

n ≥ 2: il numero di osservazioni campionarie deve essere noninferiore al numero dei coefficienti che vogliamo stimare.

Le colonne di X devono essere linearmente indipendenti.Questo significa che la variabile X non puo assumere lo stessovalore su tutte le unita statistiche osservate.



Il vettore dei residui e definito come

e = [e1, . . . , ei , . . . , en]′

E facile verificare che

e = y − Xb

= y − X(X′X)−1X′y

= (In − X(X′X)−1X′)y

= My

con M = (In − X(X′X)−1X′)



Come abbiamo gia visto, la stima della varianza e data da

s2 =1

n − 2

n∑i=1

e2i

ovvero

=1

n − 2e′e


L’incertezza, che sciagura!

Ricapitoliamo:

b = (X′X)−1X′y stim di β

e = My residui

s2 =1

n − 2e′e stima di σ2

Tutte le stime dipendono da y (oltre che da X)

y e una particolare determinazione del vettore casuale Y

Se avessimo osservato un campione diverso le stime cheavremmo ottenuto sarebbero diverse.

La statistica non serve a niente!


L’incertezza, che sciagura!

Abbiamo definito delle regole generali che ci permettono dicalcolare le stime b e s2 su qualsiasi campione osservato, sottocondizioni non particolarmente restrittive. Ma allora,

b e la determinazione campionaria di una particolare variabilecasuale bivariata:

B = (X′X)−1X′Y = [B0,B1]′

s2 e la determinazione campionaria di una particolare variabilecasuale:

S2 =1

n − 2E′E

Anche i residui sono determinazioni di n variabili casuali:

E = MY ∈ Rn


E facile verificare che

E(B) = β e ΣB = σ2(X′X)−1

E(S2) = σ2 e Var(S2) =2σ4

n − 2

Poniamo H = (X′X)−1. Avremo quindi che

Var(Bi ) = σ2hi ,i , i = 0, 1

Nel modello di regressione lineare semplice, e facile verificare che

Var(B0) =σ2

n

(1 +

x2

var (x)

)e Var(B1) =

σ2∑ni=1(xi − x)2


Possiamo stimare Var(B0) e Var(B1)?

Stime:

s2B0

=s2

n

(1 +

x2

var (x)

)e s2

B1=

s2∑ni=1(xi − x)2

Basta sostituire s2 a σ2 nelle espressioni delle rispettivevarianze.

Stimatori:

S2B0

=S2

n

(1 +

x2

var (x)

)e S2

B1=

S2∑ni=1(xi − x)2


Generalmente siamo interessati (vedremo meglio in seguito) allostandard error degli stimatori:

SEB0 =√

Var(B0) e SEB1 =√

Var(B1)

stime degli standard error:

sB0 =√

s2B0

e sB1 =√

s2B1

stimatori degli standard error:

SB0 =√

S2B0

e SB1 =√

S2B1


Correttezza

Gli stimatori dei parametri ignoti sono corretti: se potessimoripetere il campionamento all’infinito e a parita delle condizioni, lemedie delle stima ottenute coinciderebbero con i veri valori di β0,β1 e σ2


Osservazione

Dalla correttezza dello stimatore B di β segue che

E(Ei ) = E (Yi )− E(x′iB)

= x′iβ − x′E(B)

= x′iβ − xiβ = 0 ∀i = 1, . . . , n

e quindiE(E) = 0.


Efficienza

Tra tutti gli stimatori di β0 e β1, ottenibili come trasformazionilineari di Y e corretti, B0 e B1 sono i piu precisi, nel senso chehanno varianza minima.


Definizione di consistenza

Uno stimatore, Tn, del parametro θ ∈ R, si dice consistente se

limn→∞

P(|Tn − θ| > ε) = 0 ∀ε > 0

Quando la dimensione del campione e molto grande, gliscostamenti tra stime e vero valore del parametro tendono adivenire molto piccole.


B0, B1 e S2, sono stimatori consistenti per β0, β1 e σ2 sesussistono tutte le condizioni che caratterizzano il modello diregressione lineare semplice e se

limn→∞

1

nX′X = C matrice definita positiva


Ma come stimiamo i parametri con R?

Abbiamo gia importato il dataframe simulati dal filesimulati.txt

Per semplificarci la vita diamo il comando

> attach(simulati)

costruiamo ora un nuovo oggetto con il seguente comando:

> simulati.lm = lm(y~x)

Se non avessimo utilizzato il comando attach(simulati),avremmo dovuto scrivere

> simulati.lm = lm(y~x, data=simulati)

per far capire a R in quale dataframe dovesse cercare levariabili y e x


Ma come stimiamo i parametri con R?

La funzione lm (= linear model) costruisce un oggetto diclasse lm, nel nostro caso simulati.lm, che contiene tutte leinformazioni che ci servono sul modello stimato.

Una sintesi importante di tale informazione si puo ottenerecon il comando seguente:

> summary(simulati.lm)

che fornisce il seguente output



Call:lm(formula = y ~ x, data = simulati)

Residuals:1 2 3 4 5

-8.065 5.823 -6.270 5.587 2.925

Coefficients:Estimate Std. Error t value Pr(>|t|)

(Intercept) 3.014 3.445 0.875 0.446x 3.027 2.678 1.130 0.341

Residual standard error: 7.703 on 3 degrees of freedomMultiple R-Squared: 0.2986, Adjusted R-squared: 0.06482F-statistic: 1.277 on 1 and 3 DF, p-value: 0.3406


Per il momento non siamo in grado di capire tutto quello che cidice questo output, ma vediamo alcuni primi dettagli:

La prima riga rievoca il comando con cui si e chiamata lafunzione lm

Segue l’elenco dei residui (questo avviene solo quando si ha uncampione piccolo, nel nostro cason = 5)

La tabella che sta sotto la scritta Coefficients: e moltoimportante. Gli elementi sulla colonna Estimate sono lestime dei coefficienti di regressione (a sinistra si vedono leindicazioni corrispondenti all’intercetta e alla variabile X ).


Gli elementi sulla colonna Std. Error sono le stime deglistandard error degli stimatori dei coefficienti di regressione.

Nella terza riga dal basso si leggeResidual standard error: 7.703 on 3 degrees offreedomSi tratta della radice quadrata della stima della varianza deidisturbi (

√s2); il numero di gradi di liberta (degrees of

freedom) e pari a

n − numero di coefficienti di regressione stimati

nel caso del nostro esempio, 5− 2 = 3.

Il resto lo vedremo in seguito


Alcune funzioni applicabili ad un oggetto di classe lm

A titolo di esempio utilizziamo l’oggetto simulati.lm

> coef(simulati.lm) fornisce il vettore delle stime deicoefficienti di regressione (b).

> fitted(simulati.lm) fornisce il vettore delle stime dellafunzione di regressione sulle n unita statistiche (y)

> resid(simulati.lm) fornisce il vettore dei residui (e)


Alcuni grafici utili

Disegnamo il grafico di Y , facendo in modo, con l’argomentoylim, che le ordinate siano comprese tra il piu piccolo ed ilpiu grande tra gli elementi dei vettori y e y; con l’argomentomain definiamo un titolo per il grafico

> plot(y,ylim=range(c(y,fitted(simulati.lm))),main="Variabile risposta e valori stimati")

Aggiungiamo ora, in rosso, i punti corrispondenti ai valoristimati della funzione di regressione

> points(fitted(simulati.lm),col="red")

La funzione points aggiunge dei punti ad un graficopreesistente.


Avremmo potuto dare dei comandi un po’ diversi, per esempio:

> plot(y,ylim=range(c(y,fitted(simulati.lm))),main="Variabile risposta e valori stimati")

> lines(fitted(simulati.lm),col="red")

oppure

> plot(y,ylim=range(c(y,fitted(simulati.lm))),type="l",main="Variabile risposta e valori stimati")

> points(fitted(simulati.lm),col="red")


Un altro grafico utile e il seguente:

> plot(x,y,main="Dati osservati e retta di regressione stimata")

> lines(x,fitted(simulati.lm))

Visualizza quanto la retta di regressione stimata sia in grado diinterpolare i dati osservati


Varianza spiegata

Domanda

Perche utilizzare un modello di regressione quando potremmoapprossimare i valori yi con la loro media campionaria, y (Figura6)?


●

●

●

●

●

−1 0 1 2

−10

−50

510

x

y

y

ηη((x,, b))

Figura 6: Diagramma di dispersione di Y contro X , retta di regressionestimata e media campionaria di Y


Varianza spiegata

Un modo comune di valutare la variabilita campionaria di Yconsiste nel calcolare

var(Y ) =1

n

n∑i=1

(yi − y)2

Il numeratore della frazione che definisce var(Y ) si dice devianzatotale di Y :

DEVTOT =n∑

i=1

(yi − y)2


Varianza spiegata

Se nel modello di regressione compare l’intercetta, e facileverificare che

ˆy =1

n

n∑i=1

yi = y e e =1

n

n∑i=1

ei = 0

Ne segue che le varianze campionarie dei valori stimati dellaregressione e dei residui sono date da:

var(Y ) =1

n

n∑i=1

(yi − y)2 var(E ) =1

n

n∑i=1

e2i


Varianza spiegata

Possiamo quindi definire

DEVREG =n∑

i=1

(yi − y)2 devianza di regressione

e

DEVRES =n∑

i=1

e2i devianza residua


Varianza spiegata

Si puo dimostrare che

var(Y ) = var(Y ) + var(E ),

ovveroDEVTOT = DEVREG + DEVRES


Varianza spiegata

Domanda

Supponiamo di approssimare yi , i = 1 con yi , i = 1, . . . , n.Quanta parte della variabilita di Y cogliere attraverso il modello diregressione?

Risposta

R2 =DEVREG

DEVTOTCoefficiente di determinazione lineare

= 1− DEVRES

DEVTOT


Varianza spiegata

Una proprieta importante

0 ≤ R2 ≤ 1

Se R2 = 1, allora yi = yi , i = 1, . . . , n: le osservazioni su Y esu X giacciono su una retta

Se R2 = 0, allora yi = y , i = 1, . . . , n: la retta (funzione diX ) che meglio approssima le osservazioni su Y e costante(indipendente da X ) con ordinata uguale a y


Varianza spiegata

Quanto piu R2 e vicino a 1, tanto meglio i valori yi sonoapprossimati da yi .

Quanto piu R2 e vicino a 0, tanto peggio i valori yi sonoapprossimati da yi .

Si puo dimostrare che

corr (X ,Y ) = segno(b1)∣∣∣√R2

∣∣∣


Come determiniamo R2 con R?

Nel nostro esempio, troviamo il valore di R2 nell’output di


in corrispondenza di

Multiple R-Squared: 0.2986

Nel nostro esempio abbiamo quindi che R2 = 0.2986


Che delusione!

In termini di varianza spiegata, i risultati della nostra primaapplicazione sembrano deludenti, ma...

La varianza dei disturbi, εi , nell’esempio simulato era assaielevata.

Abbiamo stimato i parametri con soli 5 dati. Con uncampione piu numeroso avremmo potuto fare molto meglio.


Il modello gaussiano

Ricordiamo la formula generale del modello:

Yi = x′iβ + εi , i = 1, . . . , n

Il modello si dice gaussiano se εi ∼ N(0, σ2), ferme restandotutte le altre assunzioni.La Gaussianita del modello implica che

Yi |X = xi ∼ N(x′iβ, σ2)

Condizionatamente a X = xi , Yi si distribuisce come una v.c.normale con media x′iβ e varianza σ2


Il modello gaussiano

Conseguenze dell’ipotesi di gaussianita:

Distribuzione di B0, B1 e Yi :

B0 ∼ N

(β0,

σ2

n

(1 +

x2

var (x)

))B1 ∼ N

(β1,

σ2∑ni=1(xi − x)2

)Yi ∼ N

(β0 + β1xi , σ

2

(1

n+

(xi − x)2∑ni=1(xi − x)2

))


Tre variabili casuali importanti:

B0 − β0

SB0

∼ tn−2 (1)

B1 − β1

SB1

∼ tn−2 (2)

Yi − β0 − β1xi

S

√1

n+

(xi − x)2∑ni=1(xi − x)2

∼ tn−2 (3)

si distribuiscono come t di Student con n − 2 gradi di liberta


Una breve digressione

Una variabile casuale distribuita come una t di Student con n gradidi liberta ha una funzione di densita rappresentata in 7:

La funzione di densita e simmetrica intorno allo 0

Quando n →∞ essa tende alla funzione di densita di unanormale con media 0 e varianza 1.


Una breve digressione

x

f(t)

0

0

Figura 7: Funzione di densita di probabilita di una v.c. tn


Quantili di una v.c. T ∼ tn

Quantile di ordine α di una variabile casuale

Fissato un numero α tale che 0 < α < 1, diremo che qα e ilquantile di ordine α della variabile casuale Y se

P(Y ≤ qα) = α


Quantili di una v.c. T ∼ tn

Il quantile di ordine α di una t di student con n gradi di libertasara indicato nel seguito con il simbolo: tn,α

E facile verificare che tn,α = −tn,1−α

Supponiamo di voler trovare un intervallo,I1−α, simmetricointorno allo 0 tale che P(T ∈ I1−α) = 1− α. Tale intervallosara dato da

I1−α = [−tn,1−α/2, tn,1−α/2]

si veda la Figura 8


x

f(t)

−− tn,, 1−−αα 2 0 tn,, 1−−αα 2

0

αα 2 αα 2

1 −− αα

Figura 8: Quantili di una t di student con n gradi di liberta


Intervalli di confidenza

Supponiamo di essere interessati ad un parametro, θ ∈ R.

Sulla base dei dati campionari possiamo ottenere una stimapuntuale di θ, diciamo t.

Tenendo conto dello stato di incertezza in cui dobbiamooperare, vogliamo definire un intervallo di valori plausibili perθ.



Fissato un numero α, 0 < α < 1, supponiamo di potercostruire un intervallo avente per estremi due variabili casuali,

I1−α = [A1−α(Y),B1−α(Y)]

tale che P(A1−α(Y) ≤ θ ≤ B1−α(Y)) = 1− α.

Se potessimo estrarre infiniti campioni di fissata dimensione,n, su ciascuno di essi potremmo calcolare una particolaredeterminazione dell’intervallo aleatorio:

[a1−α(y), b1−α(y)]

una frazione 1− α di questi intervalli conterrebbe il verovalore del parametro.



Sul campione osservato y possiamo quindi costruire unparticolare intervallo, come quello indicato nel puntoprecedente, che e chiamato un intervallo di confidenza per θ allivello di fiducia 1− α.


Intervallo di confidenza per βi , i = 0, 1 al livello di fiducia1− α

Quantita pivotale (v.c. la cui distribuzione di probabilita enota e non dipende dai parametri ignoti)

T =Bi − βi

SBi

∼ tn−2



Determinazione degli intervalli aleatori Ii ,1−α

1− α = P

(−tn−2,1−α/2 ≤

Bi − βi

SBi

≤ tn−2,1−α2

)esplicitando rispetto a βi

= P(Bi − SBitn−2,1−α/2 ≤ βi ≤ Bi + SBi

tn−2,1−α/2)

Quindi gli intervalli aleatori Ii ,1−α sono dati da

Ii ,1−α = [Bi − SBitn−2,1−α/2,Bi + SBi

tn−2,1−α/2]



Determinazione degli intervalli di confidenza sul campioneosservato:

[bi − sBitn−2,1−α/2, bi + sBi

tn−2,1−α/2]

Nel caso dell’esempio sui dati simulati, fissando un livello difiducia 1− α = 0.95, abbiamo:

intervallo di confidenza per β0: [−7.950, 13.977]intervallo di confidenza per β1: [−5.496, 11.549]


Alcune considerazioni importanti

I risultati numerici che abbiamo ottenuto meritano qualchecommento:

Affermare che la probabilita che β0 sia compreso tra -7.95 e13.98 sia pari a 0.95 non ha alcun senso.E l’intervallo aleatorio (puramente teorico) I0.95 a contenere ilvero valore β0 con probabilita 0.95, non l’intervallo osservatosul campione.

L’intervallo di confidenza che abbiamo costruito e moltoampio. Questo significa che la nostra incertezza sul parametroe notevole.

Considerazioni del tutto analoghe valgono per l’intervallo diconfidenza per β1.


Altri dati simulati

Importamo ora, nel dataframe simulati1, i dati contenuti nelfile simulati1.txt

Questi dati sono simulati dallo stesso modello dei precedenti,ma ora la dimesione campionaria e n = 100.

Ripetiamo l’analisi precedente con i seguenti comandi:

> simulati1 = read.table("simulati1.txt",header=T)> attach(simulati1)> simulati1.lm = lm(y~x)> summary(simulati1.lm)

ottenendo un output simile a quello mostrato nella diapositivaseguente.


> summary(simulati1.lm)

Call:lm(formula = y ~ x)

Residuals:Min 1Q Median 3Q Max

-26.6184 -7.0684 -0.3749 6.0160 24.7992

Coefficients:Estimate Std. Error t value Pr(>|t|)

(Intercept) 6.6413 0.9923 6.693 1.37e-09x 2.7205 0.5603 4.856 4.54e-06---

Residual standard error: 9.922 on 98 degrees of freedomMultiple R-Squared: 0.1939, Adjusted R-squared: 0.1857F-statistic: 23.58 on 1 and 98 DF, p-value: 4.543e-06


Quali sono ora le stime dei parametri β0, β1 e σ2?

Quanto vale R2 e perche continua ad avere un valore basso?

Ad un livello di fiducia 1− α = 0.95, gli intervalli diconfidenza ora sono:

per β0: [4.672, 8.6105]per β1: [0.937, 4.505]

Perche?


Intervallo di confidenza per η(x0, β)

Sia x0 un fissato valore della variabile X .

Vogliamo costruire un intervallo di confidenza, di livello 1− αper

η(x0,β) = β0 + β1x0 = x′0β

con x′0 = [1, x0]

Se il modello e gaussiano, allora

η(x0,B) ∼ N(x′0β, σ2x′0(X′X)−1x0

e la quantita pivotale sara data da:

η(x0,B)− η(x0,β)√S2x′0(X

′X)−1x0

∼ tn−2


Intervallo di confidenza per η(x0, β)

Ponendo

Sη(x0,B) = S√

x′0(X′X)−1x0

Ad un livello di fiducia 1− α, quindi, l’intervallo di cofidenzaper η(x0,β) sara il seguente:[

η(x0,b)− tn−2,1−α/2Sη(x0,B), η(x0,b) + tn−2,1−α/2Sη(x0,B)

]


Come si fa con R?

Continuiamo con l’esempio sui dati simulati1: abbiamo giastimato il modello e salvato i risultati nell’oggettesimulati1.lm

Supponiamo di volere determinare un intervallo di confidenzaper la funzione di regressione calcolata in X = 0.5 e in X = 1

Costruiamo un nuovo dataframe in cui l’unica variabile checompare deve avere lo stesso nome del regressore che e statoutilizzato per costruire l’oggetto simulati1.lm:

> nuovidati=data.frame(x=c(0.5,1))


Ecco come deve comparire a video il nuovo oggetto creato:

> nuovidatix

1 0.52 1.0

Ora costruiamo un oggetto, int.regr, entro il qualepossiamo individuare i due intervalli di confidenza cheintendiamo calcolare:

int.regr=predict(simulati1.lm,nuovidati,interval="confidence",level=0.95)


L’oggetto int.regr deve apparire cosı:

> int.regrfit lwr upr

1 8.001518 5.961486 10.041552 9.361745 7.111376 11.61211

La prima riga corrisponde al primo elemento del dataframenuovidati; la seconda riga corrisponde al secondo elemento dinuovidati (X = 0.5 e X = 1).sulla prima colonna, fit, leggiamo i valori delle stime puntualidella funzione di regressione corrispondenti a X = 0.5 e aX = 1.Sulla seconda e la terza colonna leggiamo rispettivamente gliestremi inferiore e superiore degli intervalli di confidenza perη(0.5,β) e η(1,β)


La funzione predict()

Per acquisire informazioni sull’applicazione della funzionepredict() ad un oggetto di classe lm, digitare

> help(predict.lm)

il primo argomento della funzione deve essere un oggetto diclasse lm (nel nostro caso, simulati1.lm)il secondo argomento deve essere il nome del dataframe in cuidevono comparire le modalita del regressore (nel nostro caso,nuovidati)il terzo argomento, interval, deve indicare il tipo di intervalloche si vuole costruire (nel nostro caso interval=”confidence”significa che vogliamo costruire un intervallo di confidenza)il quarto argomento, level, rappresenta il livello di confidenza(nel nostro caso, level = 0.95)


Intervallo di previsione per Y in corrispondenza di X = x0

Vogliamo prevedere il valore assunto da Y in corrispondenzadi un fissato valore di X , diciamo x0, misurato su una unitastatistica non inclusa nel campione e sulla quale Y non e statorilevato.

Vogliamo costruire un intervallo di previsione, di livello 1− αper

Y0 = β0 + β1x0 + ε0 = x′0β + ε0

con x′0 = [1, x0]

Se il modello e gaussiano, allora

Y0|x0 ∼ N(x′0β, σ2

Se prevediamo Y0 con Y0 = x′0B, l’errore di previsione sara

Y0 − x′0B ∼ N(0, σ2(1 + x′0(X

′X)−1x0

)


Intervallo di previsione per Y in corrispondenza di X = x0

La quantita pivotale per l’intervallo di previsione al livello difiducia 1− α sara data da:

Y0 − x′0B√S2 (1 + x′0(X

′X)−1x0)∼ tn−2

Ponendo

SY0 = S√

1 + x′0(X′X)−1x0

Ad un livello di fiducia 1− α, quindi, l’intervallo di previsioneper Y0 sara il seguente:[

η(x0,b)− tn−2,1−α/2SY0 , η(x0,b) + tn−2,1−α/2SY0

]


Come si fa con R?

Continuiamo con l’esempio sui dati simulati1.

Supponiamo di volere determinare due intervalli di previsioneper Y in corrispondenza di X = 0.5 e in X = 1 ad un livello diconfidenza pari a 0.95

Utilizziamo il dataframe nuovidati gia introdotto inprecedenza


Ora costruiamo un oggetto, int.prev, entro il qualepossiamo individuare i due intervalli di previsione cheintendiamo calcolare:

int.prev=predict(simulati1.lm,nuovidati,interval="prediction",level=0.95)


L’oggetto int.prev deve apparire cosı:

> int.prevfit lwr upr

1 8.001518 -11.79375 27.796782 9.361745 -10.45630 29.17979

La prima riga corrisponde al primo elemento del dataframenuovidati; la seconda riga corrisponde al secondo elemento dinuovidati (X = 0.5 e X = 1).sulla prima colonna, fit, leggiamo i valori delle previsionipuntuali di Y in corrispondenza di X = 0.5 e a X = 1.Sulla seconda e la terza colonna leggiamo rispettivamente gliestremi inferiore e superiore degli intervalli di confidenza per Yin corrispondenza di X = 0.5 e a X = 1.


La funzione predict()

Per acquisire informazioni sull’applicazione della funzionepredict() ad un oggetto di classe lm, digitare

> help(predict.lm)

il primo argomento della funzione deve essere un oggetto diclasse lm (nel nostro caso simulati1.lm)il secondo argomento deve essere il nome del dataframe in cuidevono comparire le modalita del regressore (nel nostro casonuovidati)il terzo argomento, interval, deve indicare il tipo di intervalloche si vuole costruire (nel nostro caso interval=”prediction”significa che vogliamo costruire un intervallo di previsione)il quarto argomento, level rappresenta il livello di confidenza(nel nostro caso level = 0.95)


Nella Figura 9 sono rappresentati:

i dati osservati

la vera retta di regressione (rosso)

gli intervalli di confidenza per la media (rosso) e di previsione(blu) in corrispondenza di tutti i valori osservati di X al livellodi fiducia 0.95


●

●

●

●

●

●

●●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●●

●

●

●●

●

●

●

●

●

●●

●

●

●

●●

●

●

●

●

−2 0 2 4

−20

−10

010

2030

40

x

y

Figura 9: Intervalli di confidenza per η(x ,β) e di previsione per Y allivello di fiducia 0.95

Parte IIIvenus.unive.it/stone/statcomp/lucidi/seconda.pdf · 2008-02-15 · Figura 3: I dati...

Documents

Transcript of Parte IIIvenus.unive.it/stone/statcomp/lucidi/seconda.pdf · 2008-02-15 · Figura 3: I dati...