Transcript Variansanalyse

Variansanalyse
på
normalfordelte observationer
af
Jens Friis
Ensidig variansanalyse
Model enkelt normalfordelt observationsrække
Lad X1, X2, ……Xn er indbyrdes uafhængige N(μ, σ2) - fordelt stokastiske variable.
Det tilhørende observationssæt kaldes x1, x2, ……xn
Estimater
n
ˆ 
x. 
x
i 1
i
n
ˆ  s
2
 ( x  x.)
2
i
2

n 1
Kvadratsumsopspaltning
SSD
SSD1
n
 (x
i 1
SSD2
SSD
i
 x.) 2
n( x.   )2
n
 (x
i 1
i
  )2
f
n-1
1
n
Hypotesen
H0 : μ = μ0
H1 : μ ≠ μ0
med
ønskes testet.
Teststørrelsen bliver
n
Det ses, at
 (X
i 1
i
 X.)2
n 1
t
x.  0
s2
n
er en stokastisk variabel, og derfor er t ikke
n
normalfordelt. Man kan vise, at  (Xi  X.)2 er σ2χ2 - fordelt med
i 1
f=n-1 frihedsgrader.
Testoren t følger en såkaldt t-fordeling med f=n-1 frihedsgrader. t-fordelingen
konvergere mod N(0, 1) – fordelingen for n gående mod uendelig. t-fordelingens
tæthedsfunktion er også symmetrisk om 0.
Hypotesen accepteres hvis Tf-1(α/2) ≤ t ≤ Tf-1(1-α/2) , hvor Tf er fordelingdfunktionen svarende til t-fordelingen med f frihedsgrader.
Eksempel: Ved produktion af piller har man målt nicotamid-indholdet i 20 piller.
Indholdet skal være 25mg. Ved stikprøven på 20 piller fik man følgende resultater:
22,67 23,29 23,40 23,56 23,76 23,83 23,95 24,21 24,50 24,64
24,87
25,05
25,35
25,73
25,79
25,80
26,11
26,97
25,36
27,11
Model : Xi ̴ N(μ, σ2) for i=1 til 20 er uafhængige stokastiske variable.
H0 : μ = 25 ,
Parametrene estimeres
x. = 24,797
H1 : μ ≠ 25
;
s2 = 1,5187
Teststørrelsen bliver
t
24,797  25
 0,737
1,5187
20
Da 2,5%’s fraktilen er -2,093 for 19 frihedsgrader, accepters hypotesen.
Anvendelse af SPSS til analysen:
Først undersøges om observationssættet kan
anses for normalfordelt. Man får et såkaldt Q-Q plots
Det accepteres at observationssættet er normalfordelt.
Herefter testes hypotesen : klik Analyze → Compare Means → One-Sample T test
Vælg Test Value til 25
Hypotesen accepteres
Simpel lineær regression
Antag at Yi for i = 1 til k er uafhængige N(μi, σ2) -fordelte således at
i     ( xi  x.)
Man kan vise at estimaterne for parametrene
er
k
ˆ  y. ; ˆ 
(y
i 1
i
 y. )(xi  x. )
k
 (x  x )
i 1
i
.
1 k
ˆ  s 
( yi  y.  ˆ ( xi  x. )) 2

k  2 i 1
2
Bemærk at skæring med y-aksen er
2
y  ˆ  x
Man kan også vise, at estimatoren for β er
N(  ,
2
) -
k
 (x  x )
i 1
i
fordelt.
2
.
Man kan derfor teste hypotesen H0 : β = β0 med teststørrelsen
ˆ   0
t
som er t-fordelt med k-2 frihedsgrader under H0 .
s2
k
 (x  x )
i 1
i
2
.
Hvis β0 = 0 tester man uafhængighed af x og y værdierne.
Kvadratsumsopspaltning :
SSD
Omkring
linje
SSD1
f
n
 ( y  y  ˆ ( x  x ))
