Trasformazione Box-Cox: un’analisi basata sulla ...tesi.cab.unipd.it/42884/1/Rigon_Tommaso.pdf ·...

Universit degli studi di PadovaDipartimento di Scienze Statistiche

Corso di Laurea Triennale in

Statistica, Economia e Finanza

Relazione Finale

Trasformazione Box-Cox: unanalisi basata sullaverosimiglianza

Relatore Prof. Nicola SartoriDipartimento di Scienze Statistiche

Laureando Tommaso RigonMatricola N 1010510

Anno Accademico 2012/2013

Indice

Introduzione 9

1 Linferenza di verosimiglianza 111.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.2 Specificazione del modello . . . . . . . . . . . . . . . . . . . . 121.3 Assuzioni e notazioni . . . . . . . . . . . . . . . . . . . . . . . 121.4 Quantit di verosimiglianza . . . . . . . . . . . . . . . . . . . 15

1.4.1 Stima di massima verosimiglianza . . . . . . . . . . . . 151.4.2 Funzione punteggio . . . . . . . . . . . . . . . . . . . . 151.4.3 Informazione osservata e informazione attesa . . . . . . 16

1.5 Alcuni risultati asintotici . . . . . . . . . . . . . . . . . . . . . 161.5.1 Test collegati alla verosimiglianza . . . . . . . . . . . . 17

1.6 Verosimiglianza profilo . . . . . . . . . . . . . . . . . . . . . . 171.7 Metodo di Newton-Raphson . . . . . . . . . . . . . . . . . . . 18

1.7.1 Algoritmi numerici con R . . . . . . . . . . . . . . . . 19

2 Modello lineare e trasformazione di variabili 232.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.2 Il modello lineare . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2.1 La stima dei parametri . . . . . . . . . . . . . . . . . . 252.2.2 Inferenza sui parametri . . . . . . . . . . . . . . . . . . 26

2.3 Trasformazioni della variabile dipendente . . . . . . . . . . . . 262.4 Modalit dutilizzo . . . . . . . . . . . . . . . . . . . . . . . . 282.5 I dati cars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3 Box-Cox e lapproccio di verosimiglianza 333.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.2 Funzione di verosimiglianza . . . . . . . . . . . . . . . . . . . 343.3 Funzione punteggio . . . . . . . . . . . . . . . . . . . . . . . . 353.4 Matrice di informazione osservata . . . . . . . . . . . . . . . . 363.5 Stime di massima verosimiglianza . . . . . . . . . . . . . . . . 38

3.5.1 La varianza per . . . . . . . . . . . . . . . . . . . . . 393.6 Trasformazione Bickel e Docksum . . . . . . . . . . . . . . . . 403.7 Matrice di informazione attesa . . . . . . . . . . . . . . . . . . 42

3.7.1 Calcolo degli elementi di I() . . . . . . . . . . . . . . 423.8 Alcune considerazioni . . . . . . . . . . . . . . . . . . . . . . . 443.9 Test log-rapporto di verosimiglianza . . . . . . . . . . . . . . . 45

4 Verifiche e simulazioni 474.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.2 La funzione Boxcox . . . . . . . . . . . . . . . . . . . . . . . . 474.3 Verifica della correttezza del codice . . . . . . . . . . . . . . . 48

4.3.1 Stima per . . . . . . . . . . . . . . . . . . . . . . . . 504.3.2 Matrice di informazione osservata . . . . . . . . . . . . 504.3.3 Matrice di informazione attesa . . . . . . . . . . . . . . 52

4.4 Normalit degli stimatori . . . . . . . . . . . . . . . . . . . . . 534.5 Test log-rapporto di verosimiglianza profilo per . . . . . . . 574.6 Simulazioni ed intervalli di confidenza . . . . . . . . . . . . . . 574.7 Test log-rapporto di verosimiglianza . . . . . . . . . . . . . . . 594.8 Rianalisi dei dati cars . . . . . . . . . . . . . . . . . . . . . . 60

4.8.1 Intervalli di confidenza per la media e di previsione . . 624.8.2 Ulteriori analisi . . . . . . . . . . . . . . . . . . . . . . 64

5 Conclusioni 65

Bibliografia 67

A Codice R utilizzato 69

Elenco dei codici

4.1 Esempio di output per la funzione Boxcox utilizzando il data-frame cars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.2 Comandi per la stima del modello, approcci differenti . . . . . 60A.1 La funzione Boxcox . . . . . . . . . . . . . . . . . . . . . . . . 69A.2 Log-verosimiglianza profilo cambiata di segno . . . . . . . . . 74A.3 Simulazione per la correttezza della stima di . . . . . . . . . 75A.4 Log-verosimiglianza cambiata di segno . . . . . . . . . . . . . 75A.5 Simulazione per la verifica della matrice j() . . . . . . . . . . 75A.6 Test log-rapporto di verosimiglianza profilo per . . . . . . . 76A.7 Esempio di simulazione per la verifica della normalit dello

stimatore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76A.8 Simulazione per la distorsione delle deviazioni standard . . . . 77A.9 Log-verosimiglianza profilo per . . . . . . . . . . . . . . . . 78A.10 Test log-rapporto di verosimiglianza e intervalli di confidenza

per . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78A.11 Simulazione intervalli di confidenza con Wp() . . . . . . . . . 79

Elenco delle tabelle

1.1 Le prime 10 osservazioni dei dati stress . . . . . . . . . . . . 20

2.1 Sommario di R per la stima del modelloYi = 0 + 1xi + i 32

4.1 Statistiche descrittive per |1 2|, . . . . . . . . . . . . . . . 504.2 Statistiche descrittive per il vettore quadJ, . . . . . . . . . . . 514.3 Matrice j() per i dati cars . . . . . . . . . . . . . . . . . . . 534.4 Risultati della simulazione, convergenza delle matrici I() e

j() . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.5 Errore Medio per 1 . . . . . . . . . . . . . . . . . . . . . . . 554.6 Statistiche descrittive per la simulazione . . . . . . . . . . . . 594.7 Le stime del modello tramite i due metodi . . . . . . . . . . . 614.8 Output modello in cui arrotondato allintero pi vicino . . 64

Introduzione

Uno degli strumenti pi utilizzati in statistica il modello lineare. La suautilit e diffusione giustifica gli sforzi atti a comprendere il suo funzionamen-to. In questa relazione, di carattere prevalentemente teorico, viene analizzatauna particolare tecnica utilizzata nei modelli lineari: la trasformazione Box-Cox. Questo approfondimento si reso necessario poich, nelluso quotidiano,essa viene utilizzata trascurando volutamente alcuni aspetti. Lo scopo mo-strare come questo approccio porti a distorsioni rilevanti. Pur essendo gi notialcuni risultati nella letteratura statistica, questi vengono ripresi e illustratidettagliatamente, insieme ad alcuni esempi. Contestualmente, si fornisconodei codici che consentono di ripercorrere i calcoli effettuati e condurre analisisu dataset differenti.

Gli argomenti riguardanti la verosimiglianza compaiono trasversalmente eproprio per questo nel primo capitolo si presenta una rassegna delle principalinozioni connesse ad essa. Non sono informazioni direttamente collegate allarelazione ma strettamente necessarie per la sua comprensione. Viene inoltredefinita la notazione delle principali quantit di verosimiglianza.

Nel secondo capitolo vengono brevemente presentati il modello linearee la trasformazione Box-Cox. Si tratta di un introduzione volta a esporrei motivi per cui, in certi casi, necessario utilizzare una trasformazione divariabile tra quelle note. Quella indicata da Box-Cox occupa una posizionedi rilievo e viene infatti utilizzata, nel modo consueto, in un dataset, a scopoesemplificativo.

Il terzo capitolo il nocciolo della relazione. Si affrontano i calcoli analiticiche permettono di raggiungere le quantit di verosimiglianza specifiche dellatrasformazione Box-Cox. Si discute poi la sua legittimit che porter al suo

6

parziale abbandono, seguendo le argomentazioni di Bickel e Docksum.Nel quarto capitolo infine si segnalano, da un punto di vista empirico, le

distorsioni che lapproccio usuale comporta, riprendendo il dataset usato nelsecondo capitolo. Per far ci viene prima controllato, tramite simulazioni, ilcodice che ha reso possibile il calcolo di queste distorsioni.

Nel quinto capitolo si conclude riassumendo i risultati a cui si giunti.

Capitolo 1

Linferenza di verosimiglianza

1.1 IntroduzioneFare inferenza significa, letteralmente, trarre delle conclusioni sulla base

di alcune verit. In ambito statistico il punto di partenza sono i dati, da cuisi cerca di estrarre qualche informazione di carattere generale. Le quantitdinteresse sono molteplici e differiscono a seconda del contesto applicativo.Ad esempio, spesso si interessati a studiare landamento medio di un fe-nomeno. Dato che il punto di partenza sono indicazioni parziali, o meglio,quantit campionarie, naturale che le conclusioni che se ne trarranno sa-ranno soggette ad errore. V quindi la necessit di una tecnica che permettadi ottenere stime il pi precise possibile e di cui si sia in grado di determinarelaffidabilit.

Linferenza basata sulla verosimiglianza uno strumento estremamenteversatile, introdotto da Ronald Fisher (1890-1962), che presume lesistenza diun modello statistico F . In questo capitolo introduttivo verranno presentatealcune nozioni basilari e definite le notazioni per le quantit utilizzate.

La presentazione basata sui testi Azzalini (2008, cap. 2-4) e Pace eSalvan (2001, cap. 1-6).

1.2 Specificazione del modello 8

1.2 Specificazione del modelloLa prima fondamentale assunzione considerare il vettore delle osserva-

zioni y = (y1, . . . , yn) come una realizzazione di una variabile casuale Y condistribuzione di probabilit p0(y). Il modello generatore dei dati ignoto maspesso si possono esprimere delle considerazioni sul fenomeno in esame. Datala natura stocastica di Y , la stima del modello generatore dei dati tantopi accurata quanto pi si riesce a restringere la classe F alla quale p0(y)appartiene. Il modello si intende correttamente specificato se p0(y) F .

