Intervalli di confidenza Bootstrap: una veduta dkinsieme e una

Intervalli di con…denza Bootstrap:
una veduta d’insieme e una proposta per un indice di
cograduazione
ANDREA BONANOMI
Dipartimento di Scienze Statistiche
Università Cattolica S.C., Milano
ABSTRACT: This paper gathers the di¤erent bootstrap methodologies
used to determine non-standard con…dence intervals. The di¤erent formulations seem to result very useful and e¢ cient both when standard intervals
cannot be calculated and because of the advantages in robustness recorded
in the estimation phase. The di¤erent methods (Bootstrap-t, Bootstrap Percentile, BCa, ABC) are presented in detail, considering both their parametric
and non-parametric structure and evaluating the inferential properties. Finally, the di¤erent approaches to determine bootstrap con…dence intervals
were examined for a cograduation index in order to evaluate the propensity
of a measurement scale (generally an ordinal scale) to assign a mainly positive
or negative judgment compared to another scale having the same number of
response alternatives.
Keywords: Bootstrap Con…dence Intervals, Bootstrap-t, Bootstrap percentile, BCa, ABC, Cograduation Index.
1
Introduzione
Gli Intervalli di Con…denza (IC) sono strumenti abitualmente utilizzati nella
statistica metodologica ed applicata in quanto combinano la teoria puntuale
della stima e la veri…ca d’ipotesi in un singolo strumento inferenziale di
grande impatto intuitivo.
Questo lavoro si propone di discutere la costruzione di IC con tecniche
alternative di ricampionamento, tipicamente bootstrap, sia quando è noto
il modello probabilistico sottostante (caso parametrico) sia quando non lo
è (caso non parametrico). Vengono inoltre confrontate tali tecniche nella
1
costruzione di un IC di un indice di cograduazione proposto per la valutazione della propensione di una scala psicometrica ad assegnare giudizi maggiormente positivi oppure negativi rispetto ad un’altra caratterizzata dallo
stesso numero di passi di risposta.
Siano X1 ; X2 ; :::; Xn variabili casuali (v.c.) indipendenti e identicamente
distribuite alla variabile casuale di campionamento X, con incognita funzione
di ripartizione F . Sia bn = b (x1 ; :::; xn ) uno stimatore del parametro di
interesse = (F ), con varianza 2n = var bn e b2n una sua stima. Se Cn =
Cn (x1 ; x2 ; :::; xn ) è un sottoinsieme di < dipendente solo da (x1 ; x2 ; :::; xn ) e
Pr ( 2 Cn )
1
(1)
dove
(0 <
< 1) è una pre…ssata costante, allora Cn è un intervallo
di con…denza per ad un livello 1 – . Pr ( 2 Cn ) è detta “probabilità di
copertura” di Cn . La maggior parte degli IC sono approssimabili, seppur
esistano un ristretto gruppo di intervalli esatti per situazioni particolari;
l’approssimazione tipicamente più utilizzata è rappresentata dall’intervallo
standard.
Se la variabile casuale standard
bn
n
si distribuisce come una variabile casuale Y
N (0; 1), allora l’intervallo di
con…denza standard per ad un livello di con…denza 1 –2 è dato da
bn
z(
)
n
(2)
con z ( ) 100 –mo percentile della v.c. normale standardizzata.
Qualora non si conosca il vero valore n , ma solamente un suo stimatore
bn , allora
bn
(3)
bn
si distribuisce asintoticamente come una v.c. normale standardizzata (e più
in generale come una v.c. t di Student, da cui il nome di (3) di variabile
casuale studentizzata) e l’intervallo
bn
z ( ) bn
2
(4)
rappresenta l’IC approssimato di con livello di con…denza nominale 1 –2
e probabilità di copertura determinata dalla (1).
Gli intervalli standard, implementati solitamente con l’utilizzo della teoria della massima verosimiglianza, hanno indubbi vantaggi, dettati da una
metodologia semplice e completamente automatica, per lo più in ambito parametrico. I problemi con l’utilizzo di IC standard derivano dal fatto che sono
basati su approssimazioni asintotiche che nella pratica, sovente, possono rivelarsi inappropriate. In Di Ciccio –Efron (1996) si dimostra che l’intervallo
(4) può essere anche sensibilmente di¤erente dall’intervallo esatto, in quei
casi nei quali l’intervallo esatto esiste.
Gli intervalli di con…denza bootstrap1 , presentati nelle loro diverse formulazioni nel seguito, appaiono particolarmente utili per ovviare ai problemi
tipici dell’analisi statistica “tradizionale”, in particolare in assenza di distribuzione parametrica nota e nel caso in cui si abbia a disposizione un
campione di numerosità limitata, tenuto comunque presente il fatto che anche in situazioni comunemente parametriche i vantaggi di robustezza delle
stime sono evidenti.
Un confronto con gli IC standard mostra che gli intervalli bootstrap sono
non solo asintoticamente più accurati, ma anche più corretti2 . In quei problemi nei quali esistono gli esatti limiti di con…denza, essi sono tipicamente
espressi nella forma
bn
A(
bn z ( ) + p
)
n
+
B( )
+ :::
n
(5)
dove n indica la dimensione campionaria. Gli IC standard sono corretti
del primo ordine, nel senso che (5) è dominata asintoticamente dal termine
bn bn z ( ) . Tuttavia, il termine del secondo ordine bn Ap( ) ha e¤etto rilevante
n
nelle situazioni di dimensioni campionarie ridotte. Gli IC bootstrap, come
verrà mostrato nel seguito, generano automaticamente una correttezza del
( )
secondo ordine, nel senso che sono espressi nella forma bn bn z ( ) + An
e, con particolari correzioni ed accorgimenti, anche una correttezza del terzo
ordine.
Nel seguito vengono illustrati, nel dettaglio, le più importanti tipologie
di costruzione di intervalli di con…denza bootstrap.
1
Per una veduta d’insieme della tecnica bootstrap si veda Efron B, Tibshirani R. (1993),
An Introduction to the Bootstrap, Chapman and Hall, New York.
2
Per la de…nizione di correttezza di un IC si veda Piccolo (1985), pp. 755-757
3
2
Intervalli di con…denza Bootstrap–t
L’uso di tecniche di ricampionamento di tipo bootstrap porta alla costruzione
di intervalli accurati (dal punto di vista statistico) senza l’assunzione di normalità distributiva, tipica degli intervalli di con…denza tradizionali. In questo
paragrafo viene presentato un primo approccio per la determinazione di IC
bootstrap, estremamente semplice, ma che gode di proprietà ottimali di e¢ cienza.
Sia X = (x1 ; x2 ; :::; xn ) una realizzazione campionaria dall’incognita distribuzione F con componenti i.i.d.. Sia Fb una stima di F (sia parametrica
che empirica); sia, quindi, b (X) uno stimatore dell’incognito parametro di
interesse e b (X) della deviazione standard di b.
Per analogia con la statistica t-student, si de…nisca
T =
b
(6)
b
e sia T ( ) il 100 –mo percentile di T . Il limite superiore di un intervallo di
con…denza ad una coda di livello per è dato da
b
b T (1
)
:
Tale formulazione assume l’ipotesi che il percentile di T sia noto, come nel
caso t-student usuale, dove T ( ) è il percentile di una distribuzione t. Tuttavia
i percentili di T non sono noti nella maggior parte delle situazioni.
L’idea dell’approccio Bootstrap–t è di stimare i percentili di T attraverso
una procedura di ricampionamento bootstrap.
Il primo passo consiste nella stima della distribuzione di campionamento
X e della generazione di B campioni bootstrap X = (X1 ; X2 ; :::; XB ). Ogni
Xb conduce in maniera naturale alla versione bootstrap dello stimatore b ,
con deviazione standard b , portando alla quantità pivotale
T =
b
b
b
replicazione bootstrap di (6). Un elevato numero B di replicazioni (in genere
maggiore di 1000) portano alla stima dei percentili Tb( ) = B –mo valore
ordinato di fT (b) ; b = 1; 2; :::; Bg. Il percentile –mo di T (b) è quindi
4
determinato dal valore Tb(
)
tale che
n
# T (b)
Tb(
B
)
o
= :
[Es: se B=2000, =0.95, allora Tb(0:95) è il 1900–mo percentile delle quantità
pivotali T (b) ordinate].
L’intervallo di con…denza Bootstrap–t è perciò determinato
h
i
b Tb(1 ) b; b Tb( ) b
con < 21 : La procedura di determinazione degli IC bootstrap–t è semplice e
completamente automatizzata con i software statistici di programmazione più
comuni, e di intuitiva interpretazione. Presenta però una serie di problemi:
Per utilizzare tale metodo, si necessita di uno stimatore b2n della varianza. Se non è disponibile alcuno stimatore, è possibile utilizzare uno
stimatore bootstrap vbB per b. Tuttavia se sia Tb( ) che vbB devono essere
approssimati attraverso simulazione Monte Carlo, occorre concatenare
le due simulazioni bootstrap; in sostanza, se si generano B1 campioni
bootstrap per approssimare Tb, per ogni replicazione occorre generare
B2 campioni bootstrap per vbB . Il totale di replicazioni cresce pertanto
velocemente, risultando pari a B1 B2 . Il linguaggio di programmazione
scienti…co R.2.3.1 ha implementato, nei suoi pacchetti, una procedura
per la determinazione di IC bootstrap–t con stabilizzazione della varianza incognita.
L’approccio non è invariante per trasformazione, eccetto per la traslazione. Tale problema risulta, però, comune a tutti quei metodi di stima
intervallare che utilizzano quantità pivotali standardizzate o studentizzate.
I vantaggi sono comunque notevoli: gli IC standard, così come tutti
tradizionali quelli fondati su quantità pivotali studentizzate, sono simmetrici
intorno al valore stimato b. Essendo i percentili derivanti da procedura
bootstrap–t non necessariamente simmetrici, gli intervalli di con…denza ri‡ettono naturalmente tale eventuale asimmetria, particolarmente e¢ cace nel
caso di distribuzioni tendenzialmente lontane dalla Normale. Più speci…catamente, nel caso non parametrico l’approccio tradizionale si basa su
5
un’approssimazione normale della distribuzione; il bootstrap–t utilizza la generazione di data–set della distribuzione empirica Fn e il calcolo di una quantità Tb(1 ) . In questa situazione il metodo porta ad una soluzione migliore
dell’approssimazione normale, in particolare in presenza di una numerosità
campionaria ridotta, e per la determinazione di intervalli di statistiche di
locazione (media, mediana, media troncata e percentili campionari).
Un ultimo, ma non per importanza, indubbio vantaggio riguarda la proprietà di accuratezza di cui godono gli intervalli bootstrap–t.
Si supponga di avere un parametro di interesse per il quale si costruisce
un intervallo di con…denza. Per comodità si consideri un IC ad una coda
e quindi un singolo valore limite b [ ]. Un estremo dell’intervallo si dice
accurato del primo ordine se
h
i
b [ ] = + O n 1=2
Pr
e del secondo ordine se
Pr
h
i
b[ ] =
1
+O n
:
Mentre gli IC standard godono solo dell’accuratezza del primo ordine se la
distribuzione non è normale, gli intervalli bootstrap–t, come si vedrà per tutte
le procedure di determinazione di intervalli e regioni con tecniche di ricampionamento, sono accurati almeno del secondo ordine.
La nozione di correttezza si riferisce, invece, a quanto un “possibile”limite
di un intervallo si avvicina ad un “ideale” (o, qualora esista, ad un esatto)
punto estremante di con…denza.
Sia bESAT T O [ ] un limite di con…denza esatto tale che
h
i
bESAT T O [ ] = .
Pr
Un intervallo di con…denza viene detto corretto del primo ordine se
b [ ] = bESAT T O [ ] + o n 1=2 b
e del secondo ordine se
b [ ] = bESAT T O [ ] + o n
1
b:
Si può agevolmente veri…care che la correttezza di un dato ordine implica
l’accuratezza dello stesso ordine.
Gli intervalli bootstrap–t (si veda per le dimostrazioni analitiche nel dettaglio Hall, 1988 e Di Ciccio & Efron, 1992) sono corretti, e quindi accurati,
del secondo ordine.
6
3
Intervalli di con…denza bootstrap percentile
In una serie di lavori e articoli, Efron [1981, 1985, 1987, 1996] ha introdotto
e de…nito una procedura, denominata Bootstrap percentile, per determinare
intervalli di con…denza per un parametro scalare (estendibile anche nel caso
multiparametrico) utilizzando approssimazioni bootstrap. Verrà ora esaminata la versione più semplice di tale metodo, denominata Intervallo di con…denza Bootstrap percentile. Nei successivi paragra… verranno illustrati gli
ulteriori sviluppi e perfezionamenti metodologici e computazionali.
L’approccio caratteristico del metodo percentile viene naturalmente sviluppato nel contesto più semplice di un modello parametrico dipendente da
un parametro scalare e di fatto adattato con buoni risultati per applicazioni
a famiglie multiparametriche e situazioni tipicamente non parametriche.
Sia X = (x1 ; x2 ; :::; xn ) un campione da una distribuzione avente funzione
di densità f dipendente da un parametro scalare . Sia b uno stimatore di
ottenuto da X con funzione di ripartizione G (k) = Pr b k . Il limite
superiore di con…denza di livello 1 –
che soddisfa la condizione
G
(1
di
è determinato dal valore
(1
)
b =
)
con G distribuzione bootstrap per b.
Si supponga ora che esista una trasformazione monotona crescente h e
una costante tale che, 8 ;
n
o
h b
h( )
Z
(7)
dove Z è simmetricamente distribuita attorno alla media nulla con funzione
di ripartizione H.
Allora
G (s) = H f [h (s) h ( )]g ;
di conseguenza, essendo z (
G
1
(1
G
1
)
=H
1
( )
)=h
e
( )=h
1
1
z (1
h b +
z(
h b +
7
)
)
:
Dal momento che (1 ) = G 1 (1
), l’IC di può essere ottenuto direttamente dalla funzione di ripartizione bootstrap di b senza esplicita conoscenza
della trasformazione h, supponendone esclusivamente l’esistenza. La quantità Gb 1 (1
) è detta il limite di con…denza bootstrap percentile.
Tipicamente non esiste alcuna trasformazione h tale per cui si ottenga
esattamente (7), e la di¤erenza tra Gb 1 (1
) e (1
) è nell’ordine di
1
o (n ); perciò il metodo percentile, senza ulteriori accorgimenti e modi…che, può risultare debole con dimensioni campionarie ridotte. In sostanza
il metodo è facilmente riassumibile nei seguenti passi.
Sia X un campione bootstrap generato da Fb ! X , dove Fb è la stima
di F , funzione di ripartizione della variabile X di campionamento, e sia
b = s (X ) una replicazione bootstrap della stima del parametro scalare
b
di interesse, con funzione di ripartizione pari a G.
L’intervallo percentile di livello nominale 1 2 è de…nito dai percentili
b
di livello e 1 – di G.
h
i
i h
b% inf ; b% sup = G
b 1( );G
b 1 (1
) :
b
Dalla de…nizione G
( ) = b ( ) si giunge a
h
i h
b% inf ; b% sup = b ( ) ; b (1
1
i
) :
[Es: si e¤ettuano 1000 replicazioni della stima di un parametro b , si …ssa
= 0:05; l’intervallo è semplicemente de…nito ordinando le 1000 stime e
considerando il 50-mo e il 950-mo valore delle replicazioni].
Se la distribuzione bootstrap di b è tendenzialmente normale, allora la
determinazione degli IC standard e con il metodo bootstrap percentile è pressoché la medesima. Di¤erenze, anche importanti, si hanno per campioni non
numerosi con distribuzioni fortemente asimmetriche (e quindi lontane dalla
situazione di normalità). In tale situazione la discrepanza tra IC standard
e IC bootstrap percentile è notevole, ad ovvio vantaggio del metodo bootstrap che più naturalmente asseconda la non simmetria distributiva. Ciò non
consegue che tale metodo sia ottimale, in special modo se n non è “su¢ cientemente”grande, e necessita di un’importante correzione, che nella letteratura
ha portato alla determinazione di intervalli di con…denza bootstrap percentile
corretti. Data la natura approssimata degli IC bootstrap percentile, se ne
discuteranno le proprietà di accuratezza e correttezza solo per le più precise
versioni corrette.
8
4
Intervalli di con…denza bootstrap percentile
corretti e accelerati
In questa sezione del lavoro vengono presentati i metodi di determinazione
di IC bootstrap derivanti dal metodo percentile, apportati con una sorta di
“correzione”.
L’intervallo di con…denza bootstrap corretto (BC) deriva sostanzialmente
dal metodo percentile precedentemente illustrato, e da un modello che generalizza ed estende le proposte del semplice metodo bootstrap percentile.
L’approccio migliorativo considera la distorsione in (7).
Si assume l’esistenza di una trasformazione monotona crescente h, e due
costanti e z0 tali che, 8 ;
n
o
h b
h ( ) + z0 Z
dove Z è presentata nella sezione precedente. In tale contesto
G (s) = H f [h (s)
h ( )] + z0 g
con le notazioni già utilizzate, da cui segue che
G
1
(1
G
1
)=h
e
( )=h
1
1
z (1
h b +
)
+ z0
z ( ) + z0
h b +
:
Dal momento che, per costruzione, (1
) = G 1 H z (1 ) + 2z0 , quindi
l’intervallo per il parametro scalare incognito può essere ottenuto dalla
distribuzione bootstrap di b. La correzione della distorsione z0 può essere
ottenuta in via del tutto analoga da
h
i
1
b
z0 = H
G
:
In luogo della funzione H, Efron (1982) suggerisce la distribuzione normale
standardizzata e de…nisce così
b !G
1
z (1
9
)
+ 2z0
il limite di con…denza bootstrap percentile corretto (BC). E’ evidente che,
qualora z0 sia pari esattamente a zero, l’approccio BC coincide con il bootstrap percentile.
Come per il caso del bootstrap percentile, tipicamente anche i limiti
del BC di¤eriscono dai limiti esatti, qualora essi esistano, per un valore
dell’ordine di grandezza pari a o (n 1 ). Nel caso, però, di un campione ridotto
o della stima di parametri in particolari situazioni, (si veda Schenker, 1985),
il metodo BC produce approssimazioni non su¢ cientemente accurate. Nonostante sia possibile mostrare il miglioramento rispetto al metodo percentile
più semplice, tale approccio risulta in parte non adeguato. Schenker (1985)
mostra infatti che la probabilità di copertura dell’IC BC è sostanzialmente
inferiore al suo livello nominale per piccoli campioni.
Per migliorare il bootstrap corretto e ovviare i problemi emersi, Efron e
altri hanno sviluppato un’ulteriore correzione, aggiungendo una costante a
detta accelerazione.
Utilizzando le medesime notazioni, si supponga l’esistenza di una trasformazione monotona crescente h, e di tre costanti, e z0 e a tali che, 8
8
9
<h b
h ( )=
+ z0 Z
(8)
: 1 + a h( ) ;
dove Z è precedentemente descritta.
Con questa assunzione
h (s) h ( )
+ z0 ;
1 + a h( )
G (s) = H
cui segue che
G
1
(1
)=h
e
G
1
( )=h
I limiti di con…denza di
momento che
(1
)
1
1
"
"
z (1 ) + z0 [1 + a h ( )]
h b +
f1 a (z (1 ) + z0 )g
h b +
z (1
)
+ z0 [1 + a h ( )]
#
#
:
si ottengono dalla distribuzione bootstrap di b, dal
=G
1
"
H
(
z0 +
10
1
z (1 ) + z0
a (z (1 ) + z0 )
)#
h
i
e z0 = H 1 G b . Sostituendo a H la distribuzione normale standard,
si ottengono i limiti di con…denza bootstrap percentile corretti accelerati BCa
(9)
" (
)#
(1
)
z
+
z
0
(1
)
=G 1
z0 +
:
(9)
(1
)
1 a (z
+ z0 )
Se a = 0 il metodo BCa coincide con il metodo BC.
Sia un parametro di interesse; siano b una stima di basata sui dati
osservati X , generalmente stima di massima verosimiglianza, e b = b (X ) è
una replicazione bootstrap di b ottenuta ricampionando X con procedimento
bootstrap da una stima della distribuzione originaria X .
b (s) la funzione di ripartizione delle B replicazioni bootstrap delle
Sia G
stime b (b) con b = 1; :::; B,
n
o
# b (b) < s
b (s) =
:
(10)
G
B
Il limite superiore di in un intervallo ad una coda BCa di livello è de…nito
b e dalle costanti z0 (correzione di distorsione) e a (accelerazione). La
da G
de…nizione del limite BCa è quindi
" (
)#
( )
z
+
z
0
1
z0 +
BCa ( ) = G
1 a (z ( ) + z0 )
funzione di ripartizione della v.c. normale standardizzata, con z ( ) =
1
( ).
Nel caso classico della ricerca di un intervallo di con…denza a due code ad
un livello di con…denza nominale di 1 2 , l’intervallo BCa è ottenuto da
h
i h
i
bBCa inf ; bBCa sup = b
b
(
)
;
(1
)
BCa
BCa
con
dove
BCa
( )=G
1
e
BCa
(1
)=G
1
" (
z0 +
" (
z0 +
11
)#
1
z ( ) + z0
a (z ( ) + z0 )
1
z (1 ) + z0
a (z (1 ) + z0 )
)#
:
Il problema diventa pertanto ora la stima delle due costanti a e z0 caratterizzanti il modello. Si noti, in prima istanza, che se entrambe sono nulle,
allora
1
( )
BCa ( ) = G
e l’intervallo BCa coincide con l’intervallo percentile. Nel caso ulteriore in
cui G è perfettamente Normale, allora l’IC BCa coincide con l’intervallo
standard.
In generale, evidentemente, a e z0 di¤eriscono da zero (anche se tipicamente hanno un campo di variazione ridotto), producendo di fatto IC anche fortemente asimmetrici intorno al valore stimato b. Gli intervalli BCa
b e le due costanti a e z0 .
richiedono di calcolare la distribuzione bootstrap G
La distribuzione bootstrap è ottenuta direttamente da (9). Tale calcolo, come
b (s), non richiede esplicita conoscenza della funsi ricava dalla de…nizione di G
zione monotona crescente h (8), ma esclusivamente di supporne l’esistenza.
Le due costanti possono essere ottenute direttamente dalla funzione di
distribuzione bootstrap delle stime b : La correzione di distorsione z0 è facilmente calcolabile e interpretabile: alla luce della costruzione attraverso l’impiego della funzione monotona crescente h, si de…nisce la costante z0 come
h
i
1 b b
z0 =
G
(11)
b b , funzione di ripartizione delle
e può essere calcolata direttamente da G
replicazioni bootstrap della stima del parametro. Quindi, nella sua formulazione più semplice, risulta
o3
2 n
# b (b) < b
14
5;
zb0 =
(12)
B
1
funzione inversa della normale standard
della proporzione di replicazioni
bootstrap strettamente minori del valore stimato b.
La seconda costante, detta accelerazione a, misura la velocità di propagazione dell’errore standard sulla scala normalizzata dei valori. La stima di a,
seppur più complessa di quella di z0 , è facilitata da una semplice relazione
che coinvolge la score function di Fisher.
Si supponga di avere una famiglia monoparametrica con funzione di ripartizione G b , con b stimatore almeno asintoticamente di massima verosimiglianza. La costante a, con buona approssimazione, è de…nita in termini di
12
asimmetria della score function:
dove
1
a = SKEW =b l b
6
n
@
b
=
log g
l
@
b
(13)
o
@G (b)
con g b densità di frequenza de…nita @b .
In Efron (1987), Section 10, si trova la dimostrazione completa del risultato (13). Viene di seguito riportata solo una traccia della stessa dimostrazione,
seguendo l’approccio di Shao et al.
Shao, discutendo il controesempio di Schenker riguardo all’ine¢ cacia del
metodo Bootstrap percentile corretto nella sua versione più semplice, introduce un’assunzione più generale che coinvolge l’asimmetria. Con le notazione
già utilizzate nel lavoro, si de…nisce Z
2
3
h b
h( )
Z = Pr 4
+ z0 x5
(14)
1 + a h( )
Si consideri il caso parametrico in cui F = Fp ; 2 <. Si supponga che
2
n(b )
var bn = n e i primi tre cumulanti 3 di tn =
soddis…no le seguenti
condizioni, per qualche funzione j :
( )
k1 (tn ) = p1
+O n
n
k2 (tn ) = 1 + O n
1
1
( )
k3 (tn ) = p3 3 + O n
n
1
Si supponga che esistano continue le derivate almeno …no al secondo ordine
della funzione monotona crescente h. E’possibile (Pfanzagl 1985), ottenere
3
1 , 2 , 3 , indicano i coe¢ cienti rispettivamente del I, II e III termine nell’espansione
di Edgeworth del logaritmo della funzione caratteristica di tn .
13
dalla (14) l’espansione
Pr
p
n[h(b) h( )]
h0 ( )
(x)
(x2 p1)'(x)
Se h soddisfa (14), allora
n
h
p1
n
3(
6
h
)
3
1(
+
)
+
h00 ( )
2h0 ( )
h00 ( )
2h0 ( )
i
i
x =
(15)
1
+ O (n )
h0 ( )
p
= 1 + a h( );
n
poiché h(x) genericamente può essere espresso come
82 x
9
3
Z p
<
=
1
a n 5
h (x) = exp 4
d
1 :
:
;
a
0
E’possibile determinare a tale che
2
h b
h( )
Pr 4
+ z0
1 + a h( )
3
x5 =
(x) + O n
1
(16)
Ora, ponendo il secondo termine di (15) uguale a zero, si ottiene, dopo alcuni
passaggi algebrici,
1 d
3( )
a = a( ) = p
3 3
n d
e
1
h00 ( )
1
1( )
1( )
3( )
z0 = z0 ( ) = p
+
=p
:
0
2h ( )
6 3
n
n
Se b è stimatore MLE di
a=
1
@3l ( )
E
6
@ 3
e l ( ) la funzione di log-verosimiglianza, allora
@2l ( )
E
@ 2
3=2
1
@l ( )
= SKEW
6
@
(17)
.
E’possibile altresì dimostrare che gli intervalli BCa godono delle proprietà
di accuratezza e correttezza del secondo ordine.
14
La dimostrazione, molto complessa e tecnica (si veda per i dettagli Efron,
1996), si basa sostanzialmente sulla prova che la distribuzione cumulata bootb b
b b
strap di T = b di¤erisce da T = b per una quantità pari a o n 3=2 ;
analogamente la correttezza si dimostra con medesimo procedimento, confrontando i limiti di con…denza bBCa ( ) con i limiti di con…denza esatti.
La complessità dei passaggi matematici (sostanzialmente connessi allo
studio dei primi tre cumulanti di ogni distribuzione) non suggerisce in questo
lavoro un’analisi più approfondita.
Uno sguardo d’insieme dei metodi …no ad ora proposti mostra che, mentre in particolari situazioni (ad esempio, determinazione di IC ad una coda)
il bootstrap percentile e il bootstrap percentile corretto sono accurati solo del
primo ordine, il metodo bootstrap–t e bootstrap corretto accelerato sono sempre accurati …no al secondo ordine. A fronte però di una maggiore accuratezza, l’utilizzo degli ultimi due approcci è talora complesso dal punto di
vista computazionale (il bootstrap–t richiede necessariamente un valido stimatore della varianza e una trasformazione per stabilizzare la stessa, il BCa
uno stimatore del parametro d’accelerazione), mentre gli altri appaiono di
utilizzo più immediato. Lo sviluppo, però, di linguaggi di programmazione
speci…ci (in particolare il software free R) ha reso però possibile semplicemente l’utilizzo di tutti i metodi qui descritti.
5
Intervalli di con…denza ABC
In…ne viene presentata un’approssimazione del metodo BCa, detto ABC (Approximate Bootstrap Con…dence). In molte situazioni, infatti, è possibile approssimare il BCa analiticamente, senza ricorrere all’utilizzo di replicazioni
Monte Carlo, portando, ovviamente, ad un dispendio computazionale fortemente ridotto.
Il metodo ABC trova naturale impiego come versione analitica del BCa
per determinati parametri di distribuzioni appartenenti alla famiglia esponenziale, anche se l’utilizzo maggiore si ha per particolari problemi di natura
non parametrica.
L’approccio ABC richiede un ulteriore parametro da stimare rispetto ai
due tipici del BCa (zb0 e b
a), detto parametro di non linearità cq , ma non
b
richiede la stima di G, funzione di ripartizione bootstrap del parametro b.
Generalmente la funzione di densità di una famiglia esponenziale k-parame15
trica può essere scritta come
g
(X )
= exp [ 0 y
( )] ;
dove X sono i dati campionari, y = Y (X ) è un vettore k-dimensionale di
statistiche su¢ cienti, è un vettore k-dimensionale di parametri naturali,
= M (y) e ( ) è una funzione dipendente solamente dai parametri. I
vettori e sono in corrispondenza biunivoca uno-a-uno; l’algoritmo ABC
sfrutta proprio la relazione tra e
=f( )
che fortunatamente ha semplice formulazione per tutte le più comuni distribuzioni appartenenti alla famiglia esponenziale.
Lo stimatore di massima verosimiglianza (MLE) di è b = y, cosicché lo
stimatore MLE di un parametro reale = t ( ) è
b = t (b) = t (y) :
Sia ( ) = cov (y) la matrice di varianze-covarianze k k di y e sia b =
(b) la stima di massima verosimiglianza di . La stima dell’errore standard
di b = t (b) dipende da b . Sia t il vettore gradiente di = t ( ) nel punto
=b
0
@t
; :::
:
t = :::;
@ i
=b
Allora
b = t 0b t
1=2
(18)
è la stima parametrica dell’errore standard.
I valori stimati di b utilizzati negli IC standard sono trovati per di¤erenziazione numerica
@t
@ i
=
t (b + "ei )
t (b
2"
b
"ei )
(19)
per valori di " piccoli a piacere.
In genere b è semplice da calcolare per la maggior parte delle distribuzioni
più comuni della famiglia esponenziale.
16
Il primo passo dell’algoritmo ABC avviene stimando b da (18) e (19).
Quindi i parametri a, z0 e cq sono stimati con k + 2 derivate seconde.
Le prime k derivate determinano b
a:
h
i
@2
0
f
(b
+
"t)
@"2 t
"=0
b
a=
;
(20)
6b3
con b stimatore MLE del vettore di parametri naturali (nel caso di famiglia
esponenziale monoparametrica (20) è identica alla stima di b
a nell’approccio
BCa (13)).
Quindi
b
@2
t b + "b t
@"2
"=0
cbq =
2b
e cq misura la “non linearità”del parametro in funzione di .
Le k derivate …nali sono richieste per la determinazione del parametro
di correzione di distorsione z0 . La stima dell’errore di b = t (b) può essere
espressa come
k
X
@2
1=2
bb = 1
t b + "di i
2 i=1 @"2
"=0
dove di è l’i-mo autovalore e
quindi
zb0 =
1
i
l’i-mo autovettore di b . La stima di zb0 risulta
2 (b
a)
cbq
bb
b
!!
=b
a + cbq
bb
.
b
La formulazione più semplice degli intervalli ABC, denominata ABC quadratico (ABCq ), fornisce limiti di con…denza di livello come funzione di cinque
valori:
! ! = zb0 + z ( )
= (1 !ba!)2 ! = + cq 2
! bABC ( ) = b + b
o espresso nella forma esplicita
8
"
#2 9
<
=
( )
( )
zb0 + z
zb0 + z
b
b+b
(
)
=
+
c
:
q
2
ABC
: 1 a (zb0 + z ( ) )2
1 a (zb0 + z ( ) ) ;
17
L’approccio ABC, in tutte le sue di¤erenti espressioni, è completamente automatizzato in algoritmi elettronici che ne rendono estremamente semplice
la determinazione computazionale. Come il metodo BCa, dal quale in parte
deriva, gode di ottime proprietà di accuratezza e correttezza …no al secondo
ordine.
6
Intervalli di con…denza Bootstrap – Caso
non parametrico
L’estensione dei metodi di costruzione degli intervalli di con…denza bootstrap
è facilmente adattabile in ambito non parametrico, dove le tecniche di ricampionamento trovano la loro naturale applicazione. In particolare, nel seguito,
si vedrà l’adattamento degli approcci bootstrap–t, BCa e ABC in tale nuovo
contesto.
Sia X = (x1 ; x2 ; :::; xn ) un campione osservato da un’incognita distribuzione di probabilità F
x1 ; x2 ; :::; xn
F:
i:i:d:
Sia Fb la distribuzione empirica di X , che assegna probabilità pari a n1 ad
ogni determinazione xi di X . Sia, quindi, = t (F ) un parametro scalare
di interesse e b = t Fb la sua stima non parametrica (stima di massima
verosimiglianza non parametrica).
Sia X = (x1 ; x2 ; :::; xn ) un campione bootstrap non parametrico casuale
di dimensione n estratto da Fb
x1 ; x2 ; :::; xn
Fb:
X è un campione casuale di dimensione n estratto con riposizione da X =
(x1 ; x2 ; :::; xn ). Il campione bootstrap X fornisce una replicazione bootstrap
di b,
b = t Fb ;
dove Fb è la distribuzione empirica bootstrap di X pari a n1 ad ogni deterb b (s) è la probabilità che
minazione xi . La funzione cumulativa bootstrap G
F
una replicazione bootstrap sia minore di s
b b (s) = Pr b < s :
G
F
Fb
18
b b è calcolabile tipicamente con simulazione Monte Carlo; in genere è possiG
F
b b è ricostruibile mostrare che, con almeno B=2000 replicazioni bootstrap, G
F
ta ottimamente.
I limiti di con…denza nel metodo percentile e BC sono facilmente determinabili: i limiti superiori di con…denza di livello 1 – per ottenuti dai
metodi percentile e BC sono, rispettivamente,
e
b 1
G
Fb
dove
z0 =
b 1 (1
G
Fb
z (1
n
o
1
bb b
G
=
F
)
)
1
+ 2z0
i9
8 h
< # b (b) < b =
:
B
;
:
(21)
L’implementazione del metodo BCa è meno diretta dal momento che richiede
la stima non parametrica dei due parametri a e z0 .
La formula (21) fornisce la stima di zb0 , mentre la stima dell’accelerazione
b
a è ottenuta utilizzando la “funzione empirica”Ui della statistica b = t Fb
t (1
Ui = lim
"!0
") Fb + "
i
"
t Fb
, i = 1; 2; :::; n
(22)
con i che indica il punto di massa su xi ; (1 ") Fb + " i è la versione di Fb
ponderando maggiormente i punti xi e meno gli altri punti. Ui risulta quindi
essere la derivata della stima b rispetto alla massa su xi . (22) presuppone
che t Fb sia di¤erenziabile in un intorno di Fb.
La formula dell’accelerazione a è
2
3
n
P
3
Ui
7
16
7
6 i=1
a= 6
(23)
7
3=2
n
64 P
5
2
Ui
i=1
relativamente semplice da utilizzare dal momento che Ui può essere valutata
facilmente usando le di¤erenze …nite in (22).
(23) appare completamente diversa rispetto a (17); di fatto si tratta della
stessa formulazione, in ambito non parametrico, di un’appropriata situazione
multinomiale (si veda Efron, 1987, Section 8).
19
6.1
Esempio: La Media
Sia F una distribuzione de…nita nel supporto reale, = t (F ) è il valore
atteso EF (X). La funzione Ui è de…nita come (xi x) e (23) diventa
8
9
8
9
n
n
P
P
>
>
>
>
>
>
>
>
>
>
(xi x)3 >
(xi x)3 n >
<
=
<
=
i=1
i=1
1
1
= 6p n
=
a= 6
3=2 >
3=2 >
n
n
>
>
P
P
>
>
>
>
>
>
>
>
:
;
:
;
(xi x)2
(xi x)2 n
i=1
i=1
!
1
b
b3
1
b3
= p
= p
= p
3=2
3=2
6 n b2
6 n b
6 n
dove bk è il k-momento centrale campionario e b l’indice di asimmetria (skewness) di Fisher campionario.
Nel caso più generale, il metodo più semplice per la stima di b
a è data
b
sfruttando le quantità jacknife della statistica di interesse = s (X ). Sia
X(i) il campione originario di dati privo del punto i-mo xi ; sia b(i) = s X (i)
n b
P
(i)
la statistica di interesse stimata nel campione x(i) . Sia in…ne b( ) =
.
n
i=1
La più semplice formulazione dell’accelerazione in ambito non parametrico è
data da
n
3
P
b( ) b(i)
b
a=
i=1
n
P
b( )
6
i=1
b(i)
2
3=2
:
Gli intervalli di con…denza BCa non parametrici sono corretti del secondo
ordine per ogni tipologia di parametro di interesse. L’estensione del criterio ABC nel caso non parametrico è decisamente semplice e vantaggiosa in
quanto riduce sensibilmente il “peso”computazionale dell’approccio BCa. La
procedura, che richiede in tale ambito di un numero ridotto di rivalutazioni,
segue strettamente quella parametrica già presentata, ed è quindi adattata
in un contesto non parametrico.
20
7
Confronto delle tecniche per la determinazione di un intervallo di con…denza per
un indice di cograduazione
Sono stati, in…ne, confrontati i diversi approcci per la determinazione degli
intervalli di con…denza bootstrap per un indice di cograduazione per valutare
la propensione di una scala di misurazione (comunemente categoriale ordinale) ad assegnare giudizi maggiormente positivi ovvero negativi rispetto ad
un’altra caratterizzata da uno stesso numero di passi di risposta.
Un problema tipico nella statistica psicometrica è infatti la determinazione della scelta di scala di risposta ottimale. Parecchi studi (si veda
in particolare Bonanomi, 2004) hanno dimostrato che l’utilizzo di scale di
misurazione di¤erenti, pur con lo stesso numero di passi di risposta, porta
a “valutazioni” anche sensibilmente di¤erenti. In Bonanomi (2004) viene
proposto un indice di cograduazione o di propensione per misurare la tendenza ad attribuire giudizi maggiormente positivi ovvero negativi utilizzando
di¤erenti scale di risposta.
Siano S1 e S2 due scale di risposta ordinate caratterizzata da m passi
di risposta, rispettivamente s1i e s2j , i=1,. . . ,m, j=1,. . . ,m. Ad uno stesso
gruppo di n soggetti viene proposto uno stesso stimolo (o item) in un questionario di valutazione di un particolare bene/servizio. Lo scopo dell’indice
proposto è quello di valutare la propensione ad attribuire punteggi maggiormente positivi o negativi nell’utilizzo di una scala.
Determinata la tabella di contingenza (Tab 1) che raccoglie le risposte
degli soggetti all’item proposto valutato nelle due scale
Tab.1: Tabella di contingenza delle frequenze assolute delle risposte dei
soggetti con Scala 1 e Scala 2
Scala 2
Scala 1 s21
s1j
s1m Tot
s11
n 11
n 1j
n 1m n 1
..
..
..
..
..
..
..
.
.
.
.
.
.
.
s1i
..
.
n i1
..
.
s1m
Tot
n m1
n1
..
.
n ij
..
.
n mj
nj
21
..
.
n im
..
.
ni
..
.
n mm
nm
nm
n
si de…nisce l’indice di propensione rapportando la di¤erenza tra la somma
delle frequenze presenti nel triangolo superiore alla diagonale principale e la
somma delle frequenze collocate nel triangolo inferiore e per il totale delle
frequenze n.
In termini formali, si indichi tale indice di propensione con Ip
Ip =
I
J
P
P
nij
i=1 j=i+1
J
I
P
P
j=1 i=j+1
n
nij
1
Ip
1:
(24)
Un valore positivo dell’indice evidenzia la propensione sistematica ad attribuire valori più elevati quando si utilizza la scala i cui termini sono posti
per colonna della tabella (Scala 2); viceversa un valore negativo dell’indice
evidenzia la propensione sistematica ad attribuire valori più elevati quando si
utilizza la scala i cui termini sono posti per riga (Scala 1): Qualora tale indice
assume il valore esatto di 0, siamo in situazione di indi¤erenza, cioè non si
evidenziano propensioni né per una scala né per l’altra. Quando si veri…ca
perfetta concordanza tra le scale, cioè che tutte le frequenze sono concentrate sulla diagonale principale (corrispondenza perfetta tra i diversi ranghi
della scala), allora l’indice assume esattamente il valore 0. Se le discordanze
(cioè le frequenze tendono a concentrarsi sopra la diagonale principale ovvero
sotto) sono a favore della scala in colonna, allora signi…ca che è presente una
propensione all’attribuzione di ranghi maggiormente elevati per tale scala.
Il valore massimo che tale indice assume è pari a 1, quando tutte le frequenze sono concentrate nel triangolo superiore della tabella; viceversa il
valore minimo che tale indice assume è pari a 1, quando tutte le frequenze
sono concentrate nel triangolo inferiore della tabella. Tale indice risulta pertanto normalizzato (in senso improprio) tra gli estremi –1 e 1.
L’utilizzo di intervalli di con…denza bootstrap non parametrici va visto,
quindi, nell’ottica della costruzione di un test, evidentemente non parametrico, del tipo
H0 : Ip = 0
H1 : Ip 6= 0
per saggiare l’ipotesi nulla che l’indice di propensione sia uguale a zero, cioè
che sia indi¤erente valutare un item con una scala ovvero con l’altra.
A tal proposito è stato calcolato l’indice di propensione (24) per alcuni
item di un questionario sulla didattica universitaria. Sono stati confrontati gli
approcci di costruzione di intervalli di con…denza bootstrap con l’approccio
22
classico degli intervalli standard, al variare della numerosità campionaria,
per valutare l’accettazione di H0 , e valutati i miglioramenti in termini di
“ampiezza dell’intervallo” passando dall’approccio standard ai metodi non
parametrici di ricampionamento bootstrap. L’applicazione è realizzata confrontando due scale psicometriche di misurazione, una gra…ca a 5 intervalli
(Scala 1 –punteggio minimo 1, punteggio massimo 5) e una verbale a cinque
passi (Scala 2, esplicitata nelle modalità ordinate insu¢ ciente, mediocre, suf…ciente, buono, ottimo). E’stato chiesto di valutare l’aspetto didattico generale di un’università prima da 250 studenti, poi da solo 45 studenti della
facoltà di giurisprudenza, in…ne da 12 studenti di una facoltà economicobancaria.
Le tre successive tabelle di contingenza riportano le distribuzioni di frequenze congiunte rispettivamente nel caso di n=250 (Tab.2), n=45 (Tab.3),
n=12 (Tab.4).
Scala 1 insu¤
1
2
2
2
3
4
5
Tot
4
Tab. 2 : n=250, Ip=0,176
Scala 2
mediocre
su¤
buono
13
3
14
18
3
10
61
38
2
9
59
1
8
39
92
108
Scala 1 insu¤
1
3
2
3
4
5
Tot
3
Tab. 3 : n=45, Ip=0,177
Scala 2
mediocre su¤
buono
2
1
6
5
2
11
3
3
4
1
10
20
8
23
ottimo
1
6
7
ottimo
1
2
1
4
Tot
18
37
109
71
15
250
Tot
6
11
17
9
2
45
Scala 1 insu¤
1
2
3
4
5
Tot
0
Tab. 4 : n=12, Ip=0,333
Scala 2
mediocre su¤
buono
1
1
1
2
2
1
2
2
4
4
ottimo
1
1
2
Tot
1
2
4
4
1
12
Seguono poi i risultati della costruzione degli intervalli di con…denza standard e bootstrap relativi all’indice di propensione Ip ; le applicazioni sono
state sviluppate in ambiente Microsoft Excel 2000 e R 2.3.1, (free version).
Tutte le simulazioni bootstrap sono state realizzate con un numero di replicazioni pari a B=2000.
Nel caso di numerosità campionaria ridotta si può notare un comportamento eterogeneo dei metodi bootstrap rispetto all’approccio standard (che
non ha, in tale situazione, ragione d’esistenza). Al crescere della numerosità
campionaria i metodi bootstrap appaiono sistematicamente migliori in termini di ampiezza dell’intervallo), al di là delle già enunciate maggiori proprietà di accuratezza e correttezza. In particolare il Bootstrap –t con varianza
stabilizzata sembra essere il migliore. Nelle successive tabelle i risultati delle
elaborazioni nelle tre situazioni.
Tab. 5: Confronto degli IC, n=250, Ip=0,177
Lim sup Lim inf Ampiezza A da Standard
Standard
0,097
0,255
0,158
-
Boot-t
0,099
0,255
0,157
-0,69%
Boot-t varianza stabilizzata
0,107
0,247
0,140
-13,18%
Boot percentile
0,095
0,252
0,157
-0,69%
Bca
0,096
0,248
0,152
-4,00%
ABC
0,098
0,254
0,156
-1,07%
24
Tab. 6: Confronto degli IC, n=45, Ip=0,177
Lim sup Lim inf Ampiezza A da Standard
Standard
-0,017
0,373
0,390
-
Boot-t
-0,012
0,366
0,378
-3,22%
Boot-t varianza stabilizzata
-0,022
0,356
0,378
-3.37%
Boot percentile
-0,022
0,356
0,378
-3,35%
Bca
-0,022
0,333
0,356
-9,81%
ABC
-0,014
0,366
0,379
-2,95%
Tab. 7: Confronto degli IC, n=12, Ip=0,333
Lim sup Lim inf Ampiezza A da Standard
Standard
-0,035
0,702
0,737
-
Boot-t
-0,108
0,693
0,801
+7,98%
Boot-t varianza stabilizzata
0,031
0,623
0,592
-24,47%
Boot percentile
0,000
0,667
0,667
-10,58%
Bca
-0,167
0,583
0,750
+1,73%
ABC
-0,050
0,659
0,709
-3,90%
Bibliogra…a
1. Bonanomi A. (2004), Variabili ordinali e trasformazioni di scala, con
particolare riferimento alla stima dei parametri dei modelli interpretativi con variabili latenti, Tesi di dottorato di Ricerca in Statistica
Metodologica ed Applicata, sede amministrativa Università degli Studi
di Milano - Bicocca
2. Diciccio T.J., Efron B. (1996), Bootstrap Con…dence Intervals, Statistical Science, Vol. 11, No. 3, pp 189-212
25
3. Diciccio T.J., Efron B. (1992), More Accurate Con…dence Intervals in
Exponential Families, Biometrika, Vol. 79, No. 2, pp 231-245
4. DiCiccio, T.J., Romano J.P. (1988), A Review of Bootstrap Con…dence
Intervals, Journal of the Royal Statistical Society, Series B, Vol. 50,
No. 3., pp 338-354
5. Efron B. (1981), Nonparametric Standard Errors and Con…dence Intervals, The Canadian Journal of Statistics / La Revue Canadienne de
Statistique, Vol. 9, No. 2, pp. 139-158
6. Efron B. (1985), Bootstrap Con…dence Intervals for a Class of Parametric Problems, Biometrika, Vol. 72, No. 1 (Apr., 1985), pp. 45-58
7. Efron B. (1987), Better Bootstrap Con…dence Intervals, Journal of the
American Statistical Association, Vol. 82, No. 397., pp 171-185
8. Efron B, Tibshirani R. (1993), An Introduction to the Bootstrap, Chapman and Hall, New York.
9. Hall P. (1986), On the Bootstrap and Con…dence Intervals, The Annals
of Statistics, Vol. 14, No. 4, pp 1431-1452
10. Pfanzagl J. (1985), Asymptotic Expansions for General Statistical Models, Springer-Verlag, Berlin.
11. Quatto P., Salini S., Zambon A. (2002), Intervalli di Con…denza Bootstrap nel campionamento per cattura e ricattura, Studi in Onore di
Angelo Zanella, Vita e Pensiero, Milano
12. Schenker N. (1985), Qualms about Bootstrap Con…dence Intervals,
Journal of the American Statistical Association, Vol. 80, pp360-361
13. Shao J., Tu D. (1995), The Jacknife and Bootstrap, Springer, New York
26