2
i
i 1
linje
SSD2
total
SSD
ˆ 2
i
n
 (x  x)
i 1
n
(y
i 1
n-2
i
2
i
1
 y.) 2
n-1
SSD / 1
Som test for H0 : β = 0 an også anvendes SSD /(2n  2) som er F(1,n-2) fordelt.
1
Eksempel : Man for 28 patienter målt kreatininindholdet i blodet før og efter
dødens indtræden. Er der en sammenhæng? Dataene kan ses i en excelfil.
Der er en pæn lineær sammenhæng og parametrene estimeres.
28
ˆ  y  1,024 ; ˆ  1,012 ; ˆ  s  0,01200 ; SSD x   ( xi  x. ) 2  1,4285
2
2
i 1
Man vil gerne teste hypotesen H0 : β = 1
t
1,012 1,000
 0,131
0,01200
1,4285
som er t-fordelt med 26 frihedsgrader. Da 97,5%’s fraktilen er 2,056
accepteres hypotesen.
Dataene er analyseret vha. SPSS : kreatinin.sav
Analyse vha. SPSS
Først undersøges det om der er en lineær sammenhæng:
Dette accepteres.
Parametrene estimeres: Klik Analyze →Regrssion→Linear
s2
Skæring med y-aksen og
ˆ
Testet for H0 : β = 1 bliver
Spredningen på ˆ
t
1.012  1
 0.131
0.092
, som det blev vist tidligere.
Yderligere modelkontrol :
Man bør undersøge residuerne, dvs. afvigelserne fra modellen
Klik Analyze→Regression→Linear→Save og flueben som vist
Optegn de forventede mod de observerede y-værdier mod hinanden
og nogle passende plots af residuerne.
Model flere normalfordelte observationsrækker
Lad Xij , i=1,2…k, j=1,2…ni være indbyrdes uafhængige N(μi, σ2) - fordelt stokastiske
variable.
k
Det tilhørende observationssæt kaldes xij , i=1,2…k, j=1,2…ni, og lad n   ni
i 1
Estimater
ˆ i 
ni
xi . 
x
j 1
ni
ij
k ni
1
ˆ  s02 
( xij  xi .) 2

n  k i 1 j 1
Modelkontrol
Det forudsættes at for hver i er observationsrækken normalfordelt, og
at der er tale om varianshomogenitet for de k observationsrækker dvs.
ni
for
( xij  x i .)2

ˆ i2  si2  j 1
, i=1,2….k
ni  1
Man kan benytte et Barletts test eller et Levene test ( er tilgængeligt i SPSS).
Følgende hypotese ønskes testet:
H0 : μi = μ , i = 1,2…k (samme middelværdi i de k observationsrækker)
Kvadratsumsopspaltning :
SSD
Inden for
grupper
SSD0
Mellem
grupper
SSD1
Total
SSD
ni
k
 ( xij  xi .)2
f
n-k
i 1 j 1
k
 ni ( xi .  x..)2
k-1
i 1
k
ni
 ( x
i 1 j 1
ij
 x..) 2
SSD /(k  1)
n-1
Teststørrelsen for H0 er SSD 1 /(n  k ) , som er F(k-1,n-k) fordelt .
0
Store værdier er kritiske.
Hvis H0 accepteres er estimaterne følgende:
1 k ni
ˆ  x..   xij
n i 1 j 1
1 k ni
ˆ  s 
( xij  x..) 2