Il modello F pu appartenere a una delle seguenti classi: modello nonparametrico, modello semi-parametrico, modello parametrico. Nel seguito cisi occuper unicamente di questultimo, la cui definizione formale :

F = {p(y, ) : Rp},

per qualche p N. Lo spazio chiamato spazio parametrico. Si deveipotizzare che esista una relazione biunivoca tra ciascun modello ed il valoreassunto dal parametro. Questa propriet chiamata identificabilit.

1.3 Assuzioni e notazioniSe F un modello parametrico per i dati y con funzione del modello

py(y; ) e = (1, . . . , p) Rp, la funzione di verosimiglianzaL : R+ definita come:

L(; y) = c(y)p(y; ),

dove c(y) una costante che non dipende da . La funzione L(; y) ha quindila stessa forma della funzione di probabilit del modello in cui per y fissatoe libero di variare. Pur assumendo valori che corrispondono a probabilit,L(y; ) non quindi una funzione di probabilit.

La verosimiglianza in realt una classe di funzioni equivalenti che dif-feriscono per una costante moltiplicativa c(y) e si pu parlare dunque diverosimiglianze equivalenti.

1.3 Assuzioni e notazioni 9

Il sostegno empirico a favore di 1 rispetto a 2 misurato dal rapporto

L(1; y)

L(2; y),

che chiamato rapporto di verosimiglianza.Linterpretazione di L(; y) la seguente: a partire dalle osservazioni y

possibile stabilire quale tra 1 e 2 pi verosimile. Segue quindi che 1 preferibile a 1 se L(1; y) L(2; y). Ad esempio nella Figura 1.1, = 1.9 preferibile a = 1.2, infatti: L(1.9; y) L(1.2; y).

Unipotesi frequente che i dati y = (y1 . . . , yn) siano realizzazioni i.i.d,ovvero osservazioni indipendenti ed identicamente distribuite. Perci se p(yi; ) la distribuzione di probabilit marginale per la singola osservazione, lafunzione di verosimiglianza diventa

L(; y) =ni=1

p(yi; ).

Funzione di log-verosimiglianza

La funzione di log-verosimiglianza viene introdotta perch permette mag-giore facilit di calcolo pur mantenendo tutta linformazione contenuta inL(; y). definita come

l(; y) = logL(; y),

dove log() il logaritmo naturale. Se ne ricava inoltre la quantit

l(; y) = c(y) + log p(y; ),

da cui si deduce che le log-verosimiglianze si equivalgono a meno di unacostante additiva, costituendo anchesse una classe di funzioni equivalenti.Nel caso il campione sia costituito da osservazioni i.i.d., per le propriet deilogaritmi ne consegue che

l(; y) =ni=1

log p(yi; ).

1.3 Assuzioni e notazioni 10

Figura 1.1: Funzione di verosimiglianza esponenziale di media 0.5 con n = 50

1.4 Quantit di verosimiglianza 11

1.4 Quantit di verosimiglianza

1.4.1 Stima di massima verosimiglianza

Seguendo le condizioni poste in Azzalini (2008, pag. 88), si di fronte adun problema regolare di stima se:

il modello identificabile;

lo spazio parametrico un intervallo aperto di Rp;

le funzioni di probabilit (o densit) hanno tutte lo stesso supporto;

per la funzione di densit f si pu scambiare due volte il segno diintegrale con quello di derivata rispetto a .

In realt queste sono condizioni piuttosto deboli che si verificano nella granparte dei casi pratici. Nel caso oggetto di questa tesi, le condizioni di regola-rit sono soddisfatte.

Il valore detto stima di massima verosimiglianza se

L() L(), ,

ovvero se il punto di massimo assoluto per la funzione di verosimiglianza.Poich log() una trasformazione monotona crescente, ne consegue che punto di massimo per L() se e solo se lo per l().

1.4.2 Funzione punteggio

Si definisce funzione score o punteggio

l()

= l(). (1.1)

Sotto le condizioni di regolarit il punto di massima verosimiglianza va cerca-to tra le soluzioni dellequazione l() = 0, che prende il nome di equazionedi verosimiglianza. Inoltre, vale che

E[l()] = 0 e che Var[l()] = E[l()l()T ] = I().

1.5 Alcuni risultati asintotici 12

1.4.3 Informazione osservata e informazione attesa

La matrice di informazione osservata j() e quella di informazione attesaI() di dimensioni p p sono definite nel modo seguente

j() = l() =l()

T, I() = E[j()]. (1.2)

Per indicare la matrice inversa I()1 nel blocco (, ) si soliti indicarlacome I() .

1.5 Alcuni risultati asintoticiFino ad ora non si detto nulla circa i vantaggi di utilizzare un approccio

di verosimiglianza rispetto ad altri metodi. auspicabile che lo stimatorefornisca una stima sempre pi precisa dei parametri e che converga al verovalore al divergere della numerosit campionaria. Sotto le condizioni di re-golarit possibile dimostrare alcune propriet dello stimatore di massimaverosimiglianza. In particolare, (Y ) :

asintoticamente non distorto, ovvero: limnE[(Y )] = . In campio-ni finiti , in generale, distorto;

consistente, ovvero converge in probabilit a , (Y ) p e perci:limnVar((Y )) = 0;

asintoticamente efficiente.

Qualora non fosse possibile risalire alla distribuzione esatta dello stimatore, possibile ricorrere ad alcuni risultati asintotici. Per n sufficientemente grandevale che

Np(, I()1). (1.3)

La matrice I() in generale non nota e quindi la si pu sostituire sia conI() che con j() che ne costituiscono delle stime consistenti.

Questo risultato di fondamentale importanza perch ci si riconduce,almeno approssimativamente, ad una distribuzione nota. Ci permette di co-struire test dipotesi ed intervalli di confidenza anche quando la distribuzioneesatta troppo complessa da ricavare analiticamente.

1.6 Verosimiglianza profilo 13

1.5.1 Test collegati alla verosimiglianza

In un sistema dipotesi in cui si ha che H0 : = 0 e H1 : 6= 0 esistonotre test statistici asintoticamente equivalenti che, sotto H0, al divergere din, convergono in distribuzione ad una distribuzione 2p. Il primo chiamatolog-rapporto di verosimiglianza ed pari a

W (0) = 2(l() l(0)) = 2(l(0) l()). (1.4)

Il secondo chiamato test di Wald ed pari a

We(0) = ( 0)T I()( 0), (1.5)

inoltre il terzo chiamato test score o test di Rao ed pari a

Wu(0) = l(0)I(0)1l(0). (1.6)

Nel caso in cui p = 1 si ricorre alle versioni unilaterali di questi stessi te-st r(0), re(0), ru(0), che, sotto H0, si distribuiscono approssimativamentecome delle N(0, 1). Le definizioni dei test appena citati sono rispettivamente

r(0) = sgn( 0)W (0),

re(0) = ( 0)I(),

ru(0) = sgn( 0)l(0)I(0)12 .

1.6 Verosimiglianza profiloIn alcuni casi pratici, in cui p > 1, lattenzione rivolta verso un primo

gruppo di parametri detti dinteresse, indicati con . Tuttavia si deve co-munque tener conto dellesistenza dei parametri del secondo gruppo, detti didisturbo, indicati con .

Se = 0, ovvero se fosse noto, la funzione di verosimiglianza sarebbeL(, 0). Ci tipicamente non avviene e v quindi la necessit di introdurreuna sorta di surrogato per la funzione di verosimiglianza propria in cui sostituito con una sua stima. Tale funzione e la sua rispettiva trasforma-ta logaritmica sono chiamate di verosimiglianza e log-verosimiglianzaprofilo.

1.7 Metodo di Newton-Raphson 14

Lp() = L(, ) e lp() = l(, ) (1.7)

La stima per ottenuta come soluzione dellequazione l(,) = 0, ovveroponendo pari a 0 la derivata parziale in della funzione di verosimiglianzapropria considerando fissato. Per questo motivo evidente che la stima dimassima verosimiglianza profilo per coincide con quella di verosimiglianzapropria.

Viene ora presentata una propriet della verosimiglianza profilo di cui cisi servir in seguito. Definita linformazione osservata profilo come

jp() = 2

Tlp() =

2

Tl(, ),

si pu dimostrare che linversa di jp() pari allinversa della matrice diinformazione osservata complessiva, nel blocco (, ). Quindi se j(, ) unamatrice a blocchi definita come

j(, ) =

[j(, ) j (, )

j(, ) j (, )

],

allora vale che

jp()1 = [j(, )1] =

(j (, ) j(, )j(, )1j (, )

)1, (1.8)

il cui risultato deriva dalla regola di inversione a blocchi.Per sistemi dipotesi in cui si ha che H0 : = 0 e H1 : 6= 0, esiste una

statistica test detta di log-rapporto di verosimiglianza profilo. Si pudimostrare che essa sotto H0 converge in distribuzione ad una distribuzione2q, dove q = dim(0). La sua definizione

Wp(0) = 2(lp() lp(0)

)= 2(l( , ) l(0, 0)

).

1.7 Metodo di Newton-RaphsonAlcune difficolt sorgono quando si tenta di ottenere la soluzione di l() =

0. Pu accadere che non sia possibile esplicitare la stima di massima verosi-miglianza ed in tal caso si deve ricorrere ad algoritmi numerici. Uno di questi,


che non necessariamente il pi efficiente in termini computazionali, il me-todo di Newton-Raphson. un algoritmo iterativo che necessita di un puntodinizio. Ne viene data, in questo paragrafo, una sommaria spiegazione.

Lo sviluppo in serie di Taylor della funzione l() porge

l() = l(0) + l(0)( 0). (1.9)

Imponendo la soluzione l() = 0 e ricordando che l() = j() si ottieneche

l(0) j(0)( 0) = 0j(0)( 0) = l(0) = 0 + j(0)

