Cap. 19 Test statistici Cioè come verificare ipotesi statistiche utilizzando i dati campionari 1.
Test di ipotesi - old · 16/01/2014 1 Test di ipotesi In una carta - diremo che un processo è in...
Transcript of Test di ipotesi - old · 16/01/2014 1 Test di ipotesi In una carta - diremo che un processo è in...
16/01/2014
1
Test di ipotesi
In una carta - diremo che un processo è in controllo sta-
tistico se per ogni , indice dei sottogruppi,
( , ).k
X bar
k
x LCL UCL∈
Prendere questa decisione equivale ad effettuare un test di ipotesi.
Per effettuare un test di ipotesi è necessario:
1 2 selezionare un campione casuale , , nx x x• …
calcolare il valore di una statistica sul campione (nel CSQ è )k
x•
caratterizzare un intervallo (nel CSQ è ( , )) che con-
senta di rigettare l'ipotesi iniziale quando la statistica test non
vi appartiene (nel CSQ equivale a dire che non
LCL UCL
il processo è
i
•
)n controllo statistico
16/01/2014
2
In generale, si tratta di decidere se il parametro di una certa popolazio-
ne assume un determinato valore.
0se la media della popolazione assume ilAd esempio valo, re .µ µ
Nel caso del CSQ, si tratta di decidere se la media della popolazione as-
sume valore pari al target.
( )
( )0 0
1 0
: ipotesi nulla
: ipotesi alternativa
H
H
µ µ
µ µ
=
≠
Nell’esempio del contenuto dei flaconi, potremmo essere interessati
a decidere se la media della popolazione coincide con il valore target 95.
La media campionaria risulta essere 98.51. Il valore 98.51 è sufficiente-
mente vicino al valore 95 per ritenere che la media della popolazione
coincide con 95?
Per rispondere a questo quesito, abbiamo bisogno di fissare una distan-
za del valore osservato per la media campionaria dal valore target 95:
se 98.51 dista da 95 meno di questa distanza prefissata, allora riterremo
95 un valore credibile per la media della popolazione.
In sostanza si tratta di fissare un intervallo, che contenga 95 (in genere
si assume sia centrato su 95) tale che se, 98.51 appartiene a tale inter-
vallo, si ritiene plausibile, l’ipotesi nulla:
( )0 : 95 ipotesi nullaH µ =
95
98.51? REGIONE DI ACCETTAZIONE
16/01/2014
3
( )1 : 95 ipotesi alternativaH µ ≠
9598.51?REGIONE CRITICA
Per evitare che questa scelta sia soggettiva, possiamo fissare a-priori la
probabilità di commettere un errore prendendo una decisione sulla
media della popolazione.
Che tipi di errori si possono commettere?
( )0 0si rigetta H - a posteriori - quando H è vera - a prioriPα =
livello di significatività del testα ←ERRORE DI I TIPO
ERRORE DI I TIPO
( )0 01 si accetta H - a posteriori - quando H è vera - a prioriPα− =
La probabilità dell’evento complementare si riferisce a una decisione
corretta:
16/01/2014
4
( )1 1si rigetta H - a posteriori - quando H è vera - a prioriPβ =
1 potenza del testβ− ←
( )1 11 si accetta H - a posteriori - quando H è vera - a prioriPβ− =
ERRORE DI II TIPO
( )0 01 si rigetta H - a posteriori - quando H è falsa - a prioriPβ− =7
ERRORE DI II TIPO
La probabilità dell’evento complementare si riferisce a una decisione
corretta:
( )0 01 si accetta H - a posteriori - quando H è vera - a prioriPα− =
( )0 0 0 0si accetta quando = ,H xµ µ µ δ µ δ⇒ ∈ − +
Statistica osservata0µ
( )0 0 0 0uando = , ( ,...) ,Q X Nµ µ µ µ δ µ δ≈ ⇒ − +
1 α− ⇒
δ−δ
?n
σδ =
/2zn
α
σδ =
16/01/2014
5
livello di significatività del testα ←
( ) regione di accettazione 1P STATISTICA α∈ = −
0.10,0.05,0.01α =
REGIONE DI
ACCETTAZIONE
? ?
9
Quale si sceglie di fissare?
/21 P X zn
α
σα µ
− = ∈ ±
Come?
/21 P X zn
α
σα µ
− = ∈ ±
?
/2 01 P X zn
α
σα µ µ µ
− = ∈ ± =
0 /2 01 P X z
nα
σα µ µ µ
− = ∈ ± =
Regione di accettazione
Si rigetta l’ipotesi nulla…
Statistica osservata0µ
10
Pertanto fissare l’errore di I
tipo equivale a fissare la
regione di accettazione, e
viceversa.
16/01/2014
6
L’errore di II tipo
Statistica osservata
1µ µ=0µ µ=
0µ 1µ
11
La potenza del test
0µ µ=1µ µ=
0µ 1µStatistica osservata
12
E’ un valore che può
essere calcolato solo
se si conosce un valore
alternativo. Quale?
16/01/2014
7
Le statistiche test vanno scelte a secondo delle informazioni che si
hanno sul campione.
Ad esempio,
1) il campione proviene da una popolazione gaussiana (qqnorm )
2) la varianza è nota ,X N
n
σµ
≈ Ad esempio,
1) il campione proviene da una popolazione gaussiana (qqnorm )
2) la varianza non è nota 1n
XT
Sn
µ−
−≈
Ad esempio,
1) il campione non proviene da una popolazione gaussiana (qqnorm)
2) la varianza non è nota
3) la taglia supera 30 ,S
X Nn
µ
≈
Torniamo all’esempio dei flaconi, e verifichiamo se il campione proviene
da una popolazione con media 95.
Non conoscendo la varianza teorica, la statistica da utilizzare ha legge
T-Student.
La function è t.test().
> t.test(data,mu=95)
One Sample t-test
data: data
t = 4.5358, df = 119, p-value = 1.379e-05
alternative hypothesis: true mean is not equal to 95
95 percent confidence interval:
96.98146 100.05187
sample estimates:
mean of x
98.51667
A quale conclusione giungiamo?
16/01/2014
8
Statistica test=4.53
Legge T-Student
con 119 gradi di libertà
Regione
critica?
CHE COSA
E’ IL P-VALUE?
L’area a destra
della
statistica test.
0
Se
/ 2
si rigetta
p value
H
α− <
⇒
Statistica
test=4.53
Nell’esempio, il p-value
è molto inferiore a 0.025,
quindi l’ipotesi nulla
si rigetta.
0
0
Se stat reg.accett. non si rigetta
Se stat reg.accett. si rigetta
p H
p H
α
α
≥ ⇒ ∈ ⇒
< ⇒ ∉ ⇒
16/01/2014
9
Avendo stabilito che la media della popolazione non è 95, potremmo
essere interessati a sapere se la media è maggiore o minore di 95.
In quel caso si effettua un test a una coda.
REGIONE DI
ACCETTAZIONE
REGIONE
CRITICA
A 1 coda
0 0
1 0
:
:
H
H
µ µ
µ µ
=
>Si rigetta l’ipotesi nulla
perché la media
campionaria ha un valore
che ricade nella regione
critica. > t.test(data,alternative=c("greater"),mu=95)
One Sample t-test
data: data
t = 4.5358, df = 119, p-value = 6.895e-06
Si rigetta l’ipotesi nulla, in favore di
quella alternativa perché il p-value
é inferiore a 0.05.
Potremmo usare un valore del livello di significatività diverso da 0.05.
> t.test(data,alternative=c("greater"),mu=95,conf.level=0.99)
One Sample t-test
data: data
t = 4.5358, df = 119, p-value = 6.895e-06
alternative hypothesis: true mean is greater than 95
99 percent confidence interval:
96.68839 Inf
sample estimates: mean of x 98.51667
> t.test(data,alternative=c("greater"),mu=95)
One Sample t-test
data: data
t = 4.5358, df = 119, p-value = 6.895e-06
alternative hypothesis: true mean is greater than 95
95 percent confidence interval:
97.23138 Inf
sample estimates: mean of x 98.51667
Cambia solo
l’intervallo di
confidenza
16/01/2014
10
Per il caso
1) il campione non proviene da una popolazione gaussiana (qqnorm)
2) la varianza non è nota
3) la taglia supera 30 ,S
X Nn
µ
≈
si continua ad utilizzare il t.test() poiché la v.a. T-Student è approssimata da
una v.a. gaussiana
Per il caso
1) il campione proviene da una popolazione gaussiana (qqnorm )
2) la varianza è nota
,X Nn
σµ
≈
In questo caso la libreria da utilizzare è BSDA (basic statistical data analysis) che a sua
volta necessita di e1071, class, lattice. La function è z.test(). La sintassi è:
> z.test(difetti, y = NULL, alternative = "two.sided", mu = 20, sigma.x = 7.8, sigma.y = NULL,
+ conf.level = 0.99)
Il risultato per il dataset difetti.txt è
In tal caso il p-value è molto inferiore a 0.05.
Pertanto l’ipotesi nulla si rigetta.
> z.test(difetti, y = NULL, alternative = "two.sided", mu = 20, sigma.x = 7.8, sigma.y =
+ NULL,conf.level=0.99)
One-sample z-Test
data: difetti
z = -5.922, p-value = 3.181e-09
alternative hypothesis: true mean is not equal to 20
99 percent confidence interval:
7.898483 15.234850
sample estimates:
mean of x
11.56667
>
16/01/2014
11
Diremo che il processo è in controllo statistico se per ogni
, indice dei sottogruppi, ( , ).i
i x LimInf LimSup∈
Regione di accettazione
0 0
1 0
(rigettare | )
(rigettare | )
P H
P H
α µ µ
β µ µ
= =
= ≠
0
0
( ( , ) | )
( ( , ) | )
t
t
P x LCL UCL
P x LCL UCL
µ µ
µ µ
= ∉ =
= ∈ ≠
FALSO ALLARMEMANCATO ALLARME21
CURVA CARATTERISTICA OPERATIVA
Non avendo ipotesi alternative certe, immaginiamo che 1 0 kµ µ µ σ= = +
Se la popolazione è gaussiana, allora
0( ( , ) | )tP x LCL UCL kβ µ µ σ= ∈ = +
( ) ( )0 0
/ /
UCL k LCL k
n n
µ σ µ σ
σ σ
− + − + = Φ − Φ
Il plot dei valori assunti da questo parametro per un opportuno valore
di k, si chiama curva caratteristica operativa.
In R l’istruzione per effettuare questo grafico è oc.curves()
> obj<-qcc(data,type='xbar')
> oc.curves(obj)
16/01/2014
12
vari valori di k
Specificando un valore pari a TRUE per la variabile identify, è possibile
identificare interattivamente dei punti sul grafico
> beta<-oc.curves(obj)
Vari valori di beta possono essere ottenuti specificando un parametro
di output:
Se la media del processo è
95+0.4*st.dev, per n=5, la
prob. di un mancato allarme è 0.98
Per valori di k inferiori, aumenta la
probabilità di un mancato allarme.
16/01/2014
13
Nella progettazione delle carte di controllo è necessario specificare sia la dimensione
del campione che la frequenza di campionamento.
• Più grande è il campione più sensibile è il rilevamento di una variazione all’interno del
processo.
• La pratica corrente tende a diminuire la dimensione del campione e ad aumentare la
frequenza di campionamento.
Si fissa , e si cerca quel valore di tale che nβ
25
ALTRO USO DELLA CURVA OPERATIVA
( ) ( )0 0
/ /
UCL k LCL k
n n
µ σ µ σβ
σ σ
− + − + = Φ − Φ
0 0Se 3 e 3 , alloraUCL LCL
n n
σ σµ µ= + = −
( ) ( )3 3k n k nβ = Φ − − Φ − −
( ) ( )Si fissa , e si cerca quel valore di tale che
ossia, ricordando le proprietà della gaussiana...
z z zβ β ββ βΦ − Φ − =
( )2 1zβ βΦ − =1 1
2zβ
β− + ⇒ = Φ
3
è possibile ricavare
z k n
n
β⇒ = −
⇒
Ad esempio per 0.3β = > qnorm(1.3/2)
[1] 0.3853205
23 z
nk
β− ⇒ =
e per 1k = > (3-qnorm(1.3/2))^2
[1] 6.836549
n=6
16/01/2014
14
27
Cambiamento della macchina per imballaggio?
Sono statisticamente differenti le due produzioni?
IL CONFRONTO TRA DUE POPOLAZIONI
1 2 1 2( , , , ) ( , , , )n n
x x x y y y… …
POPOLAZIONE 1 POPOLAZIONE 2
• I due campioni provengono dalla stessa popolazione?
• Le due popolazioni hanno la stessa media?
• Le due popolazioni hanno la stessa varianza?
• Se provengono da popolazioni distinte, queste sono indipendenti
oppure c’è una qualche forma di relazione tra loro?
CONFRONTO TRA DUE POPOLAZIONI
16/01/2014
15
29
Le due popolazioni hanno la stessa varianza?
Supponiamo di aver monitorato la produzione di una sbarretta metallica
in una azienda in due periodi diversi dell’anno. Vogliamo verificare se la
produzione è rimasta costante nel tempo.
I dati sono in due files distinti: produzione1.txt e produzione2.txt
> prod1<-matrix(scan('C:/Programmi/R/R-3.0.2/produzione1.txt'),ncol=1, byrow=TRUE)
Read 20 items
> prod2<-matrix(scan('C:/Programmi/R/R-3.0.2/produzione2.txt'),ncol=1, byrow=TRUE)
Read 20 items
> var(prod1);var(prod2)[,1]
[1,] 1.574374
[,1]
[1,] 1.093336
>var.test(prod1, prod2, ratio = 1, alternative =
+ "two.sided", conf.level = 0.95)
F test to compare two variances
data: prod1 and prod2
F = 1.44, num df = 19, denom df = 19, p-value = 0.4341
alternative hypothesis: true ratio of variances is not equal
to 1
95 percent confidence interval: 0.5699588 3.6380212
sample estimates: ratio of variances 1.439973 Non si rigetta
( ) ( )( )
1 2
1 2
2 2
2 2
1 2 1 2
1 2
Siano e due var. al. con legge chi-quadrato con gradi
di libertà rispettivamente e . La var. al. / / /
ha legge di Fisher di gradi di libertà , .
n n
n nn n n n
n n
χ χ
χ χ
( )2
12 1 22
Fisher 1, 1S
n nS
≈ − −
221
1 12
1
222
2 22
2
( 1)
( 1)
Sn
Sn
χσ
χσ
− ≈
− ≈
2
1
12
2
2
1
1
n
n
χ
χ−
−
F-TEST DISTRIBUZIONE DISTRIBUZIONE DISTRIBUZIONE DISTRIBUZIONE DIDIDIDI FISHERFISHERFISHERFISHER
16/01/2014
16
Le due popolazioni hanno la stessa media?
La function da usare è sempre la stessa, t.test()
> t.test(prod1, prod2, mu = 0, alternative ="two.sided",paired = FALSE, var.equal =
+ TRUE, conf.level = 0.95)
Two Sample t-test
data: prod1 and prod2
t = 1.3622, df = 38, p-value = 0.1812
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval: -0.2418487 1.2368487
sample estimates: mean of x mean of y
100.2350 99.7375
Le procedure sono distinte a secondo che le varianze siano uguali oppure diverse.
PAIRED va posto uguale a TRUE quando i due campioni casuali si riferiscono alle stes-
se unità statistiche e quindi c’è una forma di dipendenza tra i due campioni. In tal caso
le taglie devono essere uguali.
Se PAIRED è posto uguale a FALSE, le taglie dei due campioni possono essere diverse.
Non si rigetta l’ipotesi nulla.
I due campioni provengono dalla stessa popolazione?
L’idea è quella di
costruire le funzioni
di ripartizioni
empiriche per i due
campioni e poi
di valutare la distanza
massima tra queste
ultime
Il test di Kolmogorov-Smirnov consente di stabilire se due campioni pro-
vengono da due popolazioni aventi la medesima legge di probabilità
In R la function da usare è ks.text(). Nella variabile di input alternative è possibile specificare
se l’ipotesi alternativa prevede che una funzione di ripartizione sia maggiore dell’altra o vice-
versa.
16/01/2014
17
> ks.test(prod1,prod2)
Two-sample Kolmogorov-Smirnov test
data: prod1 and prod2
D = 0.25, p-value = 0.5596
alternative hypothesis: two-sided
Warning message:
In ks.test(prod1, prod2) : cannot compute exact p-value with ties
> plot(ecdf(prod1),xlim=range(c(prod1,prod2)),
+ main='Confronto f.r.empiriche')
> par(new=TRUE)
> plot(ecdf(prod2),add=TRUE,col='red')
> legend(100,0.3,c("prod1","prod2"),
+ lty=c(1,1),lwd=c(2.5,2.5),col=c("black","red"))
>
• Se provengono da popolazioni distinte, queste sono indipendenti
oppure c’è una qualche forma di relazione tra loro?
Per rispondere a questo quesito, è necessario prima calcolare il coefficiente di corre-
lazione tra i due datasets:
> cor.test(prod1,prod2)
Pearson's product-moment correlation
data: prod1 and prod2
t = 0.9373, df = 18, p-value = 0.361
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2507354 0.6008810
sample estimates:
cor
0.2157209
Non si rigetta
l’ipotesi che li
coefficiente di
correlazione è
nullo.
Se le popolazioni da cui provengono i due campioni hanno legge gaus-
siana, allora risultano indipendenti. Poiché abbiamo mostrato che le due
popolazioni hanno medesima distribuzione, basta provare per un solo
campione che la popolazione da cui proviene segue una legge gaussiana.
16/01/2014
18
In questo caso, il valore del secondo parametro deve essere una stringa numerica con-
tenente il nome della funzione distribuzione del modello teorico che si vuole testare.
> ks.test(prod1,"pnorm",mean(prod1),sd(prod1))
One-sample Kolmogorov-Smirnov test
data: prod1
D = 0.1326, p-value = 0.8289
alternative hypothesis: two-sided
Siccome il p-value > 0.05
l’ipotesi nulla non si ri-
getta.
Trattandosi di campioni gaussiani, con un coefficiente di correlazione
quasi nullo, si può ritenere siano indipendenti.
Il ks-test può essere usato solo con modelli teorici continui. Nel caso
discreto, è possibile usare il test chi-quadrato.
Ad esempio, potremmo essere interessati a sapere se i dati del file
difetti.txt hanno effettivamente legge binomiale.
TEST CHI-QUADRATO
In questo test, vengono confrontati gli istogrammi costruiti sul campione con quelli che
si ottengono adoperando un modello teorico.
> setwd("C:/Programmi/R/R-3.0.2/")
> data<-read.table("difetti.txt",header=TRUE)
> attach(data)
> boundaries<-seq(0,max(difetti))
> tab.out<-table(cut(difetti,boundaries))
> appoggio<-cbind(tab.out)
>
La function cbind crea una matrice con il risultato della
function table. La function table crea una tavola di contingenza,
con il risultato della function cut, che calcola le frequenze
di occorrenza della modalità del campione. Il risultato è
Il campione va ripartito in modalità e per ciascuna di queste va contato il numero di occorrenze
nel campione.
> appoggio[1:24]
[1] 0 0 0 1 2 2 2 2 3 3 1 3 2 1 2 1 1 1 0 1 0 1 0 1
16/01/2014
19
Ad esempio, per il file difetti.txt, il modello teorico è quello di una binomiale di parametri n=50
e p pari alla linea centrale della carta di controllo p, ossia 0.23.
Il test chi-quadrato valuta le differenze tra queste due curve, normalizzate all’ordine di
grandezza delle frequenze teoriche.
La libreria di riferimento è la vcd.
Per un modello teorico discreto, l’analogo dell’istogramma è il diagramma delle
frequenze assolute.
> gf<-goodfit(difetti, type=“binomial", par = list(prob = 0.23, size = 50))
> summary(gf)
Goodness-of-fit test for binomial distribution
X^2 df P(> X^2)
Pearson NaN 50 NaN
Likelihood Ratio 34.96249 17 0.006292134
Warning message:
In summary.goodfit(gf) : Chi-squared approximation
may be incorrect
Per ottenere il grafico , usare le variabili
gf$observed; gf$count, gf$fitted che si ottengono
usando str(gf).
Il motivo per cui il test non restituisce un valore finito per il p-value è la presenza di troppi
0 nelle frequenze osservate. Vale infatti la cosiddetta regola empirica del pollice che prevede
in ciascuna classe un numero di elementi maggiori di 5.
> attach(data)
> boundaries<-seq(0,max(difetti),3)
> osservate<-table(cut(difetti,boundaries))
> osservate
(0,3] (3,6] (6,9] (9,12] (12,15] (15,18] (18,21] (21,24]
0 5 7 7 5 3 1 2
> boundaries<-c(0,9,12,15,50)
> osservate<-table(cut(difetti,boundaries))
> osservate
(0,9] (9,12] (12,15] (15,50]
12 7 5 6
In tal caso la function da usare è chisq.test() che prende in input le
frequenze osservate ed il vettore della massa di probabilità. Quest’
ultimo viene trasformato nel vettore delle frequenze attese
all’interno della function.
Regola
del
pollice
16/01/2014
20
> osservate
(0,9] (9,12] (12,15] (15,50]
12 7 5 6
REGOLA EMPIRICA DEL
POLLICE
Il valore 0 non è presente nel campione casuale. Posso dunque lasciare inalterato l’intervallo
(0,9].
> attese<-pbinom(boundaries[2],50,0.23)
> attese[2:3]<-pbinom(boundaries[3:4],50,0.23)-pbinom(boundaries[2:3],50,0.23)
> attese[4]<-1-sum(attese[1:3])
> sum(attese)
[1] 1
> chisq.test(osservate, p=attese)
Chi-squared test for given probabilities
data: osservate
X-squared = 9.0582, df = 3, p-value = 0.02853
Warning message:
In chisq.test(osservate, p = attese) :
Chi-squared approximation may be incorrect
Supponiamo di voler verificare per questo vettore se segue una legge binomiale.
14 15 13 14
14 14 13 11
12 13 14 12
12 10 6 11
12 13 12 11
> dati<-c(14,15,13,14,14,14,13,11,12,13,14,12,12,10,6,11,12,13,12,11)
I parametri n e p vengono stimati direttamente dalla function goodfit usando il metodo
di massima verosimiglianza se specifichiamo un certo parametro di input:
> gf<-goodfit(dati, type=“binomial", method='ML')
ESEMPIO:
> summary(gf)
Goodness-of-fit test for binomial distribution
X^2 df P(> X^2)
Likelihood Ratio 10.99785 5 0.0514226
16/01/2014
21
> str(gf)
List of 7
$ observed: num [1:16] 0 0 0 0 0 0 1 0 0 0 ...
$ count : int [1:16] 0 1 2 3 4 5 6 7 8 9 ...
$ fitted : num [1:16] 1.35e-10 9.22e-09 …..
$ type : chr "binomial"
$ method : chr "ML"
$ df : num 5
$ par :List of 2
..$ prob: num 0.82
..$ size: int 15
- attr(*, "class")= chr "goodfit"
> plot(gf$count,gf$observed,type='s',col='red',main='Observed vs Theoretical')
> par(new=TRUE)
> plot(gf$count,gf$fitted,type='l',lwd=4,col='green',ylim=range(0,5),ylab=' ')
>plot(gf)
Serve a capire su quale
osservazione agire per
migliorare il fitting.
Esempio: supponiamo di aver lanciato
un dado e di voler verificare se è stato
truccato. Le frequenze di occorrenza
risultano:
> osservate<-c(18,24,15,25,17,23)
> attese<-rep(1/6,6)
> chisq.test(osservate, p=attese)
Chi-squared test for given probabilities
data: osservate
X-squared = 4.2951, df = 5, p-value = 0.5078
> sum(attese)
[1] 1
> sum(osservate)
[1] 122
16/01/2014
22
Esempio: supponiamo che il numero di prove necessario prima di avere un guasto in
un determinato sistema risulti essere:
Verificare se è possibile assumere valido un modello geometrico.
> osservate<-c(24, 21, 15, 12, 10, 5, 4,2)
> boundaries<-c(1, 2, 3, 4, 5, 6, 8, 9)
> media<-sum(osservate*boundaries)/sum(osservate)
> p<-1/media
> attese<-dgeom(boundaries[1:5],p)
> attese[6]<-dgeom(6,p)+dgeom(7,p)
> attese[7]<-dgeom(8,p)
> attese[8]<-1-sum(attese[1:7])
> sum(attese)
[1] 1
> chisq.test(osservate, p=attese)
Chi-squared test for given probabilities
data: osservate
X-squared = 54.7321, df = 7, p-value = 1.685e-09
Warning message:
In chisq.test(osservate, p = attese) :
Chi-squared approximation may be incorrect
ZZZZ----TEST PER DUE CAMPIONITEST PER DUE CAMPIONITEST PER DUE CAMPIONITEST PER DUE CAMPIONI
Assumiamo che nella produzione della sbarretta metallica dell’azienda
presa in considerazione siano note le varianze nei due diversi periodi
dell’anno. In tal caso, invece del test t conviene usare lo z-test.
Per il primo dataset immaginiamo che la deviazione standard sia 1.25, per il secondo
dataset immaginiamo che la deviazione standard sia 1.
> z.test(prod1, prod2, alternative = "two.sided", sigma.x = 1.25, sigma.y=1)
Two-sample z-Test
data: prod1 and prod2
z = 1.3899, p-value = 0.1646
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval: -0.2040603 1.1990603
sample estimates: mean of x ; mean of y ; 100.2350 99.7375
16/01/2014
23
TTTT----PAIRED TESTPAIRED TESTPAIRED TESTPAIRED TEST
> t.test(prod1, prod2, mu = 0, alternative ="two.sided",paired = FALSE, var.equal =
+ TRUE, conf.level = 0.95)
Quando la variabile PAIRED è posta uguale a TRUE, le unità statistiche sulle quali viene
effettuato il test sono le stesse.
Esempio: Una scuola ha arruolato un nuovo istruttore e vuole verificare l’effica-
cia di un nuovo programma di training proposto, confrontando il tempo medio
di 10 atleti che corrono sui 100 metri.
> a<-c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3)
> b<-c(12.7, 13.6, 12.0, 15.2, 16.8, 20.0, 12.0, 15.9, 16.0, 11.1)
> t.test(a,b,paired=TRUE)
Paired t-test
data: a and b
t = -0.2133, df = 9, p-value = 0.8358
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval: -0.5802549 0.4802549
sample estimates: mean of the differences -0.05
46
Cambiamento della macchina per imballaggio?
Sono statisticamente differenti le due produzioni?
IL CONFRONTO TRA DUE POPOLAZIONI
16/01/2014
24
0 1 2
1 1 2
:
:
H p p
H p p
=
>
Per effettuare un test sulle proporzioni è necessario usare la function
prop.test() . Come parametro di input è necessario costruire una tabella
a doppia entrata contenente il numero di pezzi difettosi riscontrati nel
primo campione e il numero di pezzi difettosi riscontrati nel secondo
campione.
Difettosi Totale
primo 301 1400=28*50
secondo 133 1200=24*50
Totale 480 2600
> sum(value[labels][1:28])
[1] 301
> sum(value[labels][29:52])
[1] 133
TEST SULLE PROPORZIONI
Yates’s correction
> prop.test(x=c(301,133),n=c(1400,1200),correct=FALSE)
2-sample test for equality of proportions without continuity
correction
data: c(301, 133) out of c(1400, 1200)
X-squared = 50.4187, df = 1, p-value = 1.242e-12
alternative hypothesis: two.sided
95 percent confidence interval: 0.07626364 0.13206970
sample estimates: prop 1 prop 2
0.2150000 0.1108333
ANALISI DEI RESIDUI
Nella costruzione di un modello teorico lineare per descrivere la relazione esistente tra due
campioni casuali, era stato tralasciato un punto che rientra nella validazione del modello: ossia
provare che i residui del modello di regressione seguono una legge gaussiana di media nulla
e deviazione standard pari all’errore standard residuo.
Si definiscono residui le distanze
tra i valori della variabile dipen-
dente ottenuti mediante la retta
di regressione e quelli osservati.
ˆi i i
e y y= −
2
Perchè il modello sia valido
è necessario provare che
N(0, )ε σ∼
Per la varianza usiamo il valore del residual
standard error che fornisce la function lm di R
16/01/2014
25
49
Nella tavola è riportata la % di purezza di ossigeno,
rilasciata in un processo di distillazione chimica, e il
livello di idrocarbonio, presente nel condensa-
tore principale di unità di distillazione.
Osservazioni Liv.Idrocarbonio Purezza
1 0,99 90,01
2 1,02 89,05
3 1,15 91,43
4 1,29 93,74
5 1,46 96,73
6 1,36 94,45
7 0,87 87,59
8 1,23 91,77
9 1,55 99,42
10 1,4 93,65
11 1,19 93,54
12 1,15 92,52
13 0,98 90,56
14 1,01 89,54
15 1,11 89,85
16 1,2 90,39
17 1,26 93,25
18 1,32 93,41
19 1,43 94,98
20 0,95 87,33
L’analisi della regressione è una tecnica statistica per modellare e
investigare le relazioni tra due (o più) variabili.
Dati
salvati
in
un file
>result<-lm(datiy ~ datix)
>summary(result)
Residuals:
Min 1Q Median 3Q Max
-1.83029 -0.73334 0.04497 0.69969 1.96809
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 74.283 1.593 46.62 < 2e-16 ***
datix 14.947 1.317 11.35 1. 1.23e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.087 on 18 degrees of freedom
F-statistic: 128.9 on 1 and 18 DF, p-value: 1.227e-09 (hp_0:m=0)
16/01/2014
26
> str(result)…
$ residuals : Named num [1:20] 0.9287 -0.4797 -0.0429 0.1744 0.6234 ...
Dovendo dimostrare che i residui hanno legge gaussiana,
un test più robusto del KS-test è il test di Shapiro-Wilk
> ks.test(result$residuals,"pnorm",0,1.087)
One-sample Kolmogorov-Smirnov test
data: result$residuals
D = 0.0904, p-value = 0.9916
alternative hypothesis: two-sided
> shapiro.test(result$residuals)
Shapiro-Wilk normality test
data: result$residuals
W = 0.9796, p-value = 0.9293
Si basa sul confronto tra la varianza
campionaria e una varianza ottenuta
usando dei pesi particolari.
Il test di Shapiro-Wilks è molto usato per verificare se un campione ca-
suale ha legge gaussiana.
Infatti quando sussiste tale ipotesi, vi sono vari metodi che consentono
di effettuare i test in condizioni di robustezza.
STATISTICA PARAMETRICA
Quando tale ipotesi non è verificata, i metodi cui si ricorre rientrano
nell’ambito della statistica non parametrica.
Ad esempio si consideri un insieme di dati aventi tendenza centrale (8),
ma con legge uniforme.
x<-4+8*runif(90,min=0,max=1)
Come valore di riferimento della tendenza centrale usiamo la mediana
0 0
1 0
:
:
H M M
H M M
=
≠
Ad esempio, vogliamo verificare se la mediana del
dataset così costruito è 6.
Se fosse 6, allora circa una metà dei valori di x do-
vrebbe essere inferiore a 6.
16/01/2014
27
> install.packages('plyr')
> library('plyr')
> count(sign(x-6))
x freq
1 -1 17
2 1 73
>
Pertanto la probabilità che un valore di x sia
inferiore a 6 dovrebbe essere circa il 50%.
Ossia il numero di segni positivi e negativi
del campione casuale dovrebbe essere circa
lo stesso.
Il test che verifica se il numero dei segni negativi nel campione casuale è
circa la metà di tutto il campione è il sign test (o test binomiale).
> binom.test(17,90)
Exact binomial test
data: 17 and 90
number of successes = 17, number of trials = 90, p-value = 1.948e-09
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval: 0.1140597 0.2851131
sample estimates: probability of success 0.1888889
CONFRONTO TRA DUE CAMPIONI IN MANCANZA
DELL’IPOTESI DI POPOLAZIONE GAUSSIANA
Analogo T-test - PAIRED=FALSE Analogo T-test - PAIRED=TRUE
Mann-Whitney test Wilcoxon test
Ipotesi nulla:
i due campioni provengono da una
comune distribuzione di probabilità
Ipotesi alternativa:
c’è un ordinamento stocastico tra
le due distribuzioni di probabilità
che generano i campioni
Ipotesi nulla:
Differenze tra le osservazioni accop-
piate provengono da una distribuzio-
ne di probabilità simmetrica attorno
allo 0.
Ipotesi alternativa:
Differenze tra le osservazioni accop-
piate provengono da una distribuzio-
ne di probabilità asimmetrica.
16/01/2014
28
> wilcox.test(x,y,paired=…,alternative=‘two sided’)
Se paired=FALSE, allora si ritiene che i due campioni sono indipendenti
(Mann-Whitney).
Se paired =TRUE, allora si ritiene che i due campioni siano accoppiati
(Wilcoxon).
Esempio: Un’azienda che produce scarpe da trekking vuole lanciare una
nuova linea. Per questo motivo seleziona 13 scalatori, ad ognuno dei
quali fa indossare una scarpa da trekking della vecchia linea produttiva
ed una scarpa da trekking della nuova linea. Viene poi misurato il nume-
ro di chilometri percorsi prima che inizi il deterioramento della scarpa.
Dopo aver effettuato una analisi descrittiva per entrambi i campioni,
stabilire con un test di ipotesi se c’è differenza.
> old<-c(460,420,520,515,490,490,500,550,480,530,518,515,475)
> new<-c(530,525,500,505,520,450,495,575,474,515,490,480,493)
campione casuale old campione casuale new
Visto che è necessario verificare se entrambi i campioni hanno legge gaussiana, con-
viene usare lo Shapiro-Wilk test.
> shapiro.test(old)
Shapiro-Wilk normality test
data: old
W = 0.9553, p-value = 0.6806
> shapiro.test(new)
Shapiro-Wilk normality test
data: new
W = 0.9634, p-value = 0.8051
16/01/2014
29
Poiché vengono usati sempre gli stessi scalatori per testare le scarpe, ef-
fettuiamo un T-paired test.
> t.test(old,new, mu = 0, alternative ="two.sided", paired = TRUE) Paired t-test
data: old and new
t = -0.5824, df = 12, p-value = 0.5711
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-32.46024 18.76793
sample estimates: mean of the differences -6.846154
Se effettuando lo Shapiro-Wilk test uno dei due campioni fosse risultato
non normale, allora la sintassi da usare sarebbe stata:
> wilcox.test(old,new,paired=TRUE,alternative="two.sided")
Wilcoxon signed rank test
data: old and new
V = 45, p-value = 1
alternative hypothesis: true location shift is not equal to 0