n  1 i 1 j 1
2
2
Eksempel
To titreringsmetoder anvendes. Det ønskes undersøgt om de giver samme resultat:
T1
T2
76,35
76,23
Det skal først undersøges om de to observationsrækker kan
76,33
76,30
anses for normalfordelte, og i bekræftende fald om der er
76,45
76,33
varianshomogenitet. Dataene organiseres som liste i SPSS:
76,40
76,33
nr.
Tnr
Antag at dataene er normalfordelte.
1
76,35
76,68
76,28
Klik Analyze → Compare Means →
1
76,33
76,33
76,45
One-way Anova :
1
76,45
76,40
76,38
1
76,40
76,28
76,43
osv.
76,58
76,45
76,65
76,60
76,40
76,40
77,03
76,80
76,90
76,95
74,83
74,88
75,28
75,25
Man får
Da teststørrelsen er 0,014 og
den er F(1, 28) fordelt accepters
hypotesen om varianshomogenitet.
s12
SSD1
SSD0
SSD
s0 2
Test-størrelsen. H0 accepters
( ingen forskel på de to titreringsmetoder).
Tosidig variansanalyse
Model : Xijk ~
N(ij , 2 ) i= 1,2….r ; j=1,2….s ; k=1,2….t ; n=rst
I første omgang skal man undersøge om der er varianshomogenitet i de rs
observationsrækker. Denne hypotese kaldes H0 (arbejdshypotese).
Derefter er der flere hypoteser, som man kan opstille.
H1 :
ij    i  . j
Dvs. en rækkeeffekt plus en søjleeffekt.
H2 :
i  0
Dvs. ingen rækkeeffekt.
H2* :
j  0
Dvs. ingen søjleeffekt.
H3 :
ij  
Dvs. samme fordeling i de rs observationsrækker (fuldstændig homogenitet).
r
Der er valgt en normering således at

i 1
i
 0 og
s

j 1
j
0
.
Man kan vise, at estimaterne for middelværdiparametrene under H1 er :
1 r s t
ˆ  x... 
xijk

rst i 1 j 1 k 1
1 s t
ˆ i  xi ..  x...   xijk  x...
st j 1 k 1
1 r t
ˆ
 j  x. j.  x...   xijk  x...
rt i 1 k 1
Under H0 er estimatet for σ2 : SSD0/f0 ( se næste side)
Under H1 er estimatet for σ2 : (SSD0+SSD1)/(f0+f1 )
Kvadratsumsopspaltning:
SSD
Inden for
grupper
SSD0
Vekselvirkning
SSD1
r
s
f0=rs(t-1)
t
 ( xijk  xij .)2
i 1 j 1 k 1
r
s
 t ( x .  x ..  x. .  x...)
SSD2
2
i
ij
i 1 j 1
Rækkevirkning
f
j
r
 st ( x ..  x...)
2
i
f1=(r-1)(t-1)
f2=r-1
i 1
Søjlevirkning
SSD2*
s
 rt ( x. .  x...)
2
j
j 1
Total
SSD
r
s
t
 ( x
i 1 j 1 k 1
ijk
 x...)2