1l(0),

da cui possibile stabilire un algoritmo iterativo in cui

k+1 = k + j(k)1l(k). (1.10)

1.7.1 Algoritmi numerici con R

Nel corso della relazione verr utilizzato il metodo di Newton-Raphsoncome approssimazione numerica della stima di massima verosimiglianza. bene tener conto che, tuttavia, esistono numerosi altri algoritmi che permet-terebbero di ottenere lo stesso risultato. Molti di questi sono gi implementatinel software R. Per maggiori approfondiment riguardanti il software utilizzatosi veda R Core Team (2012).

A titolo esemplificativo, vengono calcolate le stime di massima verosi-miglianza per i dati riportati in Tabella 1.1, che rappresentano il tempo dirottura di alcune molle, sottoposte a ripetuti sforzi tramite pesi differenti. Idati hanno la seguente struttura. A questo scopo, si utilizzeranno solamentele prime 10 osservazione relative alla variabile cycles che si assumer sia-no realizzazioni indipendenti provenienti da una variabile Weibull che abbiafunzione di densit pari a

p(x; , ) =(

)(x

)1exp{

(x

)}.


Tabella 1.1: Le prime 10 osservazioni dei dati stress

cycles cens stress1 225 1 9502 171 1 9503 198 1 9504 189 1 9505 189 1 9506 135 1 9507 162 1 9508 135 1 9509 117 1 95010 162 1 950

y


library(numDeriv)

theta.cappello

Capitolo 2

Modello lineare e trasformazione divariabili

2.1 IntroduzioneNel capitolo precedente si discusso in maniera molto generale di modelli

statistici parametrici. Una loro sottoclasse data dai modelli di regressio-ne lineare la cui importanza indiscussa in ambito econometrico, sociale,medico e non soltanto. L origine degli studi riguardanti i modelli lineari abbastanza incerta: le prime fondamentali nozioni vennero introdotte daAdrien-Marie Legendre e Carl Friedrich Gauss, nei primi anni del 1800. Es-se tuttavia si svilupparono notevolmente soprattutto quando si cominci adipotizzare una distribuzione per gli errori, estensione ad opera di Karl Pear-son e George Udny Yule tra il 1897 ed il 1903. Si veda Wikipedia, ModelloLineare (2013) e i documenti ivi citati.

Oggi i modelli lineari sono ampiamente utilizzati e studiati. Ne sono stateintrodotte numerose estensioni: una fra tutte sono i modelli lineari gene-ralizzati, che per non verranno trattati in questa relazione. Per una loroesaustiva trattazione si veda McCullagh e Nelder (1989). Queste innovazionisono state possibili grazie allevoluzione della tecnologia che ha permesso losviluppo di rapidi algoritmi di calcolo numerico.

Nei modelli lineari vengono fatte alcune ipotesi circa la natura delle os-

2.2 Il modello lineare 19

servazioni e spesso esse vengono contraddette in fase di verifica. In questicasi si pu scegliere di ignorare queste violazioni, se si ritiene che siano lievio irrilevanti, oppure modificare il modello di partenza, nel tentativo di cor-reggerle. Un metodo diffuso consiste nellagire sulla variabili, dipendenti oindipendenti che siano, trasformandole con opportune funzioni. Nel caso sidecidesse di trasformare la variabile dipendente, sono state fatte numeroseproposte di cui la soluzione di Box-Cox rappresenta una generalizzazione. Inquesto capitolo verranno introdotti alcuni concetti fondamentali riguardanti imodelli lineari e verr discusso lapproccio proposto da Box-Cox, riportandoalcuni esempi.

2.2 Il modello lineareSi definisce modello lineare una qualunque relazione esprimibile nella

formaY = X + , (2.1)

in cui X una matrice di dimensione n p nella quale ciascuna colonnarappresenta una variabile, un vettore di dimensione p 1 di parametriignoti ed un vettore stocastico di dimensione p 1. Le caratteristiche delmodello lineare quindi sono

la componente stocastica viene addizionata al resto del modello;

la funzione che esprime Y in funzione di X lineare nei parametri.

Un problema di regressione lineare consiste nel cercare di stimare i pa-rametri ignoti contenuti nel vettore a partire dalle realizzazioni di y =(y1, . . . , yn) che si suppone abbiano la struttura descritta nellequazione (2.1).

Senza ulteriori restrizioni la stima di praticamente impossibile ed per questo che si suppone che valgano almeno le ipotesi del secondordine.Si richiede cio che

X sia una matrice non stocastica di rango pieno: rank(X) = p;

la media degli errori sia nulla: E[] = 0;

2.2 Il modello lineare 20

la varianza degli errori sia costante (omoschedasticit) e questi sianoincorrelati tra loro: Var() = 2In.

Utilizzando il criterio dei minimi quadrati possibile giungere a stima-tori consistenti per che godono di ottime propriet. Tuttavia per condurreanalisi inferenziali bisogna ipotizzare una distribuzione per le osservazioni.Sotto l ipotesi di normalit si ha che distribuito come una normalemultivariata e quindi

Nn(0, 2In) che implica Y Nn(X, 2In).

2.2.1 La stima dei parametri

Le stime per , tramite il criterio dei minimi quadrati, si ottengono mi-nimizzando la funzione Q() = (Y X)T (Y X) = ||Y X||2, ciola distanza tra i valori assunti da Y e quelli previsti dal modello. Si pu di-mostrare che nella funzione di massima verosimiglianza profilo per , sottolipotesi di normalit, il punto di massimo assoluto coincide con quello di mi-nimo della funzione Q(). Entrambi i metodi di stima per portano dunqueal risultato

= (XTX)1XTY.

Sotto le condizioni del secondordine lo stimatore si verifica essere non di-storto e consistente ed inoltre denominato BLUE (Best Linear UnbiasedEstimator). Ci significa che esso lo stimatore pi efficiente nella classedegli stimatori lineari non distorti. Questo risultato stato dimostrato nelnoto teorema di Gauss-Markov. Si veda Pace e Salvan (2001, pag. 296).Una stima per 2

2 =||y X||2

n,

che la stima di massima verosimiglianza. In genere per, poich 2 unostimatore distorto, si preferisce una sua lieve correzione: s2, definita come

s2 =||y X||2

n p,

che uno stimatore corretto della varianza del modello.

2.3 Trasformazioni della variabile dipendente 21

2.2.2 Inferenza sui parametri

La semplice stima dei parametri sarebbe inutile se non si riuscisse inqualche modo a quantificarne la precisione. Perci si ottengono dei test chesono connessi alla verosimiglianza e sono per propri solo dei modelli linearinormali.

Si voglia verificare un sistema dipotesi del tipo H0 : H = 0 e H1 :H 6= 0 dove H una matrice q p, con q p. Ciascuna riga rappresentauno dei vincoli lineari che si intende imporre. Nel capitolo precedente statodefinito il test statistico log-rapporto di verosimiglianza che permetterebbedi verificare H0. La distribuzione normale rende i calcoli analitici piuttostoagevoli e questo consente di arrivare, tramite opportune trasformazioni mo-notone di Wp(), ad una quantit la cui distribuzione completamente nota.Si dimostra che

F =|| 0||2/q

||y ||2/(n p) F(q,np),

in cui = X e 0 = X0. Nel caso in cui q = 1 si ha che t2 = F ,con t che ha distribuzione t di student con n p gradi di libert. Questoimplica che con un unico test possibile verificare tutti i vincoli lineari che sidesiderano. Esistono alcuni rilevanti casi particolari: spesso si vuole verificarela nullit di un parametro o di un intero gruppo di questi. Se lipotesi nullavenisse accettata, si potrebbero escludere alcune variabili dal modello con unrisparmio dal punto di vista interpretativo e, a volte, anche economico.

Queste interessanti applicazioni si reggono sulle ipotesi sopra formulatee, se queste venissero a cadere, i test statistici potrebbero diventare privi disignificato. Ecco perch verranno presentate alcune soluzioni nel paragrafosuccessivo.

2.3 Trasformazioni della variabile dipendenteLe ottime propriet del modello lineare normale decadono se qualcuna

delle ipotesi iniziali non rispettata. Ad esempio, i dati potrebbero nonavere varianza costante (eteroschedasticit), oppure gli errori potrebbero

2.3 Trasformazioni della variabile dipendente 22

non provenire da una distribuzione normale o essere correlati. Non esiste ununico approccio risolutivo e le vie percorribili in generale sono:

ignorare la violazione delle ipotesi e trattare il modello come se questefossero vere;

inserire nel modello nuove variabili. bene far notare che spesso questonon possibile per motivi economici o di contesto;

rispecificare il modello. In presenza di errori correlati, ad esempio,potrebbe essere conveniente utilizzare modelli per serie storiche;

trasformare le variabili tramite una funzione

Non esiste una ricetta che permetta di scegliere la soluzione migliore, ammes-so che ne esista una sola. La scelta dipende dalle esigenze di chi poi utilizzeril modello e dalla gravit della violazione delle ipotesi.

Anche nel caso in cui si decidesse di trasformare una variabile, le possibili-t sono numerose. Si potrebbe agire contemporaneamente su tutte le variabilio su una sola di esse. Focalizzando lattenzione sulla variabile dipendente,le funzioni utilizzate pi di frequente sono: la radice quadrata, il logaritmonaturale, il reciproco.

Una prima classe di trasformazioni per la dipendente stata proposta inTukey (1957). Essa prevedeva

y =

y se 6= 0log y se = 0 , (2.2)per valori di y positivi. A seconda del parametro si possono ottenere unagrande variet di trasformazioni. Due raffinamenti della (2.2) sono dati dalletrasformazioni Box-Cox, in Box e Cox (1964), le cui espressioni sono

y =

y1

se 6= 0

log y se = 0, (2.3)

per valori di y strettamente positivi e

y =

(y+2)11

1se 1 6= 0

log (y + 2) se 1 = 0, (2.4)