f2*=t-1
f=rst-1
Test:
H1 : aditivitet
F
SSD 1 / f1
SSD 0 / f 0
som er F ( f1 , f 2)
fordelt.
H2 : ingen rækkevirkning
F
SSD 2 / f 2
(SSD 0  SSD 1 ) /( f 0  f1 )
som er F ( f 2 , f 0  f1 )
fordelt.
H3 : fuldstændig homogenitet (heller ingen søjlevirkning )
F
SSD 2 * / f 2 *
(SSD 0  SSD 1  SSD 2 ) /( f 0  f1  f 2 )
som er F ( f 2 *, f 0  f1  f 2 ) fordelt.
Man kan også vælge at teste for ingen søjlevirkning først. Der skal så byttes rundt
på SSD2 og SSD2* og deres frihedsgrader i de to test. Hver gang man har accepteret en hypotese, er ændres estimatet for variansen. Hvis fx H2 accepteres er
Estimatet for variansen (SSD0+SSD1+SSD2)/(f0+f1+f2)
Eks. Man har testet et byggemateriale for vandgennemtrængning, målt i sekunder.
Man har derpå taget logaritmen til tiden.
Byggematerialet blev produceret på 3 forskellige maskiner 9 forskellige dage med
3 målinger pr. dag:
Først skal man lave en modelkontrol. Da der kun er
tre observationer pr. dag , er det ikke muligt at lave
dag maskine1 maskine2 maskine3
1
1,404
1,306
1,932
en fornuftig kontrol af, om der er tale om normalfordelte
1,346
1,628
1,674
1,618
1,410
1,399
observationer pr. maskine x dag. Derimod kan man
2
1,447
1,241
1,426
1,569
1,185
1,768
estimer variansen pr. maskine x dag, og teste om der er
1,820
1,516
1,859
3
1,914
1,506
1,382
varianshomogenitet. Dette gøres med enten et Bartletts
1,477
1,575
1,690
1,894
1,649
1,361
test eller Levene. I SPSS er det muligt, at foretage et
4
1,887
1,673
1,721
Levene test. For at benytte SPSS skal dataene organiseres
1,485
1,372
1,528
1,392
1,114
1,371
som en lang liste :
5
1,772
1,227
1,320
1,728
1,397
1,489
dag
maskine måling
1,545
1,531
1,336
6
1,665
1,404
1,633
1
1
1,404
1,539
1,452
1,612
1,680
1,627
1,359
1
1
1,346
7
1,918
1,229
1,328
1,931
1,508
1,802
1
1
1,618
2,129
1,436
1,385
8
1,845
1,583
1,689
1
2
1,306
1,790
1,627
2,248
2,042
1,282
1,795
1
2
1,628
9
1,540
1,636
1,703
1,428
1,067
1,370
1
2
1,410
1,704
1,384
1,839
1
3
1,932
1
3
1,674
1
3
1,399
osv.
2
1
1,447
Dette kan gøres samtidigt med den tosidige variansanalyse i SPSS:
Klik Analyze → Generel Linear Model → Univariate og udfyld som vist.
Teststørrelsen er F(26,54) fordelt. Testet er dobbeltsidigt og ikke signifikant her.
Grafisk modelkontrol for additivitet :
Der afsættes punkterne ( xi .., xij .) , i  1..r og ( x. j., xij .) , j  1..s
som skal ligge omkring en ret linje med hældningskoefficienten 1.
Herefter selve variansanalysen:
Her er r=9 , s=3(antal maskiner) og t=3
Er test for H2,
men s22/so2
SSD2
SSD1
SSD0
Test for H1 accept.
SSD
Tosidig variansanalyse med forskelligt antal observationer pr. celle
Model : Xijk ~
r
i 1 j 1
Alt er stort set som før. Man får følgende kvadratsumopspaltning.
SSD
Inden for
grupper
SSD0
Vekselvirkning
SSD1
r
nij
s

r
SSD2
SSD2*
 n ( x .  x ..  x. .  x...)
ij
SSD
2
i
ij
j
r
2
i
i
f2*=t-1
s
 n. ( x. .  x...)
2
j
r
s
j
t
 ( x
i 1 j 1 k 1
ijk
f1=(r-1)(t-1)
f2=r-1
 n .(x ..  x...)
j 1
Total
f0=n-rs
s
i 1
Søjlevirkning
f
( xijk  xij .)2
i 1 j 1 k 1
i 1 j 1
Rækkevirkning
s
N(ij , ) i= 1,2….r ; j=1,2….s ; k=1,2….nij ; n=  n
2
 x...)
2
f=n-1
ij
Lineær regression med flere observationer pr. x
Antag at Yij for i = 1 til k , j=1 til ni er uafhængige N(μij, σ2) -fordelte således at
k
i j     ( xi  x. ) , i  1,2..k , j  1,2..ni
n   ni
i 1
Man kan vise at estimaterne for parametrene
er
k
ˆ  y. ; ˆ 
ni

yij ( xi  x.)
i 1 j 1
k
 n ( x  x.)