2.4 Modalit dutilizzo 23

per valori di y maggiori di 2. Nel seguito della relazione si far riferimentoalla (2.3) a meno che non venga espressamente detto il contrario. Si noti cheessa pari alla (2.2) a meno di trasformazioni lineari. Questo implica che idue modelli sono del tutto equivalenti. Tuttavia vale lapprezzabile propriet

lim0

y 1

= log y,

che implica la continuit della funzione in .Molte altre trasformazioni sono state proposte nel corso degli anni, cia-

scuna per correggere difetti della trasformazione Box-Cox. Per una rassegnaesaustiva si veda Sakia (1992). Nonostante ci, la sua semplicit ha fatto inmodo che si diffondesse rapidamente e divenisse uno standard nellambitodella trasformazione di variabili.

Una formulazione alternativa di una certa rilevanza presente in Bickele Doksum (1981, pag. 297), riportata qui di seguito

y =sgn(y)|y| 1

, (2.5)

con > 0. Per gli y positivi la trasformazione, ed anche la funzione diverosimiglianza, coincide con la (2.3). Tuttavia ora y pu assumere anchevalori negativi. Questa trasformazione , da un punto di vista teorico, picorretta perch il campo di variazione di y ora lintero insieme reale.

2.4 Modalit dutilizzoLapproccio usuale consiste nello stimare il parametro , spesso tramite

verosimiglianza, per poi trattare la stima, o un suo intero vicino, come noto.Sono gli stessi Box e Cox a proporre questa modalit in Box e Cox (1964,pag. 239). Si deve per stabilire se dinteresse studiare il fenomeno in unadeterminata scala, suggerita dai dati, oppure lo si vuole analizzare in unascala ignota che dipende dal vero parametro 0, come loro stessi fanno nota-re. Nel secondo caso si sta quindi supponendo che il modello correttamentespecificato sia

Y 1

= X + . (2.6)

2.5 I dati cars 24

evidente che non si pu trascurare la variabilit aggiunta che la stimadi comporta. In particolare le varianze associate alle stime dei parametri saranno presumibilmente maggiori di quelle del modello in cui trattatocome noto. Si dovr quindi utilizzare un approccio alternativo per cercare diquantificare la variabilit aggiuntiva.

2.5 I dati carsVerr ora presentato un dataset per illustrare la procedura originaria

proposta da Box-Cox. Si scelto un insieme di dati in cui, ovviamente, vifosse la necessit di una trasformazione della variabile dipendente che fosse,inoltre, strettamente positiva.

I dati cars sono relativi alla distanza percorsa da un auto, che viaggiavaad una certa velocit prima di fermarsi. La distanza espressa in piedi, la ve-locit in miglia orarie. Sono 50 osservazioni che risalgono agli anni venti. Sonopresenti nel software R e possono essere facilmente richiamati col comandodata(cars). In Figura 2.1 v una loro rappresentazione. Si interessati aspiegare la distanza percorsa dalle auto in funzione della loro velocit.

Appare evidente che intercorre una relazione tra le due variabili. Si intra-vede tuttavia un problema di eteroschedasticit nel caso si scegliesse di utiliz-zare semplicemente un modello lineare, anche se in forma non eccessivamentegrave.

Si introduce allora la trasformazione Box-Cox e dunque il modello vienedefinito in questa maniera

Y i 1

= 0 + 1xi + i,

per i = 1, . . . , 50, in cui sono valide tutte le ipotesi fatte in questo capitolo.La stima di massima verosimiglianza = 0.43, ottenuta con i metodidescritti nei capitoli successivi. Al posto di utilizzare la stima di massimaverosimiglianza si assume che sia noto e posto pari a 0.5. Il modello diventa

Yi = 0 + 1xi + i, (2.7)

Utilizzando il software R il modello viene stimato e sono riportate inTabella 2.1 le stime dei coefficienti con i relativi scarti quadratici medi.

2.5 I dati cars 25

Figura 2.1: I dati cars

2.5 I dati cars 26

Figura 2.2: Log-verosimiglianza profilo per nei dati cars

2.5 I dati cars 27

Tabella 2.1: Sommario di R per la stima del modelloYi = 0 + 1xi + i

Estimate Std. Error t value Pr(>|t|)(Intercept) 1.2771 0.4844 2.64 0.0113

speed 0.3224 0.0298 10.83 0.0000

Ladattamento buono, non si riscontrano particolari problemi nei gra-fici diagnostici. Questi stessi dati verranno rianalizzati successivamente, uti-lizzando un approccio differente. Ci che per ora si fa notare che lunicavariabile esplicativa indubbiamente significativa.

Capitolo 3

Box-Cox e lapproccio di verosimiglian-za

3.1 IntroduzionePer condurre unanalisi approfondita circa gli effetti di sul resto del mo-

dello necessario ricavare le matrici di informazione osservata j() e quella diinformazione attesa I() da utilizzare, ad esempio, nel test statistico definitonella (1.5). Grazie a queste due matrici sar possibile quindi avere uneviden-za empirica a sostegno della supposizione dellaumento della variabilit nelmodello. La prima pu essere ottenuta sia analiticamente che numericamen-te. La seconda invece comporta alcune difficolt aggiuntive. Per ottenerle sidovr innanzitutto costruire un algoritmo di approssimazione numerica conil quale ottenere le stime di massima verosimiglianza. Si deciso di mostrareogni singolo passaggio effettuato per raggiungere il risultato. I risultati ine-renti alla matrice di informazione osservata e lalgoritmo di Newton-Raphsonsono parzialmente basati su Scott, 1999.

Si ipotizza che il modello abbia la forma specificata nella (2.6), che glierrori siano normali e che valgano le ipotesi del secondordine.

3.2 Funzione di verosimiglianza 29

3.2 Funzione di verosimiglianzaIl primo passo consiste nel ricavare la funzione di densit per Y , da cui poi

si otterr la funzione di verosimiglianza e tutte le quantit ad essa collegate.La distribuzione di Y non nota ma pu essere calcolata analiticamente. SiaY la variabile trasformata per la quale vale che

Y = X + .

Allora seguir che Y N(X, 2In) e la funzione di densit sar pari a

fy(y) =1

(2)n/2e

122

(yX)T (yX),

che si ricava dalla definizione di normale multivariata. Non questa tuttaviala quantit dinteresse: si vuole trovare la funzione di densit per Y . Si defi-nisca Y = g(Y ) la relazione che intercorre tra le due variabili. noto che lefunzioni di densit sono legate dalla relazione

fy(y) = fy(g(y)

)|J(y)|, (3.1)

dove |J(y)| il determinante della matrice jacobiana definita come J(y) =g(y)y

. Si ottiene che

g(yi)

yj= 0 per i 6= j e g(yi)

yj= y1j per i = j,

con i, j = 1, . . . , n. Quindi risulta che J(y) = diag(y11 , . . . , y

1n

). Trattan-

dosi di una matrice diagonale, il determinante semplicemente il prodottodi ciascun elemento:

|J(y)| =ni=1

y1i .

A questo punto, sfruttando la (3.1), si ottiene la funzione di densit e dunqueanche la funzione di verosimiglianza per Y . Viene riportata solamente laseconda che pari a

L(; y) =1

(2)n/2e

122

(yX)T (yX)ni=1

y1i ,

3.3 Funzione punteggio 30

dove = (, 2, ). Come di consueto, si preferir non trattare direttamentela funzione di verosimiglianza, ma la sua trasformazione logaritmica. Percisi proseguir utilizzando la funzione di log-verosimiglianza che

l(; y) = n2

log (22) 122

(y X)T (y X) + ( 1)ni=1

log yi.

3.3 Funzione punteggioLa funzione punteggio, definita nella (1.1), verr calcolata a gruppi di

parametri. Si procede innanzitutto alla derivazione per e per 2, il cuicalcolo abbastanza agevole. Derivando per si ottiene che

l(; y)

= 1

22(y X)T (y X)

= 122

(yT y (X)Ty yTX + TXTX)

= 122

(2XTy + 2XTX)

= 12

(XTy +XTX).

Inoltre la derivazione per 2 fornisce

l(; y)

2= n

22+

(y X)T (y X)2(2)2

.

Si nota subito che queste due funzioni sono molto simili alle corrispondentiquantit di verosimiglianza in cui, al posto della variabile y, c la sua tra-sformazione y. Ci comporta che si otterranno le consuete stime di massimaverosimiglianza per la variabile trasformata, per fissato.

Infine si deriva rispetto a , ottenendo

l(; y)

= 1

22

(y X)T (y X)

+ni=1

log yi

= 2(y X)T

22y

+ni=1

log yi

= uT (y X)

2+

ni=1

log yi,

3.4 Matrice di informazione osservata 31

dove u = y

. Il generico elemento ui per 6= 0 dato da

ui =

yi 1

=yi(log yi 1

)+1

2.

Per = 0 invece la funzione ricavata come soluzione del limite lim0 ui.Nel complesso risulta quindi che

ui =

yi (log y

i 1)+12

se 6= 0log2 (yi)

2se = 0

.

Riassumendo, la funzione score il vettore seguente

l(; y) =

l(;y)

l(;y)2

l(;y)

= 12

(XTy +XTX)

n22

+ (yX)T (yX)

2(2)2

uT (yX)

2+n

i=1 log yi

.

3.4 Matrice di informazione osservataLa matrice j() ricavata prima ancora di impostare le equazioni di vero-

simiglianza perch sar necessaria in seguito per il loro calcolo. Questa inver-sione nella procedura quindi, in questo caso, concettualmente pi linearedi quella usuale.

Come per la funzione l(; y), anche j() verr calcolata per gruppi diparametri. Si suddivide la matrice simmetrica j() a blocchi, come segue

j() =

j11() j12() j13() j22() j23() j33()

= l(;y)T

l(;y)2

l(;y)

l(;y)(2)2

l(;y)2

l(;y)2

. (3.2)Innanzitutto si ottiene il blocco per la derivata seconda di

j11() = l(; y)

T=

( 12

(XTy +XTX))

=XTX

2,

3.4 Matrice di informazione osservata 32

ovvero una matrice di dimensione p p. Il vettore misto in e 2

j12() = l(; y)

2=

2

( 12

(XTy +XTX))

=XTy XTX

(2)2,

di dimensione p 1. Il vettore misto in e

j13() = l(; y)

=

( 12

(XTy +XTX))

= XTu

2,

ricordando che si definito u = y

. Quindi j13() ha dimensione p 1. Loscalare relativo alla derivata seconda per 2 dato da

j22() = l(; y)

22=

2

( n

22+

(y X)T (y X)2(2)2

)= n

2(2)2+

(y X)T (y X)(2)3

.

Lo scalare relativo alla derivata mista in e 2 pari a

j23() = l(; y)

2=

((y X)T (y X)2(2)2

)= 2(y X)

T

2(2)2y

= uT (y X)

(2)2= (y X)

Tu

(2)2

Infine va calcolata la derivata seconda in

j33() = l(; y)

2=

uT (y X)2

=uTu+ (y X)Tv

2,

in cui v = u. La sua derivata allelemento i-esimo per 6= 0 data da:

3.5 Stime di massima verosimiglianza 33

vi =ui

=yi [(log (y

i )

2 2 log (yi ) + 2] 23

.

Come con ui il caso per = 0 definito come il lim0 vi = log (yi)3

3. Que-

sto risultato mostrato in appendice. Riassumendo, quindi, la matrice diinformazione osservata pari a

j() =

XTX2

XT yXTX(2)2

XTu2

n2(2)2

+ (yX)T (yX)

(2)3 (yX)

Tu(2)2

uTu+(yX)T v

2

.

3.5 Stime di massima verosimiglianzaNella ricerca della stima di massima verosimiglianza, che va cercata tra

le soluzioni dellequazione l() = 0, si incorre in un problema: non tuttele soluzioni sono ricavabili esplicitamente. In particolare, lequazione relativaalla derivata in non ammette soluzione esplicita. Per le altre si ricava che

= (XTX)1XTY, (3.3)

2 =Y T (In P )Y

n=

(Y X)T (Y X)n

. (3.4)

La matrice P di proiezione ed definita come P = X(XTX)1XT . Questosuggerisce che per trovare il vettore sufficiente massimizzare la funzionedi verosimiglianza profilo per e sostituire il punto di massimo nella (3.3) enella (3.4). La log-verosimiglianza profilo pari a

lp(; y) = n

2log (22)

1

22(y X)T (y X) + ( 1)

ni=1

log yi

che equivalente a

lp(; y) = n

2log (2) +

ni=1

log yi (3.5)

in cui la stima vincolata 2 data dalla (3.4). Per massimizzare la (3.5) sipossono utilizzare numerosi metodi numerici. Verr utilizzato lalgoritmo di

3.5 Stime di massima verosimiglianza 34

Newton-Raphson che verr poi confrontato con altri algoritmi numerici, perverificarne la correttezza.

Innanzitutto viene calcolata la derivata di lp(; y) che coincide con l(;y)in cui per i parametri e 2 sono sostituiti con la loro stima di massimaverosimiglianza. Si ottiene quindi che:

g() =lp(; y)

= u

T (y X)2

+ni=1

log yi (3.6)

La soluzione di g() = 0 non disponibile in forma chiusa, perci lalgoritmoiterativo sar:

k+1 = k + jp(k)1g(k) (3.7)

seguendo le indicazioni date nella (1.10).

3.5.1 La varianza per

La funzione jp() andrebbe calcolata come jp() =lp(;y)

2, ma poich

j() nota, essa si pu ricavare abbastanza velocemente usando la (1.8).Innanzitutto si fa notare che j12(, ) = j21(, )

T = 0, con = (, 2).Infatti:

j12(, ) =XTy XTX

(2)2=XTy XTX(XTX)1XTy

(2)2= 0

. Anche j22(, ) pu essere semplificata e risulta che

j22(, ) = n

2(2)2+

(y X)T (y X)(2)3

= n2(2)2

+n2

(2)3= n

2(2)2+

n

(2)2

=n

2(2)2

Perci una nuova suddivisione in blocchi d

j(, ) =

[j() j()

j() j()

]=

XTX2

0 XTu2

0T n2(2)2

(yX)Tu

(2)2

uTX2

(yX)Tu

(2)2uTu+(yX)T v

2,

3.6 Trasformazione Bickel e Docksum 35

e quindi complessivamente si ha che

jp()1 =

[j(, )

1]

=(j() j()j()1j()

)1.

Si vuole verificare il risultato di:

2j()j()1j() =

[uTX y

T (InP )u

2

] [(XTX)1 00T 2

2

n

][XTuy

T (InP )u

2

]

=[uTX(XTX)1 2y

T (InP )u

n

] [ XTuy

T (InP )u

2

]

= uTPu+ 2(yT (In P )u)2

n2.

Quindi nel complesso si ottiene che:

jp() =uTu+ (y X)Tv uTPu

2 2n

(yT (In P )u2

)2=uT (In P )u+ yT (In P )v

2 2n

(yT (In P )u2

)2. (3.8)

Ottenuto questo risultato bisogna affidarsi ad un calcolatore per ottenere lastima per . Un esempio di codice per il software R dato in appendice.

3.6 Trasformazione Bickel e DocksumFinora si scelto di ignorare una particolarit della trasformazione Box-

Cox: non si tenuto conto del fatto che gli errori variano nellintervalloreale mentre la trasformazione definita solo nei reali positivi. Nella maggiorparte dei casi questo ha poco senso, a meno che non si voglia sostenere che lavariabile Y possa assumere anche valori complessi, cosa che nei casi praticiavviene di rado.

Condizionandosi al caso in cui ciascuna realizzazione di Y positiva, siottengono risultati ragionevoli sia nelle stime di massima verosimiglianza,sia nella matrice di informazione osservata, pur trattandosi di un approccioformalmente poco corretto.

Qualora si intendesse calcolare il valore atteso E(Y ), per, ci si scontre-rebbe con questa incogruenza formale. Il valore che si otterrebbe sarebbe

3.6 Trasformazione Bickel e Docksum 36

complesso, E(Y ) C, e non sarebbe pi n utilizzabile n interpretabile.Per risolvere questa difficolt si sceglie di abbandonare la trasformazio-ne Box-Cox e di utilizzare al suo posto quella proposta da Bickel e Doksum,1981, definita nella (2.5).

Quanto descritto nei capitoli precedenti, tuttavia, non perde di validitnel caso in cui y > 0. La trasformazione Bickel-Docksum non deve esserevista come un approccio completamente differente, quanto come un esten-sione della trasformazione Box-Cox dato che, per y > 0, queste coincidono.Perci sia le stime di massima verosimiglianza che la stima della matrice di in-formazione osservata saranno le stesse. Nel caso generale, invece, la funzionedi log-verosimiglianza

lBD(; y) = n

2log (22) 1

22(y X)T (y X) + ( 1)

ni=1

log |yi|,

che quindi differisce da l(; y) solamente per un modulo e per il fatto chey ora rappresenta la nuova trasformazione. Il determinante della matricejacobiana

|J(y)| =ni=1

|yi|1.

I primi due blocchi della funzione score sono praticamente identici a quelliprecedentemente ottenuti. Vi sono alcune distinzioni per nel termine l(;y)

.

Si ricava infatti che

lBD(; y)

= u

T (y X)

2+

ni=1

log |yi|,

in cui, con conti analoghi ai precedenti, si definisce lelemento i-esimo delvettore u

ui =yi

=sign(yi)|yi|(log (|yi|) 1) + 1

2.

Non necessario ricavare la matrice di informazione osservata ex novo, ma sufficiente un accorgimento: jBD() si ottiene sostituendo u al vettore ue v al vettore v allinterno della matrice j(). Si definisce v = u , il cuielemento i-esimo

vi =ui

=sign(yi)|yi|[(log (|yi|)

2 2 log (|yi|) + 2] 23

.

3.7 Matrice di informazione attesa 37

3.7 Matrice di informazione attesaAnche se jBD() ha propriet apprezzabili, esistono validi motivi che spin-

gono alla ricerca della matrice I(). Innanzitutto I() definita positiva inqualunque punto essa venga calcolata, mentre questo non vero per jBD()che lo , in generale, solo nel punto di massima verosimiglianza. Inoltre I()1

raggiunge il limite inferiore di Cramer-Rao che implica, sostanzialmente,una maggiore efficienza dello stimatore. La dimostrazione e lillustrazione diquesto risultato sono date in Azzalini, 2008, pag. 79.

Come si visto, per, la trasformazione Box-Cox non si presta al calcolodei valori attesi. Qui di seguito quindi verr mostrato come ricavare la matricedi informazione attesa della trasformazione (2.5) basandosi parzialmente suirisultati gi ottenuti in Bickel e Doksum, 1981.

3.7.1 Calcolo degli elementi di I()

La matrice I() non pu essere calcolata tramite un processo numericoautomatico come avviene con jBD(). Si suppone anche questa volta che I()sia una matrice decomposta a blocchi

I() =

i11() i12() i13() i22() i23() i33()

, (3.9)e verr calcolata come valore atteso di ciascun elemento di jBD(). Per primacosa si otterranno i valori attesi degli elementi che non presentano grosseproblematiche computazionali. Lelemento per i coefficienti pari a

i11() = E[XTX2

]=XTX

2, (3.10)

che quindi coincide con quello della matrice di informazione osservata, i11() =j11(). Lelemento in e 2 pari a

i12() = E[XTY XTX

(2)2

]=XTX XTX

(2)2= 0, (3.11)

3.7 Matrice di informazione attesa 38

essendo E[Y] = X. In questo caso i12() = j12() ed inoltre anche i21() =j21(). Lelemento per 2 pari a

i22() = E[ n

2(2)2+

(Y X)T (Y X)(2)3

].

Per definizione, posto Q = ZTZ con Z Nn(0, In), si ha che Q 2n. Daci segue che se U Nn(0, 2In) allora UTU 22n. noto che E[Q] =n. Poich immediato notare che (Y X) Nn(0, 2In), si ottiene cheE[(Y X)T (Y X)] = n2. Perci

i22() = n

2(2)2+

n2

(2)3=

n

2(2)2. (3.12)

Quindi anche per questo caso si ha che i22() = j22(). Per gli elementisuccessivi la questione diventa pi complessa perch non si riesce a ricon-dursi a distribuzioni note di cui si conosce il valore atteso. Si ricorrer adalgoritmi numerici per approssimare i valori attesi che non sono ricavabilianaliticamente.

Si introduce un risultato preliminare, dimostrato in Ross, 2007, pag. 204.Sia Z una variabile casuale continua con funzione di densit f(z) e g() unafunzione a valori reali; allora vale che

E[g(Z)] =

+

g(z)f(z) dz. (3.13)

I valori attesi di cui non si dispone di un risultato analitico sono, sostanzial-mente, tre. In particolare si vogliono conoscere

i31() = E[u

TX

2

],

i32() = E[(Y X)

Tu(2)2

],

i33() = E[uT u + (Y X)Tv

2

].

(3.14)

Si procede quindi per integrazione numerica tramite la funzione integratedel software R, sfruttatando la (3.13). Si scelto di utilizzare V , definitaqui di seguito, come variabile di riferimento, la cui distribuzione nota.

3.8 Alcune considerazioni 39

necessario quindi esprimere ciascun termine in funzione di V . Si definiscequindi

V = sign(Y )|Y | che implica che |V | = |Y |. (3.15)

Ne consegue inoltre che V Nn(1+X, 22In), grazie alle propriet delladistribuzione normale. I vettori u e v possono essere riscritti come segue

u =1 + V (log |V | 1)

2,

v =V(log(|V |)2 2 log |V |+ 2

)2

3.

Ponendo pari a X si ottiene anche che

Y X =V 1 . (3.16)

Unendo queste quantit nella maniera indicata nella (3.14), si perviene alrisultato finale. Non vengono riportate le formule complessive per evitareinutili appesantimenti nella lettura. Per i32() pu essere utile ricordare che

E[u

T (y X)

2+

ni=1

log |yi|]= 0,

per le propriet del vettore score. Si ricava quindi che

i32() = E[(Y X)

Tu(2)2

]=

1

2

ni=1

E[log |Vi|]. (3.17)

3.8 Alcune considerazioniIn questo capitolo sono state analizzate due differenti trasformazioni: quel-

la proposta da Box-Cox e quella suggerita da Bickel e Docksum. Esse nonsono intercambiabili e non possono essere ipotizzate contemporaneamente.Un approccio sensato, per dati positivi, potrebbe essere: ipotizzare la tra-sformazione (2.5) per la variabile dipendente tenendo conto che i risultatiinferenziali che se ne otterrebbero utilizzando invece la (2.4), sarebbero i me-desimi. Infatti coinciderebbero sia le stime di massima verosimiglianza che lamatrice di informazione osservata. Tuttavia ora lecito utilizzare la matrice

3.9 Test log-rapporto di verosimiglianza 40

di informazione attesa senza incorrere in valori complessi, che sarebbero didifficile interpretazione ed inoltre pi complicati da gestire a livello di calcolonumerico.

Nonostante le dimensioni trascurabili di questi valori, non si vede luti-lit di rocambolesche approssimazioni quando una soluzione pi corretta facilmente disponibile.

3.9 Test log-rapporto di verosimiglianzaUn metodo alternativo per condurre verifiche dipotesi sui coefficienti

il test log-rapporto di verosimiglianza. La sua definizione

Wp(0) = 2(lp() lp(0)

). (3.18)

Il suo utilizzo richiede la risoluzione di un problema di ottimo vincolato inquanto nella funzione lp(0) i termini di disturbo devono essere massimizzati.Per questa ragione necessario il calcolo analitico della seguente quantit

min l(; y) sotto il vincolo H = C.

Il teorema dei moltiplicatori di Lagrange garantisce che le soluzioni del pro-blema di ottimo vanno ricercate tra le soluzioni del seguente sistema{

l(; y) = HT

H = C,

dove il vettore dei moltiplicatori di Lagrange. Si sono omesse la stime di2 e perch esse coincidono con quelle usuali, sostituendovi le nuove stimedi . Come si visto in precedenza, il risultato che se ne ottiene {

12

(XTy +XTX) = HTH = C

.

Focalizzando lattenzione sul primo termine si ottiene che

12

(XTX)1(XTy)1

2(XTX)1XTX = (XTX)1HT,

= (XTX)1(XTy) 2(XTX)1HT, = 2(XTX)1HT.


Sostituendo la quantit ottenuta nel vincolo si ottiene che

H( 2(XTX)1HT

)= C,

2H(XTX)1HT = C H

Definendo K =(H(XTX)1HT

)1, si ottiene quindi che = 1

2K(C H)

Risostituendo il termine nellequazione precedente si ottiene in conclusioneche

0 = + (XTX)1HTK(C H).

Le stime per 2 e per derivano di conseguenza e vengono calcolate la primaanaliticamente, la seconda numericamente in maniera analoga a quanto vistofinora.

Capitolo 4

Verifiche e simulazioni

4.1 IntroduzioneNel capitolo precedente sono state trattate le trasformazioni dal punto

di vista teorico. Scopo di questo capitolo invece presentare la funzioneBoxcox per il software R, con la quale si otterranno tutte le quantit diverosimiglianza viste finora. Questa funzione verr poi sfruttata per condurredelle simulazioni.

Una prima parte del capitolo sar dedicata a verificare la correttezza diquesta funzione, confrontandola, ove possibile, con altre librerie gi presen-ti allinterno del software. Una seconda parte sar volta a verificare sia laconsistenza che la normalit degli stimatori, tramite simulazioni. Infine ci sioccuper di mettere in luce quali possono essere le conseguenze che deriva-no dallassumere il parametro come noto. Inoltre, per fornire un esempiopratico, verranno ripresi i dati cars.

4.2 La funzione BoxcoxLa funzione, come riportata nel Codice A.1 a pagina 69, richiede come

argomento il modello di cui si vuole effettuare la trasformazione di variabile.Il suo ouput costituito da una lista che comprende:

4.3 Verifica della correttezza del codice 43

un dataframe, denonominato Coefficienti, contenente le informazionirelative ai coefficienti : stime, deviazioni standard, test statistici;

un dataframe, denonominato lambda, contenente le informazioni rela-tive al parametro : stima, deviazione standard, test statistico per laverifica di = 1;

un intervallo di confidenza per con livello di significativit , deno-minato Intervallo_lambda;

la stima di massima verosimiglianza per 2, denominata Stima_varianza;

un oggetto matrix contenente la stima della matrice di informazioneattesa, I(), denonimata Informazione_attesa;

un oggetto matrix contenente la stima della matrice di informazioneosservata, j(), denominata Informazione_osservata;

un grafico della log-verosimiglianza profilo per in cui segnalata lasua stima e un intervallo di confidenza, disattivabile ponendo comeargomento plot=FALSE.

Tutti i test statistici e gli intervalli di confidenza si basano sul test di Wald incui si utilizzata la matrice di informazione attesa I(). Alloccorrenza si puselezionare un livello di significativit diverso da quello predefinito, ovveroalpha=0.05, inserendolo come argomento. Un esempio di output riportatonel Codice 4.1.

4.3 Verifica della correttezza del codicePrima di iniziare a trarre una qualunque conclusione necessario since-

rarsi della bont del codice utilizzato. Per fare ci si scelto, in alcuni casi,di confrontare i risultati prodotti dalla funzione Boxcox con quelli che si ot-terrebbero utilizzando altre funzioni. Altre volte questo non sar possibileperch potrebbero non esistere in R delle funzioni analoghe. In questi casi sisceglier un approccio differente.


Codice 4.1: Esempio di output per la funzione Boxcox utilizzando il dataframe carsBoxcox(lm(dist~speed))

#$Coefficienti

# Stima StdError TestZ Pvalue

#(Intercept) 1.0466220 1.0536595 0.993321 0.32055357

#speed 0.5064258 0.2203969 2.297790 0.02157372

#$lambda

# Stima StdError_Wald Test Pvalue

#1 0.4305987 0.1208956 -4.70986 2.478867e-06

#$Intervallo_lambda

#[1] 0.1936476 0.6675497

#$Stima_varianza

#[1] 2.8362

#$Informazione_Attesa

# 17.62922 271.490 0.000000 -374.90744

# 271.49000 4663.987 0.000000 -6635.20490

# 0.00000 0.000 3.107894 -62.40306

# -374.90744 -6635.205 -62.403057 10831.22406

#$Informazione_Osservata

# 17.62922 271.490 0.000000 -375.18496

# 271.49000 4663.987 0.000000 -6624.82784

# 0.00000 0.000 3.107894 -62.33529

# -375.18496 -6624.828 -62.335291 10797.92233


4.3.1 Stima per

La stima di massima verosimiglianza, nella funzione Boxcox, ottenutatramite lalgoritmo di Newton-Raphson, seguendo il metodo illustrato nelcapitolo precedente. Si deve confrontare quindi la stima per , calcolata conquesto metodo, con quella che si ottiene tramite la funzione nlminb, nellaquale si massimizza la log-verosimiglianza profilo. Se le due stime per fos-sero uguali, conciderebbero di conseguenza anche quelle per e 2 che sono,come si visto, calcolate esplicitamente. Per far ci si definita la funzio-ne nlogprofilo, che la funzione di log-verosimiglianza profilo cambiata disegno: lp(; y), presentata nel codice Codice A.2 a pagina 74.

Si condotta una simulazione: sono stati generati 1000 differenti modelli(supponendo = 0), in cui n = 50. Per ciascuno di questi modelli si stimato tramite i due differenti metodi allocando i risultati nei vettori lambda1 elambda2. Viene quindi analizzato il vettore |12|, di cui si forniscono alcunestatistiche descrittive. I comandi per la simulazione sono dati nel Codice A.3a pagina 75. La Tabella 4.1 riporta i risultati ottenuti. Eccezion fatta per

Tabella 4.1: Statistiche descrittive per |1 2|,

Min. 1st Qu. Median Mean 3rd Qu. Max.0.00000 0.00000 0.00000 0.00003 0.00000 0.00709

alcuni valori anomali, i due metodi portano agli stessi risultati. Lindicatorepi rilevante il terzo quartile di |1 2|, pari esattamente a 0 che indicaun errore praticamente nullo per la gran parte dei casi.

Utilizzando i dati cars, invece, si possono confrontare i grafici in Figu-ra 4.1, prodotti dalla funzione Boxcox e dalla funzione boxcox della libreriaMASS. Si tenga presente che la prima costruisce intervalli di confidenza ba-sandosi sul test di Wald, la seconda sul test log-rapporto di verosimiglianzae perci questi non coincidono, pur essendo molto simili.

4.3.2 Matrice di informazione osservata

Di interesse maggiore la verifica della correttezza della versione anali-tica della matrice j(). A tale scopo, si pu utilizzare il comando hessian


Figura 4.1: Grafico funzioni Boxcox e boxcox, utilizzando il dataframe cars

della libreria numDeriv per un confronto numerico. necessaria, per, qual-che operazione aggiuntiva perch deve essere definita una nuova funzione,quella di verosimiglianza, che sar poi utilizzata come argomento della fun-zione hessian. Si deciso di definirla di segno invertito per poter ottenereautomaticamente la matrice di informazione osservata. La funzione di vero-simiglianza cambiata di segno, nlog.ver, definita nel Codice A.4 a pagina75. Le stime di massima verosimiglianza sono ottenute anchesse tramite cal-colo numerico. Per confrontare i risultati si scelto di utilizzare la seguentequantit

tr(

(j1() j2())2), (4.1)

che rappresenta la somma di ciascun elemento della matrice differenza presoal quadrato. Vengono generati 1000 modelli, a cui viene applicata la (4.1), edil risultato viene allocato nel vettore quadJ. Questo processo eseguito nelCodice A.5 a pagina 75. Le statistiche descrittive, riportate in Tabella 4.2,

Tabella 4.2: Statistiche descrittive per il vettore quadJ,

Min. 1st Qu. Median Mean 3rd Qu. Max.0.00000 0.00000 0.00000 0.02030 0.00000 16.46000


confermano che, ad eccezione di qualche errore di approssimazione numeri-ca, le due stime portano circa allo stesso risultato. Anche in questo caso ilterzo quartile nullo ed indica che nella gran parte dei casi le stime sonopraticamente coincidenti.

Confronto matrici dei dati cars

Viene ora riportato un esempio, in cui si confrontano le matrici ottenutenei due metodi, utilizzando i dati cars. La prima si ottiene semplicementecon i comandi seguenti

attach(cars)

Boxcox(lm(dist~speed))$Informazione_Osservata

La seconda invece si ottiene tramite i comandi

theta.cappello

4.4 Normalit degli stimatori 48

Tabella 4.3: Matrice j() per i dati cars

(a) Metodo analitico

(Intercept) 1 2 17.62922 271.49000 0.00000 -375.18496271.49000 4663.98658 0.00000 -6624.827840.00000 0.00000 3.10789 -62.33529

-375.18496 -6624.82784 -62.33529 10797.92233

(b) Metodo numerico

(Intercept) 1 2 17.62922 271.49002 0.00000 -375.18505271.49002 4663.98702 0.00008 -6624.829420.00000 0.00008 3.10790 -62.33543

-375.18505 -6624.82942 -62.33543 10797.92750

con Nn(0, In), ovvero si sta implicitamente assumendo = 0. Si sono ge-nerate 2000 osservazioni indipendenti ed identicamente distribuite. I comandiper questa simulazione sono

set.seed(11)

x


(c) Matrice di informazione attesa

(Intercept) 1 2 1854.58 4669.65 0.00 -9165.764669.65 15592.38 0.00 -32989.32

0.00 0.00 859.87 -4704.75-9165.76 -32989.32 -4704.75 109022.59

(d) Matrice di informazione osservata

(Intercept) 1 2 1854.58 4669.65 0.00 -9165.354669.65 15592.38 0.00 -33055.40

0.00 0.00 859.87 -4704.77-9165.35 -33055.40 -4704.77 109668.21

Tabella 4.4: Risultati della simulazione, convergenza delle matrici I() e j()

Il processo generatore dei dati

y 1

= 0.5x+ , (4.3)

con Nn(0, In) e = 0.5. La quantit pivotale di cui nota la distribuzione la seguente

1 0.5I()11

N(0, 1). (4.4)

Il coefficiente stimato per N = 1000 volte, con n = 200, allocando ciascunvalore nel vettore beta.stima mentre gli standard error vengono conservatinei vettori SeJ e SeI. Nel Codice A.7 a pagina 76 si effettua questa specifi-ca simulazione. In Figura 4.2 vi una diagnostica grafica per verifica dellanormalit della quantit pivotale oggetto danalisi. Si applica inoltre il testdi Shapiro-Wilks, il cui p-value risulta essere circa 0.06. Alla luce di ci, siaccetta lipotesi di normalit dello stimatore nonostante lasimmetria mani-festata dai dati. Si deve infatti tenere conto che la normalit garantita soloasintoticamente e non per campioni finiti.

Per meglio comprendere quanto sia plausibile assumere la normalit dellostimatore, questa simulazione viene replicata cambiando alcuni parametri.


(a) Il grafico quantile contro quan-tile

(b) Frequenze e approssimazionegaussiana

Figura 4.2: Grafici diagnostici per la normalit dello stimatore con N = 1000 en = 200

Si assumer n = 50, 100, 150, 200 e = 1, 0, 0.5, 1. Nella Tabella 4.5 vieneriportato lerrore medio per 1 definito come

Errore Medio =N

i=1 1iN

1. (4.5)

I grafici diagnostici per queste simulazioni sono riportati in Figura 4.3. Ciche se ne deduce che, nella maggior parte dei casi, la convergenza allanormalit non avviene. Lo stimatore pare essere distorto anche se convergelentamente al vero valore, allaumentare della numerosit campionaria. Lavelocit di convergenza dipende dal valore dei parametri e va valutata casoper caso. Ci suggerisce una certa cautela al momento della verifica dipotesisoprattutto quando non si dispone di una grande quantit di dati.

= 1 = 0 = 0.5 = 1n = 50 -0.36 0.023 0.0015 0.0444n = 100 -0.19 0.011 0.0002 0.0264n = 150 -0.09 0.009 0.0001 0.0204n = 200 -0.06 0.009 0.0001 0.0112

Tabella 4.5: Errore Medio per 1


(a)=1

(b)=

0

(c)=

0.5

(d)=

1

Figu

ra4.

3:Grafic

idiagn

ositicip

erla

norm

alitdi1pe

rn=

50,100,150,200

insensoorario

apa

rtireda

llaprim

ain

alto

asinistra

4.5 Test log-rapporto di verosimiglianza profilo per 52

4.5 Test log-rapporto di verosimiglianza profilo per

Per ulteriore conferma della bont dei test, si implementato anche iltest log-rapporto di verosimiglianza profilo per che necessita delle funzioniBoxcox e nlogprofilo. La funzione data nel Codice A.6 a pagina 76. Unsemplice esempio dato dal confronto tra gli intervalli di confidenza ottenuticol test Wp() e quelli alla Wald. I risultati non dovrebbero coincidere esat-tamente, ma essere quantomeno simili. Nei dati cars, infatti, per = 0.05,si ha che

Test Wald IC1 = [0.1936; 0.6675],

Test Wp IC2 = [0.2203; 0.6696].

Si noti anche che il primo risulta essere simmetrico rispetto a per costru-zione, ma ci non si verifica anche per il secondo.

4.6 Simulazioni ed intervalli di confidenzaVengono ora fatte delle simulazioni per mostrare lentit della distorsione

sulle deviazioni standard dei coefficienti nel caso si assuma noto. Se laspecificazione corretta, infatti, lintervallo di confidenza con = 0.05 peril coefficiente j dovrebbe contenere il vero valore il 95% dei casi.

Vengono simulati N = 1000 modelli in cui si utilizza unestensione dellavariabile speed che assume numerosit campionaria pari a 400. Si artifi-ciosamente aumentato lampiezza del campione per permettere una miglioreconvergenza degli stimatori alla distribuzione normale. Si assume che i verivalori siano pari alle stime di massima verosimiglianza, ottenute in prece-denza per i dati cars, e che la specificazione sia corretta. In altre parole siassume che

dist 1

= 0 + 1speed + . (4.6)

Nel Codice A.8 a pagina 77 sono mostrati i comandi che generano i modelli,calcolano la stima per 1 e per la sua deviazione standard.

4.6 Simulazioni ed intervalli di confidenza 53

La differenza tra le deviazioni standard sistematica: la variabilit au-menta mediamente di 4 volte. In Figura 4.4 data una rappresentazione diqueste quantit. Vengono ricavati gli intervalli di confidenza per entrambi gli

(a) Approccio con supposto noto (b) Approccio con ignoto

Figura 4.4: Istogramma della deviazioni standard di 1

approcci: supponendo noto o lasciandolo ignoto. Si controlla, poi, quantevolte effettivamente il vero valore 1 contenuto nellintervallo. I comandiutilizzati sono

Conf.int1


sd(Coeff) mean(Se2) mean(Se1)

0.0668779 0.0670876 0.0137689

Tabella 4.6: Statistiche descrittive per la simulazione

deviazioni standard, mean(Se2). Seguendo vie differenti si giunti al mede-simo risultato: sd(Coeff) rappresenta una stima bootstrap della deviazionestandard del coefficiente. un metodo altrettanto valido ma che comportasforzi computazionali maggiori.

Si gi detto che nella gran parte dei casi si preferisce stimare il modellotramite unapprossimazione di allintero pi vicino. Verr ricondotta lastessa simulazione aggiungendo il comando

lambda.stima

4.8 Rianalisi dei dati cars 55

Wp.beta riportata nel Codice A.10 a pagina 78 che consente il calcolo deltest in questione e gli intervalli di confidenza. A sua volta, questa funzionenecessita della logprofilo.beta, riportata nel Codice A.9 a pagina 78 che,come il nome suggerisce, la log-verosimiglianza profilo per .

Si condotta una simulazione assolutamente analoga alla precedente alloscopo di evidenziare le eventuali differenze da quanto fatto pocanzi. Il Codi-ce A.11 a pagina 79 contiene i passaggi per esteso. Non sorprendentemente,i risultati portano a considerazioni molto simili: i nuovi intervalli di confi-denza comprendono il vero valore di il 94.6% delle volte. Migliore invece la distribuzione della statistica sotto lipotesi nulla gi utilizzando n = 30.Sotto H0 infatti vale che

Wp() 21,

poich un solo parametro stato vincolato. Nella Figura 4.5 v una diagno-stica grafica per questa ipotesi. Complessivamente, quindi, entrambi i testconcordano nel suggerire un aumento della variabilit dei coefficienti. Tutta-via il test di Wald pi facilmente interpretabile mentre il test di log-rapportodi verosimiglianza converge alla distribuzione nulla molto pi rapidamenteed quindi pi affidabile per i piccoli campioni.

4.8 Rianalisi dei dati carsVengono ora confrontati i due differenti approcci: quello utilizzato di con-

sueto, in cui si assume noto e pari alla stima di massima verosimiglianza,e quello in cui invece ignoto.

Codice 4.2: Comandi per la stima del modello, approcci differentilambda


Figura 4.5: Diagnostica grafica per la verifica della distribuzione 21 con N = 500,n = 30

(a) Output in cui supposto noto

Estimate Std. Error t value Pr(>|t|)(Intercept) 1.0466 0.7553 1.39 0.1723

speed 0.5064 0.0464 10.91 0.0000

(b) Output in ignoto

Estimate Std. Error z value Pr(>|z|)(Intercept) 1.0466 1.0536 0.9933 0.3205

speed 0.5064 0.2204 2.2978 0.0216

Tabella 4.7: Le stime del modello tramite i due metodi


di confidenza e soprattutto una differente statistica test. In generaleci che si verifica un aumento della variabilit delle stime. In questo casoparticolare, tuttavia, le conclusioni che se ne traggono sono le medesime. En-trambi gli approcci suggeriscono di includere la variabile speed nel modello,anche se a differenti livelli di confidenza.

4.8.1 Intervalli di confidenza per la media e di previsione

I due approcci conducono alla stessa stima puntuale ma sono differenti,tuttavia, anche gli intervalli di confidenza per i valori medi e gli intervalli diprevisione. Nei modelli lineari, si ha che

Y = PY Nn(X, 2P ), che implica che Yi N(i, 2hi),

da cui si ricavano gli intervalli di confidenza per il valor medio, calcolatipunto per punto. Il vettore h la diagonale della matrice di proiezione ed il vettore dei punti leva. Per gli intervalli di previsione si ha che

Y = Y + Nn(Y , 2In), che implica che Yi N(i, 2(1 + hi)).

Se si suppone che il modello dipenda da un valore ignoto di si ha innanzi-tutto che

Np(, I()

),

che deriva dalla teoria della verosimiglianza. Perci poich Y = X, vale larelazione

Y Nn(X,XI()XT

),

che permette di costruire intervalli di confidenza che tengano conto dellavariabilit di . Analogamente a prima, gli intervalli di previsione si costrui-scono a partire da

Y = Y + Nn(X, 2Ip +XI()

XT).

Sostituendo le quantit ignote con le opportune stime, si sono ottenuti in-tervalli di confidenza e di previsione per entrambi gli approcci, riportati inFigura 4.6. I risultati ottenuti sono completamente diversi e perci se ne con-clude che non indifferente lutilizzo di un metodo rispetto allaltro, comeillustrato in questo esempio. Si fa notare, comunque, che gli intervalli sonostati calcolati per la variabile trasformata e non nella scala originaria.


(a) Approccio in cui supposto noto

(b) Approccio in cui ignoto

Figura 4.6: Intervalli di confidenza e di previsione per i dati cars, = 0.05


4.8.2 Ulteriori analisi

Nella gran parte dei casi la trasformazione non avviene scegliendo parialla stima di massima verosimiglianza ma si sceglie, invece, un valore facil-mente interpretabile. Nella (2.7) si scelto infatti la radice quadrata e nonla stima di massima verosimiglianza. In questi casi il calcolo della matrice diinformazione osservata in un punto differente da quello di massimo pu por-tare a risultati privi di senso come ad esempio varianze negative. La matricedi informazione attesa, invece, definita positiva in qualunque punto essavenga calcolata. Perci si utilizza il comando

Boxcox(lm(dist~speed),lambdanoto=0.5)

che consente di calcolare il modello in cui = 0.5 ma utilizzando la matricedi informazione attesa che tenga conto della sua variabilit. Il modello ri-sultante equivalente, a meno di trasformazioni lineari, a quello stimato nelparagrafo 2.5 a pagina 29. Il risultato, per quel che riguarda i coefficienti, riportato in Tabella 4.8. I risultati, anche in questo caso, differiscono da

Tabella 4.8: Output modello in cui arrotondato allintero pi vicino

Estimate Std. Error z value Pr(>|z|)(Intercept) 0.55410 1.39444 0.39736 0.69110

speed 0.64483 0.26926 2.39483 0.01663

quelli che si otterrebbero supponendo = 0.5. Le deviazioni standard sonomaggiori e viene nuovamente confermato il fatto che la variabilit di non un dettaglio trascurabile. Non un caso, inoltre, che la stima di 1 coincidacon quella calcolata nella (4.7). Se il modello generatore dei dati fosse quellodescritto nella (4.6), si tratterebbe di una stima distorta.

Capitolo 5

Conclusioni

Lattenzione di questa relazione era posta sulla modalit con cui, normal-mente, si utilizza la trasformazione Box-Cox. Trattare il parametro comenoto un procedimento formalmente corretto solo nel caso in cui si stia cer-cando una differente scala da usare come riferimento. La stima del parametro solo un suggerimento che indica quale sia la trasformazione pi opportuna.

Nella maggioranza dei casi, tuttavia, si sta effettivamente supponendoche i dati provengano da una scala ignota che dipende da . Perci nontener conto del fatto che questo parametro stato stimato comporta delledistorsioni. La dimostrazione teorica di questo fatto era gi stata presentatain precedenza ed stata, in questa relazione, solamente ripresa.

Una percezione diffusa che, se anche fossero presenti delle distorsioni,queste sarebbero irrilevanti. Tramite alcune simulazioni, confrontando gli in-tervalli di confidenza per le stime dei parametri, si giunti alla conclusioneche queste sono invece piuttosto consistenti. Le conseguenze possono esse-re anche drammatiche. Spesso inoltre si utilizza un arrotondamento di aldecimale vicino e questo comporta in aggiunta una distorsione nelle stime.

Comunque si deve far notare che lentit del danno dipende dallo speci-fico caso in analisi e che le tecniche adottate per affrontare il problema sonodifficilmente utilizzabili nella quotidianit. A volte infatti pu essere utilesacrificare la correttezza formale per poter ottenere dei risultati immedia-tamente interpretabili e di facile comprensione. Ci avviene, ovviamente, incontesti diversi da quello accademico.

Bibliografia

Azzalini, A. Inferenza Statistica. Milano: Springer Verlag, 2008.Bickel, P.J. e A. Doksum. An analysis of transformations revisited. In:

Journal of the American Statistical Association 76 (1981), pp. 296311.Box, G.E.P. e D.R. Cox. An analysis of transformation. In: Journal of the

Royal Statistical Society. B (1964), pp. 211243.McCullagh, P. e J. A. Nelder. Generalized linear models. London: Chapman

e Hall, 1989.Pace, L. e A. Salvan. Introduzione alla statistica II. Milano: Cedam, 2001.R Core Team. R: A Language and Environment for Statistical Computing.

ISBN 3-900051-07-0. R Foundation for Statistical Computing. Vienna,Austria, 2012. url: http://www.R-project.org/.

Ross, S. M. Calcolo delle probabilita. Milano: Apogeo, 2007.Sakia, R. M. The Box-Cox Transformation Technique: A Review. In: Jour-

nal of the Royal Statistical Society. B (1992), pp. 169178.Scott, H. A likelihood based inference on the Box-Cox family of transforma-

tion: SAS and Matlab programs. Montana State University, 1999.Tukey, J. W. On the comparative anatomy of transformation diagnostics.

In: Annals of Mathematical Statistics (1957), pp. 602632.Wikipedia, Modello Lineare. 2013. url: http://it.wikipedia.org/wiki/

Regressione_lineare.

http://www.R-project.org/http://it.wikipedia.org/wiki/Regressione_linearehttp://it.wikipedia.org/wiki/Regressione_lineare

Appendice A

Codice R utilizzato

Codice A.1: La funzione BoxcoxBoxcox

63

sigma.stima

64

J12

65

}

for(j in 1:n) {

up

66

}

I11

67

if(sigma.stima.profilo

68

nlogprofilo

70

lambda

71

Coeff[i]

72

Wptest

Trasformazione Box-Cox: un’analisi basata sulla ...tesi.cab.unipd.it/42884/1/Rigon_Tommaso.pdf ·...

Documents

Transcript of Trasformazione Box-Cox: un’analisi basata sulla ...tesi.cab.unipd.it/42884/1/Rigon_Tommaso.pdf ·...