i 1
i
x. 
Bemærk at
2
1 k
 ni xi
n i 1
i
Bemærk igen at skæring med y-aksen er y  ˆ  x
Man kan også vise, at estimatoren for β er N(  ,
2
n
 n (x  x )
i 1
i
i
2
) -
fordelt.
.
Man kan derfor teste hypotesen H2 : β = β0 med teststørrelsen
ˆ   0
t
s01
k
2
 ni ( xi  x. ) 2
som er t-fordelt med f0+1 frihedsgrader under H0 .
i 1
Hvis β0 = 0 tester man uafhængighed af x og y værdierne. Vedr. s012 se følgende.
Kvadratsumsopspaltning :
SSD
Inden for
grupper
SSD0
Omkring
linjen
SSD1
Regressionslinjen
SSD2
Total
SSD
k
f
f0=n-k
ni
 ( yij  yi .)2
i 1 j 1
k
 n ( y .  y..  ˆ ( x  x.))
2
i
i 1
ˆ 2
i
i
k
 ni ( xi  x.)2
f1=k-2
f2=1
i 1
k
ni
 ( y
i 1 j 1
Testet for H1 : lineær regression er
ij
 y..) 2
f=n-1
SSD 1 /(k  2)
som er F(k-2,n-k) fordelt.
SSD 0 /(n  k )
Bemærk, at hvis H1 accepteres er estimatet for variansen s012=(SSD0+SSD1)/(f0+f1)
SSD 2 / 1
Testet for H2: β = 0 fuldstændig homogenitet er (SSD  SSD ) /(n  2)
0
1
som er F(1, n-2) fordelt.
Modelkontrol:
Det skal undersøges, at for hvert k kan observarionsrækken yij, j=1,2..ni
anses for normalfordelt
Eksempel:
Nedenstående tabel viser logaritmen til trækstyrken (kg/cm2) og den reciprokke
hærdningstid ( dage) for nogle cementstykker:
dage måling nr. Træk.styrke log
reciprok dag
1
1
13,00
1,114
1,000
1
2
13,30
1,124
1,000
1
3
11,80
1,072
1,000
2
1
21,90
1,340
0,500
2
2
24,50
1,389
0,500
2
3
24,70
1,393
0,500
3
1
29,80
1,474
0,333
3
2
28,00
1,447
0,333
3
3
24,10
1,382
0,333
3
4
24,20
1,384
0,333
3
5
26,20
1,418
0,333
7
1
32,40
1,511
0,143
7
2
30,40
1,483
0,143
7
3
34,50
1,538
0,143
7
4
33,10
1,520
0,143
7
5
35,70
1,553
0,143
28
1
41,80
1,621
0,036
28
2
42,60
1,629
0,036
28
3
40,30
1,605
0,036
28
4
35,70
1,553
0,036
Først en grafisk undersøgelse:
28
5
37,30
1,572
0,036
Som det ses er der tale om en pæn lineær
Sammenhæng.
Lad yij betegne log(trækstyrke) og xi den
reciprokke hærdningstid. n = 21, k = 5 0
Klik Analyze → Compare Means →
One-Way Anova →
Accept af varianshomogenitet.
SSD0
Herefter skal der foretages en lineær regression.
Tast Analyze → Regression → Linear og man får
SSD0+SSD1
ˆ
Test for linearitet
F
Skæring med y-aksen
SSD 1/f1 (0,020466 0,016808) / 3

 1,16
SSD 0 /f 0
0,016808/ 16
som accepteres.
Videregående regressionsanalyse :
Model:
Antag at Yi for i = 1 til k er uafhængige N(μi, σ2) -fordelte således at
p
i   xij  j ,hvor xij’erne er kendte værdier og βj’erne ukendte parametre.
j 1
Dette kan formuleres med matricer:
 1   x11
  
  2   x21
 .  .
  
 .   .
  x
 k   k1
Og lad
x12
x22
.
.
xk 2

. . x1 p  1 
 
. . x2 p   2 
. . .  .   Xβ
 
. . .  . 
. . xkp   p 
L1  Xβ  R k : β  R p
βˆ  (X' X)1 X' y
og lad
 y1 
 
 y2 
y . 
 
 . 
y 
 k
betegne observationerne.
 være et underrum. Estimaterne bliver
ˆ  s 
2
2
y  Xβˆ
n  dim L1
Ofte sættes første søjle i X til 1-taller således, at β1 er det generelle niveau.
Eksempel : Indianere i Peru
Ændringer i menneskers livsbetingelser kan give sig udslag i fysiologiske ændringer,
eksempelvis i ændret blodtryk.
En gruppe antropologer undersøgte hvordan blodtrykket ændrer sig hos peruvianske
indianere der flyttes fra deres oprindelige primitive samfund i de høje Andesbjerge til den
såkaldte civilisation, dvs. storbyen, der i øvrigt ligger i langt mindre højde over havets
overflade end deres oprindelig bopæl (Davin (1975), her citeret e er Ryan et al. (1976)).
Antropologerne udvalgte en stikprøve på 39 mænd over 21 år der havde undergået en
sådan flytning. På hver af disse måltes blodtrykket (det systoliske og det diastoliske) samt
en række baggrundsvariable, heriblandt alder, antal år siden flytningen, højde, vægt og
puls. Desuden har man udregnet endnu en baggrundsvariabel, nemlig »brøkdel af livet
levet i de nye omgivelser«, dvs. antal år siden flytning divideret med nuværende alder.
Man forestillede sig at denne baggrundsvariabel kunne have stor »forklaringsevne«.
Her vil vi ikke se på hele talmaterialet, men kun på blodtrykket (det systoliske) der
skal optræde som y-variabel, og på de to x-variable brøkdel af livet i de nye omgivelser og
vægt. Disse er angivet i tabel 11.8 (fra Ryan et al. (1976)).
1. Antropologerne mente at x2, brøkdel levet i de nye omgivelser, var et godt mål for
hvor længe personerne havde levet i de civiliserede omgivelser, og at det derfor
måtte være interessant at se om x2 kunne forklare variationen i blodtrykket y. Første
skridt kunne derfor være at estimere en simpel lineær regressionsmodel med x2
som forklarende variabel. Gør det!
2. Hvis man i et koordinatsystem afsætter y mod x2, viser det sig imidlertid at det faktisk
ikke virker særlig rimeligt at hævde at (middelværdien af) y afhænger lineært
af x2. Derfor må man give sig til at overveje om andre af de målte baggrundsvariable
med fordel kan inddrages.
Nu ved man at en persons vægt har betydning for den pågældendes blodtryk,
så næste modelforslag kunne være en multipel regressionsmodel med både x2 og
x3 som forklarende variable.
I SPSS indtastes dataene således: (hvis man ikke havde 1-tallene vil SPSS give det samme)
y
170
120
125
148
140
x1
1
1
1
1
1
x2
0,048
0,273
0,208
0,042
0,040
x3
71,0
56,5
56,0
61,0
65,0
Tast Analyze → Regression → Linear
Osv.
Eksempel : Indianerne i Peru ( se opgaveark)
s2
 ' erene
Alle test for βi = 0 er signifikante.
test for lig 0
Modelkontrol : Der laves først simple grafer over sammenhæng mellem y’erne
og x2’erne og derpå x3’erne. Der er ikke overbevisende lineær sammenhæng.
Parametrene i den multiple regression estimeres og de forventede værdier og
residuerene beregnes :klik yderligere på Save og sæt flueben somvist.
Sammenhænget mellem
forventet og observeret
er ikke overbevisende
men acceptabelt.
Residuerene undersøges:
Det accepteres, at
residuerne kan anses
for normalfordelte,
men det er ikke flot.