ALGEMENE STATISTIEK. A.W. van der Vaart en anderen

ALGEMENE STATISTIEK A.W. van der Vaart en anderen

VOORWOORD Dit diktaat wordt gebruikt bij het vak Statistiek voor Natuurkunde. Het is een uittreksel van het boek ”Algemene Statistiek” geschreven door A.W. van der Vaart en anderen. Enkele voorbeelden en paragrafen zijn later toegevoegd.

INHOUD 1. Inleiding . . . . . . . . . . . . . . . . . . . 1.1. Wat is statistiek? . . . . . . . . . . . . . 1.2. Beschrijvende versus Mathematische Statistiek 1.3. Indeling van het boek . . . . . . . . . . . 2. Statistische Modellen . . . . . . . . . . . . . 2.1. Introductie . . . . . . . . . . . . . . . . 2.2. Enkele voorbeelden . . . . . . . . . . . . Opgaven . . . . . . . . . . . . . . . . . 3. Verdelingsonderzoek . . . . . . . . . . . . . . 3.1. Introductie . . . . . . . . . . . . . . . . 3.2. Univariate Steekproeven . . . . . . . . . . 3.3. Samenhang . . . . . . . . . . . . . . . . Opgaven . . . . . . . . . . . . . . . . . 4. Schatters . . . . . . . . . . . . . . . . . . 4.1. Introductie . . . . . . . . . . . . . . . . 4.2. Maximum Likelihood-Schatters . . . . . . . 4.3. Momentenschatters . . . . . . . . . . . . Opgaven . . . . . . . . . . . . . . . . . 5. Toetsen . . . . . . . . . . . . . . . . . . . 5.1. Nulhypothese en Alternatieve Hypothese . . . 5.2. Toetsingsgrootheid en Kritiek Gebied . . . . 5.3. Statistische Significantie . . . . . . . . . . 5.4. Overschrijdingskansen . . . . . . . . . . . 5.5. Enkele Standaard Toetsen . . . . . . . . . Opgaven . . . . . . . . . . . . . . . . . 6. Betrouwbaarheidsgebieden . . . . . . . . . . . 6.1. Introductie . . . . . . . . . . . . . . . . 6.2. Pivots en Bijna-Pivots . . . . . . . . . . . 6.3. Maximum Likelihood-Schatters als Bijna-Pivots Opgaven . . . . . . . . . . . . . . . . . 7. Enkele Regressiemodellen . . . . . . . . . . . 7.1. Lineaire Regressie . . . . . . . . . . . . . 7.2. Niet-lineaire en niet-parametrische regressie . Opgaven . . . . . . . . . . . . . . . . . 8. Appendix A: Elementen uit de Kansrekening . . . 8.1. Verdelingen . . . . . . . . . . . . . . . 8.2. Verwachting en variantie . . . . . . . . . . 8.3. Standaard verdelingen . . . . . . . . . . . 8.4. Multivariate en marginale verdelingen . . . . 8.5. Onafhankelijkheid en conditionering . . . . . 8.6. Limietstellingen en de normale benadering . . Opgaven . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

1 1 1 2 3 3 3 9 10 10 10 13 15 16 16 17 22 23 24 24 25 32 33 35 45 48 48 50 53 56 57 58 68 69 71 71 73 74 76 77 79 80

1 Inleiding

1.1

Wat is statistiek?

Statistiek is de kunst van het modelleren van situaties waarin toeval een rol speelt, en van het trekken van conclusies op basis van data waargenomen in dergelijke situaties. Enkele typerende vragen die met behulp van statistiek kunnen worden beantwoord zijn: (i) Wat is de kans dat de Maas komend jaar buiten zijn oevers treedt? (ii) Is de nieuwe medische behandeling significant beter dan de oude? (iii) Wat is de onzekerheidsmarge in de voorspelling van het aantal zetels voor politieke partij A? Het beantwoorden van dergelijke vragen is verre van eenvoudig. De mathematische statistiek levert een algemeen kader waarmee de onderzoeksvraag beantwoord kan worden op basis van een opgesteld statistisch model. Binnen dit kader geeft het ook een oordeel over de kwaliteit van een gegeven antwoord. Om een geschikt statistisch model voor beschikbare data op te stellen, moet inzicht verkregen worden in de manier waarop de data verzameld zijn. Wanneer er nog geen data beschikbaar zijn, zullen die moeten worden verzameld. Het verkrijgen van relevante data vereist een goede, doordachte opzet. Zo zal bij een onderzoeksvraag die een bepaalde populatie betreft (bijvoorbeeld de populatie van patiënten met een hoge bloeddruk, stemgerechtigden, of eindprodukten van een productieproces), data verzameld moeten worden van een groep “mensen” die representatief is voor de gehele populatie. Ten slotte moet dan een geschikt statistisch model worden opgesteld voor de data. De vragen (i)–(ii)–(iii) corresponderen met de drie basis concepten in de statistiek: schatten, toetsen en betrouwbaarheidsgebieden, welke uitgebreid aan de orde komen in dit boek. De nadruk ligt in dit boek op de mathematische statistiek; het verzamelen van data, het vervolgens modelleren van de data, en beschrijvende statistiek komen slechts summier aan bod.

1.2

Beschrijvende versus Mathematische Statistiek

Waarnemingen, meestal rijen getallen, kan men middelen, tabelleren, grafisch weergeven, of anderszins bewerken. De beschrijvende statistiek houdt zich bezig met het verzamelen en op inzichtelijke wijze samenvatten van data. Zulke beschrijvende statistiek, op grote schaal beoefend door bijvoorbeeld het Centraal Bureau voor de Statistiek, is van groot belang en kan heel interessant zijn. Beschrijvende statistiek wordt ook veel gebruikt bij het opstellen van statistische modellen (zie Hoofdstuk 2) en het controleren van modelaannames (zie Hoofdstuk 3). In dit boek komt zij echter nauwelijks aan de orde.

2

1: Inleiding

De mathematische statistiek ontwikkelt en bestudeert methoden voor het analyseren van waarnemingen, die gebaseerd zijn op kansmodellen. Waarneming x wordt opgevat als een realisatie van een stochastische grootheid of vector X. In de waarschijnlijkheidsrekening wordt een precieze definitie gegeven van stochastische vectoren. Voor de statistiek is vooral van belang dat een stochastische vector een kansverdeling bezit. Deze kan worden vastgelegd door een verdelingsfunctie of kansdichtheid. In de statistiek willen we op grond van de realisatie x de “ware” kansverdeling van X bepalen. Op grond van kennis van die ware kansverdeling kunnen we vervolgens nieuwe uitkomsten voorspellen, of oude uitkomsten verklaren.

1.3

Indeling van het boek

De drie kernpunten van de mathematische statistiek zijn schatten, toetsen en het construeren van betrouwbaarheidsgebieden. Deze onderwerpen komen achtereenvolgens aan de orde in de Hoofdstukken 4, 5 en 6. Deze concepten maken gebruik van een statistisch model voor de data, waarvan in Hoofdstuk 2 de definitie en een aantal voorbeelden worden gegeven. Enkele technieken uit de beschrijvende statistiek die hulp kunnen bieden bij het opstellen en valideren van statistische modellen worden besproken in Hoofdstuk 3. In Hoofdstuk 7 worden enkele regressiemodellen die in de praktijk veel gebruikt worden beschreven. De theorie uit de voorgaande hoofdstukken wordt hierin toegepast om onbekende modelparameters te schatten, te toetsen en betrouwbaarheidsintervallen voor deze parameters op te stellen. In Appendix 8 wordt een aantal elementen uit de kansrekening behandeld die van belang zijn voor het begrip van de stof in het boek.

2 Statistische Modellen

2.1

Introductie

In zekere zin is de richting van de statistiek precies de omgekeerde van die van de waarschijnlijkheidsrekening: de uitslagen van een experiment zijn waargenomen, maar het onderliggende kansmodel is (deels) onbekend en dient uit de uitslagen te worden afgeleid. Uiteraard is de experimentele situatie niet geheel onbekend. Alle bekende informatie wordt gebruikt om een zo goed mogelijk statistisch model te construeren. Een formele definitie van een “statistisch model” is als volgt. Definitie 2.1. Een statistisch model is een collectie van kansverdelingen op een gegeven uitkomstenruimte. De interpretatie van een statistisch model is: de collectie van alle mogelijk geachte kansverdelingen voor de waarneming X. Hierin is X het geheel van de waarnemingen. Meestal is deze totale waarneming opgebouwd uit “deelwaarnemingen” en is X = (X1 , . . ., Xn ) een stochastische vector. Wanneer de variabelen X1 , . . ., Xn corresponderen met onderling onafhankelijke replicaties van hetzelfde experiment, dan spreken we van een steekproef. De variabelen X1 , . . ., Xn zijn dan onderling onafhankelijk en identiek verdeeld en hun simultane verdeling wordt volledig bepaald door de marginale verdeling, die voor alle Xi ’s gelijk is. In dat geval kan het statistische model voor X = (X1 , . . ., Xn ) worden beschreven door een collectie van (marginale) kansverdelingen voor de deelwaarnemingen X1 , . . ., Xn .

2.2

Enkele voorbeelden

Het begrip “statistisch model” wordt pas echt duidelijk door voorbeelden. Zo eenvoudig als het wiskundige begrip “statistisch model” is uitgedrukt in de voorgaande definitie, zo ingewikkeld is het proces van statistisch modelleren van een gegeven praktijksituatie. Het resultaat van een statistisch onderzoek staat of valt echter met het construeren van een goed model. Voorbeeld 2.2 (Steekproef ). Van een grote populatie bestaande uit N personen heeft een onbekende fractie p een bepaalde eigenschap A; deze fractie p willen we “schatten”. Het wordt als te veel moeite beschouwd om alle personen uit de populatie op eigenschap A te onderzoeken. In plaats daarvan kiest men aselect n personen uit de populatie, met teruglegging. Men neemt

4

2: Statistische Modellen

(een realisatie van) de stochastische grootheden X1 , . . ., Xn waar, waarbij 0 als de ie persoon A niet heeft, Xi = 1 als de ie persoon A wel heeft. Vanwege de manier waarop het experiment is ingericht (trekken met teruglegging) weten we a priori dat X1 , . . ., Xn onderling onafhankelijk en alternatief verdeeld zijn. Dat laatste wil zeggen dat P(Xi = 1) = 1 − P(Xi = 0) = p voor i = 1, . . ., n. Over de parameter p is geen a priori kennis beschikbaar, anders dan dat 0 ≤ p ≤ 1. De totale waarneming is hier de vector X = (X1 , . . ., Xn ). Het statistische model voor X bestaat uit alle mogelijk geachte (simultane) kansverdelingen van X waarvan de co¨ ordinaten, X1 , . . ., Xn , onderling onafhankelijk en alternatief verdeeld zijn. Voor iedere mogelijke waarde van p bevat het statistische model precies één kansverdeling voor X. Het ligt voor de hand de onbekende p te “schatten” met de fractie van het aantal personen Pn met eigenschap A; dus met n−1 i=1 xi waarbij xi gelijk is aan 1 of 0 al naar gelang de persoon eigenschap A wel of niet heeft. In Hoofdstuk 4 geven we een precieze betekenis aan “schatten”. In Hoofdstuk 6 gebruiken we het zojuist beschreven model om te kwantificeren hoever deze schatter van p zal afwijken, met behulp van een “betrouwbaarheidsinterval”. Bijna nooit zullen de populatie- en steekproeffractie immers precies gelijk zijn. Een betrouwbaarheidsinterval geeft een precieze betekenis aan de “foutenmarge” die vaak bij de uitslag van een opiniepeiling wordt vermeld. We zullen ook berekenen hoe groot die marge is wanneer we bijvoorbeeld 1000 personen uit de populatie onderzoeken, een gebruikelijk aantal bij een opiniepeiling onder de Nederlandse bevolking. Voorbeeld 2.3 (Meetfouten). Als een fysicus middels een experiment herhaaldelijk de waarde van een constante µ bepaalt, vindt hij niet steeds dezelfde waarde. Zie bijvoorbeeld Figuur 2.1, waarin de 23 bepalingen van de lichtsnelheid door Michelson in 1882 zijn weergegeven. De vraag is hoe de onbekende constante µ op grond van de waarnemingen, een rij getallen x1 , . . ., xn , “geschat” kan worden. Voor de waarnemingen in Figuur 2.1 zal deze schatting in de range 700–900 liggen, maar de vraag is waar. Een statistisch model verleent houvast bij het beantwoorden van deze vraag. Kansmodellen zijn in deze context voor het eerst toegepast aan het eind van de 18e eeuw, en de normale verdeling werd door Gauss rond 1810 “ontdekt” precies met het doel inzicht te geven in deze situatie.

600

700

800

900

1000

Figuur 2.1. Grafische weergave van de resultaten van de 23 metingen van de lichtsnelheid door Michelson in 1882. De schaal op de horizontale as geeft de gemeten lichtsnelheid (in km per seconde) min 299000 km/sec.

Als de metingen steeds onder dezelfde omstandigheden worden verricht, steeds onafhankelijk van het verleden, dan is het redelijk in het model op te nemen dat deze getallen realisaties zijn van onderling onafhankelijke en identiek verdeelde stochastische variabelen X 1 , . . ., Xn .

2.2: Enkele voorbeelden

5

De meetfouten ei = Xi − µ zijn dan eveneens stochastische variabelen. Een gebruikelijke aanname is dat de verwachte meetfout gelijk is aan 0, met andere woorden Eei = 0, en dus is EXi = E(ei + µ) = µ. Aangezien wordt aangenomen dat X1 , . . ., Xn onafhankelijke stochastische variabelen zijn en dezelfde kansverdeling hebben, ligt het model voor X = (X1 , . . ., Xn ) vast als we een statistisch model voor Xi kiezen. Voor Xi postuleren we het model: alle kansverdelingen met eindige verwachting µ. Het statistische model voor X is dan: alle mogelijke kansverdelingen van X = (X1 , . . ., Xn ) zodanig dat de co¨ ordinaten X1 , . . ., Xn onderling onafhankelijk en identiek verdeeld zijn met verwachting µ. Fysici menen vaak meer a priori informatie te bezitten, en doen dan meer modelaannames. Ze veronderstellen bijvoorbeeld dat de meetfouten normaal verdeeld zijn met verwachting 0 en variantie σ 2 , ofwel dat de waarnemingen X1 , . . ., Xn normaal verdeeld zijn met verwachting µ en variantie σ 2 . Het statistische model is dan gelijk aan: alle kansverdelingen van X = (X1 , . . ., Xn ) zodanig dat de co¨ ordinaten onderling onafhankelijk en N (µ, σ 2 )-verdeeld zijn. Het uiteindelijke doel is iets te zeggen over µ. Bij het tweede model is meer bekend, dus moet het mogelijk zijn met meer “zekerheid” iets over µ te zeggen. Anderzijds is er natuurlijk meer “kans” dat het tweede model onjuist is, in welk geval de winst aan zekerheid slechts een schijnzekerheid is. In de praktijk blijken meetfouten vaak, maar niet altijd, bij benadering normaal verdeeld te zijn. Dergelijke normaliteit is te motiveren met behulp van de Centrale Limietstelling (zie Stelling 8.21) indien een meetfout kan worden opgevat als de som van een groot aantal onafhankelijke kleine meetfouten (met eindige varianties), maar kan niet op theoretische gronden worden bewezen. In Hoofdstuk 3 bespreken we technieken om normaliteit aan de data zelf te onderzoeken. Het belang van een precies omschreven model is onder andere dat het mogelijk maakt te bepalen wat een zinvolle manier is om µ uit de waarnemingen te “schatten”. Het middelen van x1 , . . ., xn ligt voor de hand. We kunnen laten zien dat dit het beste is (volgens een bepaald criterium) als de meetfouten inderdaad een normale verdeling volgen met verwachting 0. Zouden de meetfouten echter Cauchy-verdeeld zijn, dan is middelenP desastreus. Dit blijkt uit Figuur 2.2. n Deze toont voor n = 1, 2, . . ., 1000 het gemiddelde n−1 i=1 xi van de eerste n realisaties x1 , . . ., x1000 van een steekproef uit een standaard Cauchy-verdeling. De gemiddelden gedragen zich chaotisch en komen niet steeds dichter bij P 0. Dit kan worden verklaard uit het opmerkelijke n theoretische resultaat dat het gemiddelde n−1 i=1 Xi van onderling onafhankelijke standaard Cauchy-verdeelde stochastische grootheden X1 , . . ., Xn zelf ook standaard Cauchy-verdeeld is. Middelen doet hier niets!

2

1

0

-1

-2

-100

100

300

500

700

900

1100

Figuur 2.2. Cumulatieve gemiddelden (verticale as) van n = 1, 2, . . ., 1000 (horizontale as) realisaties uit de standaard Cauchy-verdeling.

Voorbeeld 2.4 (Gepaarde en ongepaarde waarnemingen). De laatste jaren is het aantal verschillende diëten op de markt sterk toegenomen. Om de effectiviteit van diëten A en B met elkaar te vergelijken wordt een aselecte groep zwaarlijvige mensen geheel willekeurig in twee groepen ter grootte n en m verdeeld. De mensen in de eerste groep volgen dieet A en

6


de mensen in de tweede groep dieet B. Na een halfjaar tijd wordt genoteerd hoeveel elke deelnemer is afgevallen. Voor de groep mensen die dieet A volgden, geeft dat de waarnemingen x1 , . . ., xn , waarbij xi de gewichtsafname van de ie persoon in de eerste groep voorstelt. Voor de tweede groep worden de gewichtsafnames genoteerd met y1 , . . ., ym . De waarden x1 , . . ., xn kunnen worden gezien als de realisaties van n onderling onafhankelijke en identiek verdeelde stochastische grootheden X1 , . . ., Xn . Als statistisch model van Xi nemen we alle mogelijke continue kansverdelingen op R. Daarmee sluiten we bij voorbaat een eventuele toename in gewicht niet uit. Het statistische model voor X = (X1 , . . ., Xn ) ligt nu vast. Eveneens kunnen y1 , . . ., ym worden gezien als realisaties van stochastische variabelen Y1 , . . ., Ym welke onderling onafhankelijk en gelijk verdeeld zijn. Het statistische model voor Y = (Y1 , . . ., Ym ) nemen we analoog aan het model voor X. Om de twee diëten met elkaar te vergelijken kunnen de gemiddelde gewichtsafnames in de twee groepen met elkaar worden vergeleken. Met deze opzet van het onderzoek worden twee datasets die op geen enkele manier afhankelijk van elkaar zijn vergeleken; immers de groep zwaarlijvigen was aselect getrokken en geheel willekeurig in twee groepen verdeeld. Soms heeft het zin om de data opzettelijk afhankelijk van elkaar te maken, bijvoorbeeld door mensen te paren. Een reden om waarnemingen te paren kan zijn dat er meer factoren zijn die mogelijk invloed hebben op de uitkomst, gewichtsafname in dit voorbeeld. Corrigeren voor het effect van deze factoren kan de onderzoeksresultaten betrouwbaarder maken. In dit voorbeeld hebben geslacht en begingewicht mogelijk invloed op de gewichtsafname. Om hier rekening mee te houden bij het vergelijken van de twee diëten, worden de personen in de steekproef in n groepjes van twee gedeeld; de mensen worden gepaard. De twee personen in elk paar zijn van hetzelfde geslacht en hebben ongeveer hetzelfde (begin)gewicht. Van elk groepje volgt één persoon dieet A en de andere persoon dieet B; wie welk dieet volgt, wordt geheel willekeurig bepaald. Na een halfjaar wordt gekeken hoeveel elke persoon is afgevallen; dit geeft de waarnemingsparen (x1 , y1 ), . . ., (xn , yn ) waarbij xi de gewichtsafname van de persoon in het ie paar is die dieet A volgde en yi de gewichtsafname van de persoon in hetzelfde paar die dieet B volgde. Omdat we ge¨ınteresseerd zijn in verschil in effectiviteit tussen de twee diëten, ligt het voor de hand om naar de verschillen z1 = x1 − y1 , . . ., zn = xn − yn te kijken en hier een statistisch model voor op te stellen in plaats van voor de gehele dataset. De verschillen z 1 , . . ., zn worden weer gezien als realisaties van onafhankelijke en gelijk verdeelde stochastische grootheden Z1 , . . ., Zn . We nemen als (marginaal) statistisch model voor Zi alle mogelijke continue kansverdelingen op R. Omdat Z1 , . . ., Zn onderling onafhankelijk en identiek verdeeld zijn, ligt het statistische model van Z = (Z1 , . . ., Zn ) daarmee vast. Met deze tweede onderzoeksmethode worden personen gepaard op geslacht en begingewicht; we spreken dan van gepaarde waarnemingen. Bij de eerste methode was er geen sprake van paren en hadden we te maken met ongepaarde data. Een alternatief statistisch model dat ook rekening kan houden met het effect van geslacht en begingewicht is een zogenaamd regressiemodel. Een regressiemodel kan eenvoudig worden uitgebreid, zodat met nog meer factoren rekening kan worden gehouden. Het regressiemodel komt in Voorbeeld 2.5 en in Hoofdstuk 7 aan de orde. Voorbeeld 2.5 (Regressie). Lange ouders krijgen over het algemeen lange kinderen en korte ouders, korte kinderen. De lengte van de ouders hebben een grote voorspellende waarde voor de zogenaamde eindlengte van hun kinderen, de lengte als kinderen zijn uitgegroeid. Er zijn meer factoren die invloed hebben. Het geslacht van het kind speelt natuurlijk een belangrijke rol. Ook omgevingsfactoren als gezonde voeding en hygiëne zijn van belang. Door verbeterde voeding en een toegenomen hygiëne in de afgelopen 150 jaar hebben factoren die de lengtegroei belemmeren, als infectieziekten en ondervoeding, minder kans gekregen in de meeste Westerse landen. Hierdoor is de gemiddelde lichaamslengte toegenomen en worden kinderen elke generatie langer. De streeflengte (of “target height”) van een kind is de eindlengte die kan worden verwacht op basis van de lengte van de ouders, het geslacht van het kind en de toename van lichaamslengte over generaties. De vraag luidt op welke manier de streeflengte afhangt van deze factoren. Definieer Y als de eindlengte die een kind zal bereiken, x1 en x2 als de lengte van de

2.2: Enkele voorbeelden

7

biologische vader en moeder, en x3 als een indicator voor het geslacht (x3 = −1 voor een meisje en x3 = 1 voor een jongen). De streeflengte EY wordt gemodelleerd met een zogenaamd lineair regressiemodel EY = β0 + β1 x1 + β2 x2 + β3 x3 , waarbij β0 de toename van de gemiddelde lichaamslengte per generatie is, β1 en β2 de mate waarin de lengte van de ouders invloed hebben op de streeflengte van hun nageslacht en β 3 is de afwijking van de streeflengte tot de gemiddelde volwassen lengte die wordt veroorzaakt door het geslacht van het kind. Aangezien mannen gemiddeld langer zijn dan vrouwen zal β 3 positief zijn. Bovenstaand lineair model zegt niets over individuele lengtes, maar enkel over dat van het nageslacht van ouders met een bepaalde lengte. Zo hebben twee broers dezelfde streeflengte; ze hebben immers dezelfde biologische ouders, hetzelfde geslacht en zijn geboren in dezelfde generatie. De werkelijke eindlengte Y kan geschreven worden als Y = β0 + β1 x1 + β2 x2 + β3 x3 + e, waarbij e = Y − EY de afwijking is van de werkelijke eindlengte Y ten opzichte van de streeflengte EY . De waarneming Y wordt ook wel de afhankelijke variabele genoemd en de variabelen x1 , x2 en x3 de onafhankelijke of verklarende variabelen. Veelal wordt verondersteld dat e normaal verdeeld is met verwachting 0 en onbekende variantie σ 2 . De eindlengte Y heeft dan een normale verdeling met verwachting β0 + β1 x1 + β2 x2 + β3 x3 en variantie σ 2 . In Nederland wordt periodiek de lengtegroei van de jeugd in kaart gebracht. In 1997 vond de Vierde Landelijke Groeistudie plaats. Een onderdeel van het onderzoek betrof de relatie tussen de eindlengte van kinderen en de lengte van hun ouders. Om deze relatie te bepalen waren gegevens verzameld van jongvolwassenen en hun ouders. Dit leverde de volgende waarnemingen: (y1 , x1,1 , x1,2 , x1,3 ), . . .(yn , xn,1 , xn,2 , xn,3 ) op, waar yi de lichaamslengte van de ie jongvolwassene is, xi,1 en xi,2 de lengte van de biologische ouders, en xi,3 een indicator voor het geslacht van de ie jongvolwassene. Veronderstel dat de waarnemingen onafhankelijke replicaties zijn uit bovenstaand lineair regressiemodel; dat wil zeggen dat gegeven x i,1 , xi,2 , en xi,3 , Yi verwachting β0 + β1 xi,1 + β2 xi,2 + β3 xi,3 en variantie σ 2 heeft. De parameters (β0 , β1 , β2 , β3 ) zijn onbekend en kunnen geschat worden op basis van de waarnemingen. Voor een eenvoudige interpretatie van het model is er voor gekozen om β1 = β2 = 1/2 te nemen, zodat de streeflengte gelijk is aan de gemiddelde ouderlengte gecorrigeerd voor het geslacht van het kind en de invloed van de tijd. De parameters β0 en β3 zijn gelijk aan de toename van de lichaamslengte in de afgelopen generatie en de helft van het gemiddelde lengteverschil tussen mannen en vrouwen. Deze parameters werden geschat met behulp van de kleinste kwadratenmethode. De parameter β0 is geschat met 4.5 centimeter en β3 met 6.5 centimeter. Het geschatte regressiemodel is dan gelijk aan (2.1)

1 Y = 4.5 + (x1 + x2 ) + 6.5x3 + e. 2

In Figuur 2.3 is de lichaamslengte van 44 jongvolwassen mannen (links) en 67 jongvolwassen vrouwen (rechts) uitgezet tegen de gemiddelde lichaamslengte van hun ouders. † De lijn is gelijk aan de geschatte regressielijn gevonden in de Vierde Landelijke Groeistudie. Het geschatte regressiemodel dat gevonden werd in de Vierde Landelijke Groeistudie, kunnen we gebruiken voor het voorspellen van de eindlengte van kinderen die nu geboren worden. We moeten dan wel veronderstellen dat de lengtetoename de komende generatie opnieuw 4.5 centimeter is en het gemiddelde lengteverschil tussen mannen en vrouwen gelijk aan 13 centimeter blijft. Op basis van het bovenstaande model zijn de streeflengten voor zonen en dochters van een man met een lengte van 180 centimeter en een vrouw van 172 centimeter gelijk aan 4.5 + (180 + 172)/2 + 6.5 = 187 centimeter en 4.5 + (180 + 172)/2 − 6.5 = 174 centimeter. † Bron: De data zijn verzameld door de afdeling Biologische Psychologie van de Vrije Universiteit in het kader van een onderzoek naar gezondheid, levensstijl en persoonlijkheid.


170

160

175

165

180

170

185

175

190

180

195

185

200

8

165

170

175

180

185

165

170

175

180

185

Figuur 2.3. Lengte van zonen (links) en dochters (rechts) uitgezet tegen de gemiddelde lichaamslengte van hun ouders. De lijn is de regressielijn gevonden in de Vierde Landelijke Groeistudie.

In andere Europese landen worden andere modellen gebruikt. In Zwitserland, bijvoorbeeld, is de streeflengte gelijk aan EY = 51.1 + 0.718

x1 + x 2 + 6.5x3 . 2

Nu is de streeflengte van de zonen en dochters van ouders met dezelfde lengte als in het voorbeeld hiervoor gelijk aan 184 en 171 centimeter. In het voorgaande voorbeeld bestaat er een lineair verband tussen de respons Y en de onbekende parameters β0 , . . ., β3 . In dat geval spreken we van een lineair regressiemodel. Het meest eenvoudige lineaire regressiemodel is het model waarbij er slechts één verklarende variabele is: Y = β0 + β1 x + e; het enkelvoudige lineaire regressiemodel (in tegenstelling tot meervoudige lineaire regressie als er meerdere verklarende variabelen zijn). In het algemeen spreken we van een regressiemodel als er een specifieke samenhang bestaat tussen de respons Y en waarnemingen x1 , . . ., xp : Y = fθ (x1 , . . ., xp ) + e waarbij fθ de relatie tussen de waarnemingen x1 , . . ., xp en de respons Y beschrijft, en de stochastische variabele e een niet-waarneembare meetfout is met verwachting nul en onbekende variantie σ 2 . Indien de functie fθ bekend is op de eindig-dimensionale parameter θ na, dan spreken we van een parametrisch model. Het lineaire regressiemodel is hier een voorbeeld van; in dit model is θ = (β0 , . . ., βp ) ∈ Rp+1 en fθ (x1 , . . ., xp ) = β0 + β1 x1 + . . . + βp xp . Het regressiemodel ligt dan vast als waarden voor θ en σ 2 bekend zijn. De functie fθ kan echter ook onbekend zijn op de eindig dimensionale parameter θ en een oneindig dimensionale parameter na. We spreken dan van een semi-parametrisch model. Een voorbeeld van een semi-parametrisch model is het Cox-regressiemodel. In Hoofdstuk 7 komen verschillende regressiemodellen, waaronder het lineaire regressiemodel, uitvoerig aan de orde.

2: Opgaven

9

Opgaven 1. Veronderstel dat aselect n mensen uit een populatie worden gevraagd naar hun politieke voorkeur. Noteer het aantal personen in de steekproef met politieke voorkeur voor partij A met X. De fractie personen in de populatie met politieke voorkeur voor partij A is de onbekende kans p. Beschrijf een bijbehorend statistisch model. Bedenk een intu¨ıtief redelijke “schatting” voor p. 2. Veronderstel dat aselect m + n patiënten met een hoge bloeddruk worden gekozen en geheel willekeurig worden verdeeld in twee groepen ter grootte van m en n. De eerste groep, de “treatment group”, krijgt een bepaald bloeddrukverlagend medicijn toegediend; de tweede groep, de “control group”, ontvangt een placebo. De bloeddruk van iedere patiënt wordt één week na het toedienen van het medicijn of de placebo gemeten. Dit geeft waarnemingen x1 , . . ., xm en y1 , . . ., yn . (i) Formuleer een geschikt statistisch model. (ii) Geef een intu¨ıtief redelijke “schatting” voor het effect van het medicijn op de hoogte van de bloeddruk (meerdere antwoorden zijn mogelijk!). 3. Bij het keuren van een partij goederen gaat men door tot men 3 afgekeurde exemplaren heeft aangetroffen. (i) Formuleer een geschikt statistisch model. (ii) Het derde afgekeurde exemplaar blijkt het 50ste exemplaar te zijn dat men onderzoekt. Geef een schatting van het percentage defecte artikelen in de partij. Beargumenteer je keuze. 4. Het vermoeden bestaat dat er een lineair verband is tussen het inkomen van een persoon en zijn leeftijd en opleidingsniveau (laag, midden, hoog). (i) Beschrijf een lineair regressiemodel met “inkomen” als afhankelijke variabele en “leeftijd” en “opleiding” als onafhankelijke variabelen. Bedenk goed hoe je de variabele “opleiding” in het model opneemt. (ii) Men wil onderzoeken of het geslacht van een persoon invloed heeft op het inkomen. Pas het lineaire regressiemodel aan, zodat dit onderzocht kan worden. 5. Men wil een schatting maken van de gemiddelde lengte van wolvezels in een grote bak. Hiertoe wordt de bak eerst goed geschud, waarna met gesloten ogen een tevoren vastgesteld aantal vezels één voor één uit de bak wordt genomen. Men schat de gemiddelde lengte van de vezels in de bak met de gemiddelde lengte van de wolvezels in de steekproef. Is de geschatte lengte systematisch te groot, systematisch te klein of juist goed?

3 Verdelingsonderzoek

3.1

Introductie

Een statistisch model is een uitdrukking van onze a priori kennis van het kansexperiment waaruit de waargenomen data is voortgekomen. Het model postuleert dat de waarneming X is gegenereerd volgens één van de kansmaten in het model. Hoe vinden we een goed model? In sommige gevallen is het model duidelijk uit de manier waarop het kansexperiment is opgezet. Als bij een opininiepeiling de steekproef inderdaad aselect en zonder teruglegging uit een goed omschreven populatie wordt genomen, dan is de hypergeometrische verdeling onvermijdelijk. Betreffen de waarnemingen aantallen uitgezonden radio-actieve deeltjes, dan is de Poisson-verdeling de juiste keus vanwege de natuurkundige theorie van radioactiviteit. Het is ook mogelijk dat het uitgevoerde experiment sterk lijkt op eerdere experimenten, en dat een bepaald model wordt gesuggereerd door de ervaring in het verleden. Lang niet altijd is een bepaald statistisch model echter geheel onomstreden. Het is dan op z’n minst nodig om het gekozen model te valideren. Soms vinden controles plaats na het schatten van de parameters van het model. Een aantal eenvoudige controles kan ook vooraf worden uitgevoerd. In dit hoofdstuk bespreken we enkele grafische technieken om univariate en multivariate steekproeven te onderzoeken. Deze technieken worden, naast op de data zelf, ook veelvuldig toegepast op “residuen” na het fitten van, bijvoorbeeld, een regressiemodel.

3.2

Univariate Steekproeven

Veronderstel dat de getallen x1 , . . ., xn de resultaten zijn van een herhaaldelijk uitgevoerd experiment. Uit de manier waarop de n experimenten zijn uitgevoerd (steeds vanuit dezelfde beginsituatie, zonder “herinnering” van de voorgaande experimenten) leiden we af dat het redelijk is de n getallen op te vatten als realisaties van onderling onafhankelijke, identiek verdeelde stochastische grootheden X1 , . . ., Xn . Dit legt het statistische model al voor een belangrijk deel vast. De overgebleven vraag is: welke (marginale) verdeling gebruiken we? 3.2.1

Histogrammen

Een kansverdeling beschrijft de “verdeling” van de totale kansmassa 1 over de verschillende mogelijke waarden x. We kunnen een kansverdeling vastleggen door zowel de bijbehorende verdelingsfunctie als de bijbehorende kansdichtheid. Een kansdichtheid is een ingewikkelder object dan de verdelingsfunctie, maar geeft een betere visuele indruk van de verdeling van kansen: de verdeling legt veel kansmassa in punten x waar de waarde van de kansdichtheid f (x) groot is, en weinig in x voor welke f (x) ≈ 0.

3.2: Univariate Steekproeven

11

Een eenvoudige techniek om een indruk te krijgen van een kansdichtheid waaruit data x1 , . . ., xn afkomstig zijn is het histogram. Voor een gegeven partitie a0 < a1 < · · · < am die het bereik van de data x1 , . . ., xn overdekt is dit de functie die op het interval (aj−1 , aj ] een waarde aanneemt die gelijk is aan het aantal datapunten xi die in het interval valt, gedeeld door de lengte van het interval. Als de lengten van alle intervallen (aj−1 , aj ] gelijk zijn, dan wordt het histogram ook wel gedefinieerd zonder door de intervallengten te delen. In dat geval zijn de hoogten van de staven van het histogram gelijk aan de totale aantallen waarnemingen in de verschillende intervallen. De keuze van de intervallen is een kwestie van smaak. Als de intervallen te smal gekozen worden, dan is het histogram over het algemeen te piekerig om kenmerken van de ware kansdicht op te merken. Als de intervallen te breed gekozen worden, gaat daarentegen elk detail verloren en is er nog maar weinig te zeggen over de ware kansdichtheid op basis van het histogram. Om een indruk te krijgen uit welke kansdichtheid data afkomstig zouden kunnen zijn, is het handig het histogram en mogelijke kansdichtheden in één plaatje weer te geven. Dit kan door het histogram te schalen met 1/n, waarbij n het totaal aantal datapunten is. De oppervlakte onder het histogram is dan gelijk aan 1, net zoals dat het geval is bij een kansdichtheid. In x ∈ (aj−1 , aj ] is het geschaalde histogram gelijk aan n X #(1 ≤ i ≤ n: xi ∈ (aj−1 , aj ] 1 = hn (x) = 1a <x ≤a , n(aj − aj−1 ) n(aj − aj−1 ) i=1 j−1 i j

waarbij de indicatorfunctie 1aj−1 <xi ≤aj gelijk is aan 1 als aj−1 < xi ≤ aj en 0 als dit niet het geval is. Een alternatieve schrijfwijze voor deze indicatorfunctie is 1(aj−1 ,aj ] (xi ). Een histogram geeft een goede indruk van de dichtheid waaruit de data x1 , . . ., xn afkomstig zijn, mits de partitie a0 < a1 < · · · < am geschikt gekozen is en het aantal datapunten n niet te klein is. Om dit in te zien beschouwen we x1 , . . ., xn als realisaties van de stochastische variabelen met een dichtheid f en berekenen we de verwachte waarde van het geschaalde histogram hn in termen van X1 , . . ., Xn in een willekeurig punt x waar f (x) > 0. Veronderstel dat voor zekere 1 < j ≤ m geldt dat aj−1 < x ≤ aj dan is deze verwachte waarde gelijk aan n X 1 1 1a <X ≤a = E1aj−1 <X1 ≤aj n(aj − aj−1 ) i=1 j−1 i j aj − aj−1 R aj f (s) ds 1 a = P(aj−1 < X1 ≤ aj ) = j−1 . aj − aj−1 aj − aj−1

Ehn (x) = E

Als f niet te veel varieert over het interval (aj−1 , aj ], dan is de uitdrukking aan de rechterkant ongeveer gelijk aan de waarde van f in dit interval. De berekening leert dat de verwachte waarde van hn (x) bij benadering gelijk is aan f (x). Vanwege de Wet van de Grote Aantallen hebben we bovendien dat de waarde hn (x) in kans naar deze verwachte waarde convergeert. Een histogram geeft dus een indruk van de kansverdeling waaruit een steekproef is gegenereerd. Helaas wordt een goede indruk pas verkregen als een voldoend grote steekproef beschikbaar is (bijvoorbeeld n = 100 of nog liever n = 500). We mogen daarom niet meer dan een eerste indruk van een histogram verwachten. Andere, meer gecompliceerde technieken, kunnen betere resultaten geven. Voorbeeld 3.1. In Figuur 3.1 zijn histogrammen getekend van de lichaamslengte (in cm) van 100 mannen (links) en 110 vrouwen (rechts).‡ De histogrammen zijn zo geschaald dat de oppervlaktes onder de histogrammen gelijk aan 1 zijn. In beide figuren is eveneens de dichtheid van een normale verdeling getekend. De verwachting en variantie van deze normale verdelingen zijn gelijk aan het steekproefgemiddelde en de steekproefvariantie van de bijbehorende data (zie Hoofdstuk 4).

‡ Bron: De data zijn verzameld door de afdeling Biologische Psychologie van de Vrije Universiteit in het kader van een onderzoek naar gezondheid, levensstijl en persoonlijkheid.

3: Verdelingsonderzoek

0.00

0.00

0.02

0.02

0.04

0.04

0.06

0.06

12

165

175

185

195

155

165

175

185

Figuur 3.1. Histogram van de lichaamslengte van 100 mannen (links) en 110 vrouwen (rechts), tezamen met de kansdichtheden van de normale verdeling met de verwachtingen gelijk aan de steekproefgemiddelden en de varianties gelijk aan de steekproefvarianties van de data.

0.4 0.3 0.2 0.1 0.0

0.0

0.1

0.2

0.3

0.4

Voorbeeld 3.2 (Normale verdeling). Figuur 3.2 geeft de dichtheid van de standaard normale verdeling tezamen met vier realisaties van het histogram, gebaseerd op 30, 30, 100 en 100 waarnemingen, waarbij de partities gekozen werden door het statistische softwarepakket R. De figuren linksboven en rechtsonder vertonen duidelijke afwijkingen van symmetrie. Omdat de data uit de normale verdeling werden gegenereerd is dit slechts te wijten aan toevalsvariatie.

-1

0

1

2

3

-3

-2

-1

0

1

2

3

-3

-2

-1

0

1

2

3

-3

-2

-1

0

1

2

3

0.3 0.2 0.1 0.0

0.0

0.1

0.2

0.3

0.4

-2

0.4

-3

Figuur 3.2. Histogrammen van steekproeven van 30, 30, 100 en 100 (boven, boven, onder, onder) waarnemingen uit de standaard normale dichtheid en de ware dichtheid.

3.2.2

Boxplots

Een boxplot is een grafische weergave van de data die een indruk geeft van de locatie en de spreiding van de data, eventuele extreme waarden in de waarnemingen en de symmetrie van de verdeling waaruit de waarnemingen afkomstig zijn. In de boxplot staan de waarnemingen uitgezet langs de verticale as. De onderkant van de “box” staat getekend ter hoogte van het

3.3: Samenhang

13

kleinste kwartiel, en de bovenkant ter hoogte van het grootste kwartiel van de data. Het kleinste (respectievelijk grootste) kwartiel van de data is die waarde x zodanig dat een kwart van de waargenomen data kleiner (respectievelijk groter) is dan x. De breedte van de box is willekeurig. In de box staat ter hoogte van de mediaan van de data een horizontale lijn. De mediaan is de middelste waarde in de rij gesorteerde waarnemingen. Aan de boven- en onderkant van de box staan zogeheten whiskers getekend. De whisker aan de bovenkant verbindt de box met de grootste waarneming die binnen 1.5 maal de interkwartiel afstand boven het grootste kwartiel ligt. De interkwartiel afstand is de afstand tussen het bovenste en het onderste kwartiel, ofwel de hoogte van de box. De whisker aan de onderkant wordt op analoge wijze getekend. Waarnemingen die buiten de whiskers vallen worden apart aangegeven, bijvoorbeeld met een sterretje, rondje of streepje.

−6

−4

−2

0

2

4

Voorbeeld 3.3. In Figuur 3.3 staan boxplots getekend van steekproeven uit de exponentiële verdeling met parameter 1, de standaard normale verdeling en de standaard Cauchy-verdeling. De steekproeven uit de exponentiële en de Cauchy-verdeling bevatten extreme waarden, te zien aan de rondjes buiten de whiskers. De boxplot in het midden geeft aan dat de data uit de standaard normale verdeling aardig symmetrisch rond de mediaan liggen en geen extreme waarden bevatten.

exp(1)

N(0,1)

Cauchy

Figuur 3.3. Boxplots van steekproeven ter grootte 20 uit de standaard exponenti¨ ele verdeling (links), de standaard normale verdeling (midden) en de standaard Cauchy-verdeling (rechts).

3.3

Samenhang

In veel gevallen zijn de waarnemingen xi geen getallen, maar vectoren xi = (xi,1 , . . ., xi,d ). We zijn dan vaak ge¨ınteresseerd in de relatie tussen de verschillende co¨ ordinaten. We zullen ons in deze paragraaf beperken tot vectoren met twee co¨ ordinaten en noteren deze met (x i , yi ) (in plaats van (xi,1 , xi,2 )). Een scatterplot van een steekproef van tweedimensionale data (x1 , y1 ), . . ., (xn , yn ) is een plot van deze punten in het platte vlak. Is er een duidelijk verband tussen de x- en y-co¨ ordinaten van de data, dan is dit op het oog onmiddellijk zichtbaar. De variabelen in het rechterplaatje van Figuur 3.4 vertonen bijvoorbeeld een duidelijk lineair verband, terwijl in het linkerplaatje geen samenhang is te ontdekken.

14

-1.5

-1

-1.0

-0.5

0

0.0

1

0.5

1.0

3: Verdelingsonderzoek

-3

-2

-1

0

1

2

3

-3

-2

-1

0

1

2

3

Figuur 3.4. Scatterplots van twee steekproeven van 50 punten. Links met onafhankelijke co¨ ordinaten (r x,y = −0.05), rechts met co¨ ordinaten die een lineair verband vertonen (rx,y = 0.87).

Het lineaire verband in het rechterplaatje van Figuur 3.4 is onmiskenbaar, maar niet perfect. De punten liggen niet exact op een rechte lijn, maar variëren rond een (denkbeeldige) rechte. Notatie 3.4. Het steekproefgemiddelde van een steekproef X1 , . . ., Xn is de stochastische grootheid n 1X X= Xi . n i=1 De steekproefvariantie van een steekproef X1 , . . ., Xn is de stochastische grootheid n

2 SX

1 X = (Xi − X)2 . n−1 i=1

De steekproefcorrelatiecoëfficiënt van een steekproef van paren (X 1 , Y1 ), . . ., (Xn , Yn ) is Pn (Xi − X)(Yi − Y ) p p rX,Y = i=1 . 2 (n − 1) SX SY2

De steekproefcorrelatiecoëfficiënt rx,y van de waargenomen paren (x1 , y1 ), . . ., (xn , yn ) is een getalsmaat voor de kracht van het lineaire verband en ligt tussen −1 en 1. De waarde kan als volgt worden ge¨ınterpreteerd: (i) Als rx,y = 1, dan liggen de n punten in de scatterplot precies op de lijn y = y +(sy /sx ) (x− x) (perfect positief verband). (ii) Als rx,y = −1, dan liggen de n punten in de scatterplot precies op de lijn y = y − (sy /sx ) (x − x) (perfect negatief verband). (iii) Zijn X1 , . . ., Xn en Y1 , . . ., Yn onafhankelijke steekproeven, dan zal de gerealiseerde rx,y waarden dicht bij 0 aannemen. De eerste twee beweringen en dat |rx,y | ≤ 1 zijn een gevolg van de ongelijkheid van CauchySchwarz uit de lineaire algebra.[ De derde bewering is een gevolg van het feit dat onafhankelijke stochastische grootheden ongecorreleerd zijn, gecombineerd met het intu¨ıtief aannemelijke feit dat de steekproefcorrelatiecoëfficiënt de populatiecorrelatiecoëfficiënt ρ= √

cov(X, Y ) E(X − EX)(Y − EY ) √ p =p var X var Y E(X − EX)2 E(Y − EY )2

zal benaderen voor n groot. Omdat cov(X, Y ) = E(X − EX)(Y − EY ) = E(XY ) − EXEY is ρ gelijk aan 0 voor onafhankelijk stochastische grootheden X en Y : onafhankelijke stochastische [

Het inwendig product van vectoren a en b in Rn voldoet aan |ha, bi| ≤ kak kbk voor k · k de Euclidische norm.

3: Opgaven

15

0

-2

2

-1

0

4

1

6

2

8

3

grootheden zijn ongecorreleerd. Een verdere interpretatie van de steekproefcorrelatiecoëfficiënt wordt gegeven bij de behandeling van het lineaire regressiemodel in Hoofdstuk 7. We mogen bewering (iii) niet omdraaien in de zin dat een correlatie dicht bij 0 zou impliceren dat de twee co¨ ordinaten onafhankelijk zijn. Dit wordt ge¨ıllustreerd in Figuur 3.5. In het linkerplaatje is een duidelijk lineair verband waarneembaar, corresponderend met een correlatiecoëfficiënt van 0.98. Het rechterplaatje is een scatterplot van de punten (x i , yi2 ) voor de punten (xi , yi ) uit het linkerplaatje. Het kwadratische verband is duidelijk zichtbaar. De “sterkte van het verband” tussen de twee co¨ ordinaten in het rechterplaatje doet niet onder voor de sterkte in het linkerplaatje. De steekproefcorrelatiecoëfficiënt voor de punten in het rechterplaatje is echter gelijk aan −0.05. Blijkbaar is deze getalsmaat blind voor het aanwezige kwadratische verband.

-2

-1

0

1

2

3

-2

-1

0

1

2

3

Figuur 3.5. Scatterplots van twee steekproeven van 50 punten, met steekproefcorrelatieco¨ effici¨ enten, respectievelijk, 0.98 en -0.05. Het rechterplaatje geeft de punten (xi , yi2 ) voor de punten (xi , yi ) uit het linkerplaatje.

Opgaven 1. Veronderstel dat hn het geschaalde histogram van een steekproef X1 , . . ., Xn uit een dichtheid f is. De partitie van het histogram wordt gegeven R a door a0 < a1 < . . . < am . Bewijs dat voor aj−1 < x ≤ aj geldt dat hn (x) → (aj − aj−1 )−1 a j f (s) ds met kans 1, als n → ∞. j−1

2. Zij X een standaard normaal verdeelde stochastische grootheid. Bereken de correlatiecoëfficiënt tussen de stochastische grootheden X en Y = X 2 .

3. Veronderstel dat X en Y onderling onafhankelijk zijn en beide standaard normaal verdeeld. Bereken de correlatiecoëfficiënt tussen X en Z waar Z = X + Y .

4 Schatters

4.1

Introductie

Een statistisch model bestaat uit alle kansverdelingen welke a priori mogelijk worden geacht voor de gegeven data. Gegeven een correct opgesteld model gaan we ervan uit dat de data volgens één van de kansverdelingen in het model is gegenereerd. Na het opstellen van een geschikt statistisch model is de volgende stap het bepalen welke kansverdeling binnen het model het best aansluit bij de gegevens. Als het model wordt gegeven door een parameter, dan is dit equivalent met het bepalen van de best passende parameterwaarde, vaak aangeduid als de “ware” parameter. In de statistiek heet dit proces “schatten”. Andere namen zijn “fitten” en “leren”. Veronderstel dat de kansverdeling van X afhangt van een onbekende parameter θ, zodat het statistische model de vorm {Pθ : θ ∈ Θ} heeft, voor Pθ de kansverdeling van X als θ de “ware” parameterwaarde is. Op grond van een waarneming x willen we de ware waarde van θ schatten, of wellicht de waarde van een functie g(θ) van θ, bijvoorbeeld de eerste co¨ ordinaat θ1 als θ = (θ1 , θ2 ). “Schatten” betekent hier het doen van een uitspraak over θ of g(θ) van de vorm: “ik denk dat g(θ) bij benadering gelijk is aan T (x)”, voor zekere waarde T (x) die van de waargenomen waarde x afhangt. Definitie 4.1. Een schatter (Engels: estimator) of statistiek (Engels: statistic) is een stochastische vector T (X) die alleen van de waarneming X afhangt. De bijbehorende schatting (Engels: estimate), bij gerealiseerde waarneming x, is T (x). Volgens deze definitie zijn heel veel objecten schatters. Waar het om gaat is dat T (X) een functie van X is die niet van de parameter θ mag afhangen: we moeten T (x) kunnen uitrekenen op grond van de data x. Na verrichting van de waarneming krijgt T een gerealiseerde waarde t = T (x), waarmee we θ (of g(θ)) schatten. We korten T (X) heel vaak af tot T . Hoewel iedere functie van de waarneming een schatter is, is niet iedere schatter een goede schatter. Een goede schatter voor g(θ) is een functie T van de waarneming zodanig dat T “dichtbij” de te schatten waarde g(θ) ligt. Een maat die wiskundig relatief eenvoudig is te hanteren is de verwachte kwadratische fout (Engels: mean square error of MSE). Voor een schatter T voor de waarde g(θ) wordt deze gedefinieerd als

2 MSE(θ; T ) = Eθ T − g(θ) .

We geven de voorkeur aan een schatter met een kleine verwachte kwadratische fout (MSE) voor alle parameterwaarden van θ tegelijk. De verwachte kwadratische fout van een reëelwaardige schatter T kan worden ontbonden in twee termen: 2 MSE(θ; T ) = varθ T + Eθ T − g(θ)

4.2: Maximum Likelihood-Schatters

17

(ga na). Beide termen in deze decompositie zijn niet-negatief. Dus de verwachte kwadratische fout kan alleen klein zijn als beide termen klein zijn. Als de tweede term gelijk aan 0 is, dan heet de schatter zuiver. Definitie 4.2. Een schatter T heet zuiver (Engels: unbiased) voor het schatten van g(θ) als Eθ T = g(θ) voor alle θ ∈ Θ. De onzuiverheid (Engels: bias) is gedefinieerd als E θ T − g(θ). ˆ Het dakje geeft Zowel schatters als schattingen van θ worden vaak aangegeven met θ. aan dat θˆ een functie van de waarneming is, maar deze notatie maakt geen verschil tussen de ˆ ˆ stochastische vector of zijn realisatie: θˆ kan zowel θ(X) als θ(x) betekenen.

4.2

Maximum Likelihood-Schatters

De “methode van de maximum likelihood-schatters” (Nederlands: methode van de meest aannemelijke schatters) is de meest gebruikte methode om schatters voor een onbekende parameter te vinden. Voordat deze methode in het algemeen wordt gepresenteerd, wordt voor het (eenvoudige) geval van de binomiale verdeling de maximum likelihood-schatter afgeleid in het volgende voorbeeld. Voorbeeld 4.3 (Binomiale verdeling). Veronderstel dat we 10 keer met een onzuivere munt gooien. De kans p op “kop” is bij deze munt niet noodzakelijkerwijze 1/2. Definieer X als het aantal malen “kop” in de 10 worpen. De stochastische variabele X heeft dan een binomiale verdeling met parameters 10 en onbekende p ∈ [0, 1]. Stel dat we 3 maal “kop” werpen. De kans op deze uitkomst is gelijk aan 10 3 Pp (X = 3) = p (1 − p)7 . 3

0.00

0.05

0.10

0.15

0.20

0.25

De kans p is onbekend en moet geschat worden. Welke waarde voor p is nu meest waarschijnlijk?

0.0

0.2

0.4

0.6

0.8

1.0

p Figuur 4.1. De kans Pp (X = 3) als functie van p waar de stochast X binomiaal verdeeld is met parameters 10 en p.

In Figuur 4.1 is de kans Pp (X = 3) getekend als functie van p. We zien dat er precies één waarde voor p is die deze kans maximaliseert, namelijk de waarde 0.3. Deze waarde voor p kent de grootste kans toe aan de waarneming “3 maal kop”. De schatting pˆ = 0.3 blijkt in deze situatie de maximum likelihood-schatting te zijn.

18

4: Schatters

De maximum likelihood-methode vereist de specificatie van de likelihood-functie, welke wordt afgeleid uit de kansdichtheid van de waarneming. Hierbij verstaan we onder een kansdichtheid pθ van een stochastische vector X deR functie x 7→ Pθ (X = x) als X discreet verdeeld is en de functie pθ zodanig dat Pθ (X ∈ B) = B pθ (x) dx als X continu verdeeld is.

Definitie 4.4. Zij X een stochastische vector met een kansdichtheid pθ die van een parameter θ ∈ Θ afhangt. De functie θ 7→ L(θ; x): = pθ (x) opgevat als functie van θ ∈ Θ voor vaste x heet de likelihood-functie (Nederlands: aannemelijkheidsfunctie). Vaak is X = (X1 , . . ., Xn ) een vector met onderling onafhankelijke identiekQverdeelde co¨ ordinaten Xi . Dan is de dichtheid van X in (x1 , . . ., xn ) gelijk aan het product ni=1 pθ (xi ) van marginale dichtheden van de X1 , . . ., Xn . Voor waargenomen waarden (x1 , . . ., xn ) is de likelihood-functie gelijk aan θ 7→ L(θ; x1 , . . ., xn ) =

n Y

pθ (xi ),

i=1

waarin nu pθ de (marginale) dichtheid van een enkele Xi weergeeft. De algemene definitie van maximum likelihood-schatters is echter geldig voor een waarnemingsvector van willekeurige vorm, en we geven er daarom de voorkeur aan de waarneming als x te schrijven, in plaats van (x1 , . . ., xn ), en de likelihood-functie als L(θ; x) ≡ pθ (x) . Definitie 4.5. De maximum likelihood-schatting voor θ is die waarde T (x) ∈ Θ die de functie θ 7→ L(θ; x) maximaliseert. De maximum likelihood-schatter (Nederlands: meest aannemelijke schatter) is de bijbehorende schatter T (X). In het geval van een discrete kansverdeling kan de maximum likelihood-schatting worden omschreven als: die waarde van de parameter die de grootste waarschijnlijkheid toekent aan de waargenomen waarde x. We maximaliseren in dat geval immers de kansdichtheid pθ (x) = Pθ (X = x) naar θ voor vaste x (zie Voorbeeld 4.3). Dit is een intu¨ıtief redelijk schattingsprincipe en verklaart de naam. Dit principe moet echter alleen beschouwd worden als een schattingsmethode: maximum likelihood-schatters zijn niet noodzakelijkerwijze de beste schatters, ondanks de mooie naam. Als g: Θ → H een 1 − 1-duidige functie is met een verzameling H als bereik, dan zouden we het model ook door de parameter η = g(θ) ∈ H kunnen parametriseren in plaats van door ˆ de maximum likelihood-schatter voor η is, als θ ∈ Θ. Het volgt direct uit de definitie dat g(θ) ˆ θ de maximum likelihood-schatter voor θ is. In overeenstemming hiermee definiëren we voor ˆ iedere willekeurige functie g de maximum likelihood-schatter voor g(θ) simpelweg als g( θ). Bij een gegeven model is het uitrekenen van de maximum likelihood-schatter een kwestie van calculus. Vaak geschiedt dit door de likelihood-functie te differentiëren en de afgeleiden gelijk aan nul te stellen. (Het geval van de homogene verdeling in is hier echter een uitzondering op.) Een trucje dat het rekenwerk beperkt (vooral bij onafhankelijke waarnemingen) is om eerst de logaritme van de likelihood te nemen. Omdat de logaritme een monotone functie is, geldt dat de waarde θˆ de functie θ 7→ L(θ; x) maximaliseert dan en slechts dan als deze waarde de functie θ 7→ log L(θ, x) maximaliseert. (Het gaat om de plaats waar het maximum wordt aangenomen, niet de grootte van het maximum!) Voor vaste x wordt de log likelihood-functie gegeven door θ 7→ log L(θ; x) = log pθ (x).

Als L differentieerbaar is in θ ∈ Θ ⊂ Rk en zijn maximum in een inwendig punt van Θ aanneemt, dan geldt ∂ log L(θ; x)|θ=θˆ = 0, j = 1, . . ., k. ∂θj


19

Dit stelsel van likelihood-vergelijkingen is lang niet altijd expliciet oplosbaar. Zonodig gebruikt men iteratietechnieken om stapsgewijs een steeds betere benadering van de oplossing te verkrijgen, uitgaande van een geschikte startwaarde. De vector van partiële afgeleiden (gradiënt) van θ 7→ log L(θ; x) wordt de score-functie van het statistische model genoemd. Als de waarneming X = (X1 , . . ., Xn ) is opgebouwd uit onafhankelijke, identiek verdeelde deelwaarnemingen Xi , dan bezit de likelihood L(θ; x) voor waargenomen x de productstructuur Q L(θ; x) = i pθ (xi ). De log likelihood is dan θ 7→ log L(θ; x1 , . . ., xn ) = log

n Y

pθ (xi ) =

i=1

n X

log pθ (xi ),

i=1

waarin pθ de (marginale) dichtheid van een enkele Xi weergeeft. De afgeleide van log L, de score-functie, is de som van de score-functies voor de individuele waarnemingen. De likelihoodvergelijkingen hebben dan de vorm n X i=1

`˙θ (xi )|θ=θˆ = 0,

met `˙θ (xi ) = ∇θ `θ (xi )

en

`θ (xi ) = log pθ (xi ).

De gradiënt `˙θ is de “score-functie voor één waarneming”. In verschillende voorbeelden wordt het maximum van de likelihood-functie niet in het inwendige van de parameterverzameling aangenomen. Dan is de maximum likelihood-schatter θˆ meestal niet een stationair punt van de afgeleide van de likelihood-functie maar een randmaximum, en gelden de voorgaande vergelijkingen niet. In weer andere voorbeelden is de likelihoodfunctie niet overal differentieerbaar (of zelfs continu), en voldoet de maximum likelihoodschatter evenmin aan de likelihood-vergelijkingen. Voorts is het mogelijk dat de likelihoodfunctie meerdere (locale) maxima en ook minima bezit. Dan kunnen de likelihood-vergelijkingen meer dan één oplossing bezitten. De maximum likelihood-schatter is per definitie het globale maximum van de likelihood-functie. In Definitie 4.5 wordt de maximum likelihood-schatter gebaseerd op de maximum likelihood-schatting. In de praktijk schrijft men echter vaak direct de (log) likelihood-functie in termen van de stochastische grootheid X in plaats van de realisatie x en leidt op die manier direct de schatter af door deze functie te maximaliseren naar θ. Deze verkorte notatie wordt gehanteerd in de volgende voorbeelden van toepassingen van de maximum likelihood-methode. Voorbeelden waarin de methode wordt toegepast op regressiemodellen zijn te vinden in Hoofdstuk 7. Voorbeeld 4.6 (Alternatieve verdeling). De kansdichtheid van de alternatieve verdeling alt(p) kan worden geschreven als x 7→ px (1 − p)1−x ; voor x = 0 staat hier 1 − p en voor x = 1 staat er p. Voor een steekproef X1 , . . ., Xn uit de alt(p)-verdeling is de log likelihood-functie derhalve p 7→ log L(p; X1 , . . ., Xn ) = log

n Y

i=1

pXi (1 − p)1−Xi =

n X i=1

n X Xi log p + n − Xi log(1 − p). i=1

Pn Neem de parameterverzameling gelijk aan het interval [0, 1]. Als 0 < i=1 Xi < n, dan geldt log L(p; X) → −∞ als p ↓ 0 of p ↑ 1, zodat L(p; X) zijn maximum aanneemt op (0, 1). Nulstellen van de afgeleide naar p geeft één oplossing; de maximum likelihood-schatter pˆ = X. P Als ni=1 Xi gelijk is aan 0 of n, dan heeft L(p; X) een randmaximum in 0 of 1. Ook in deze gevallen kan de maximum likelihood-schatter worden geschreven als pˆ = X.

4: Schatters

0

2*10^-6

6*10^-6

10^-5

20

0.0

0.2

0.4

0.6

0.8

1.0

Figuur 4.2. Een realisatie van de likelihood-functie als functie van p voor een steekproef uit de alternatieve verdeling. De waargenomen waarde is Σn ˆ = 0.25. i=1 xi = 5 voor n = 20, en de maximum likelihood-schatting is p

Voorbeeld 4.7 (Exponenti¨ ele verdeling). Veronderstel dat X1 , . . ., Xn een steekproef is uit de exponentiële verdeling met onbekende parameter λ > 0. Dan is de log likelihood-functie voor X1 , . . ., Xn gelijk aan λ 7→ log

n Y

i=1

λe−λXi = n log λ − λ

n X

Xi .

i=1

De parameterruimte voor λ is (0, ∞). Nulstellen van de afgeleide van de log likelihood-functie ˆ = 1/X. De tweede afgeleide van naar λ en de gevonden vergelijking oplossen naar λ geeft λ ˆ heeft de likelihood-functie de log likelihood-functie naar λ is negatief voor alle λ > 0, dus in λ ook daadwerkelijk een maximum. De maximum likelihood-schatter van Eθ Xi kunnen we hieruit afleiden. Definieer de functie g als g(λ) = 1/λ voor λ > 0. Dan geldt EXi = g(λ). De maximum likelihood-schatter voor ˆ = 1/λ ˆ = X. EXi = 1/λ = g(λ) is daarom gelijk aan g(λ) Voorbeeld 4.8 (Normale verdeling). De log likelihood-functie voor een steekproef X 1 , . . ., Xn uit de N (µ, σ 2 )-verdeling wordt gegeven door (µ, σ 2 ) 7→ log

n Y

i=1

√

1 2πσ 2

1

e− 2 (Xi −µ)

2

/σ 2

= − 21 n log 2π − 21 n log σ 2 −

n 1 X (Xi − µ)2 . 2σ 2 i=1

We nemen de natuurlijke parameterruimte voor de parameter θ = (µ, σ 2 ): Θ = R × (0, ∞). De partiële afgeleiden van de log likelihood naar µ en σ 2 zijn n 1X ∂ log L(µ, σ 2 ; X) = 2 (Xi − µ) ∂µ σ i=1

n n 1 X ∂ 2 (Xi − µ)2 . log L(µ, σ ; X) = − + ∂σ 2 2σ 2 2σ 4 i=1

Nulstellen van de eerste vergelijking geeft één oplossing: µ ˆ = X. In deze waarde voor µ heeft de log likelihood inderdaad een globaal maximum voor iedere σ 2 > 0 aangezien de waarde van de log likelihood naar −∞ gaat voor µ → ±∞. Vervolgens substitueren we µ = µ ˆ in de tweede partiële afgeleide, stellen deze gelijk aan 0 en lossen de likelihood-vergelijking vervolgens Pn op naar σ 2 . Dit geeft opnieuw één oplossing: σ ˆ 2 = n−1 i=1 (Xi − X)2 . Om gelijke reden als hiervoor heeft de log likelihood functie in deze waarde inderdaad een maximum. (Overigens levert het maximaliseren van de log likelihood-functie naar σ in plaats van σ 2 de wortel uit σ ˆ 2 als maximum likelihood-schatter voor σ op.) Om te controleren of de (differentieerbare) log likelihood-functie een maximum heeft in de gevonden oplossing van de likelihood-vergelijkingen,


21

kan ook de Hessiaan-matrix van de log likelihood-functie in het punt (ˆ µ, σ ˆ 2 ) berekend worden, welke hier gelijk is aan 1 −nˆ σ2 0 . 0 −n/2 σ ˆ4

Beide eigenwaarden van deze matrix zijn negatief en daarmee heeft de log likelihood een maximum in het punt (ˆ µ, σ ˆ 2 ). De gevonden maximum likelihood-schatter voor (µ, σ 2 ) is gelijk aan n n−1 1X 2 X, (Xi − X)2 = X, SX n i=1 n

met

2 SX =

n 1 X (Xi − X)2 . n − 1 i=1

Het steekproefgemiddelde is zuiver voor µ, maar de maximum likelihood-schatter σ ˆ 2 heeft een lichte onzuiverheid. Vanwege de kleine onzuiverheid wordt vaak de voorkeur gegeven aan de 2 2 steekproefvariantie SX = (n/(n−1))ˆ σ 2 . De verwachte kwadratische 2 fout van SX is echter groter 2 dan die van σ ˆ , en beide verliezen het van (n − 1)/(n + 1) SX in termen van de verwachte kwadratische fout.] Omdat het verschil klein is voor grote aantallen waarnemingen, maakt het meestal niet veel uit welke van deze schatters wordt gebruikt. Voorbeeld 4.9 (Normale verdeling met restrictie). Veronderstel dat de waarnemingen X1 , . . ., Xn onafhankelijk en normaal verdeeld zijn met verwachting µ en variantie 1, waarbij bekend is dat µ ≥ 0. Voor x1 , . . ., xn een realisatie van X1 , . . ., Xn , neemt de likelihood-functie op R een absoluut maximum aan in x. Omdat x negatief kan zijn en bekend is dat µ ≥ 0, is x niet de maximum likelihood-schatting. In het geval dat x ≤ 0, neemt de likelihood-functie op de parameterverzameling [0, ∞) een randmaximum aan in 0. De maximum likelihood-schatting is x als deze niet-negatief is en anders 0. De bijbehorende maximum likelihood-schatter is dan X1X≥0 ; X als X ≥ 0 en 0 anders. Een statistisch model en de maximum likelihood-schatter worden bepaald door zowel de vorm van de dichtheid van de waarneming als de definitie van de parameterverzameling! Voorbeeld 4.10 (Gamma-verdeling). Stel dat X1 , . . ., Xn een steekproef is uit de Gammaverdeling met kansdichtheid xα−1 λα e−λx pα,λ (x) = . Γ(α) Hierin zijn α > 0 en λ > 0 de onbekende vormen inverse schaalparameter, en Γ de Gammafunctie Z ∞ Γ(α) = sα−1 e−s ds. 0

De log likelihood-functie voor X1 , . . ., Xn is dan gelijk aan (α, λ) 7→ log

n Y X α−1 λα e−λXi i

i=1

Γ(α)

= (α − 1)

n X i=1

log Xi + nα log λ − λ

n X i=1

Xi − n log Γ(α).

De parameterruimte voor θ = (α, λ) nemen we gelijk aan Θ = [0, ∞) × [0, ∞). Om de maximum likelihood-schatters voor α en λ te vinden, bepalen we de partiële afgeleiden van de log likelihood-functie naar λ en α n

∂ nα X log L(α, λ; X1 , . . ., Xn ) = − Xi , ∂λ λ i=1 n

X ∂ log L(α, λ; X1 , . . ., Xn ) = log Xi + n log λ − n ∂α i=1

R ∞ α−1 s log s e−s ds 0 R . ∞ α−1 −s s e ds 0

] Het vereist enig rekenwerk om deze bewering te staven. Stelling 5.15 kan worden gebruikt om dit rekenwerk te vergemakkelijken.

22

4: Schatters

(In de afgeleide naar α hebben we de Gamma-functie α 7→ Γ(α) onder het integraalteken gedifferentieerd en gebruikt dat (∂/∂α)sα = sα log s.) De partiële afgeleiden zijn gelijk aan 0 in ˆ dit geeft twee likelihood-vergelijkingen. Uit de eerste de maximum likelihood-schatters (ˆ α, λ); ˆ vergelijking volgt onmiddellijk dat λ = α ˆ /X. Dit substitueren we in de tweede likelihoodvergelijking. Dit geeft R ∞ α−1 n X s ˆ log s e−s ds log Xi + n log α ˆ − n log X − n 0 R ∞ α−1 = 0. s ˆ e−s ds 0 i=1

Deze vergelijking heeft geen expliciete oplossing voor α ˆ , maar kan numeriek, met een iteratieve methode, worden opgelost wanneer een realisatie van X1 , . . ., Xn is waargenomen. Voor de meeste numerieke algoritmen zijn startwaarden nodig van waaruit gezocht wordt naar een oplossing van de vergelijking. De momentenschattingen kunnen als startwaarden dienen (zie Paragraaf 4.3). ˆ te ˆ =α De gevonden waarde α ˆ substitueren we vervolgens in de vergelijking λ ˆ /X om λ vinden. Om te controleren of de log likelihood-functie in de oplossing ook daadwerkelijk een ˆ berekenen. Als maximum heeft, moeten we de eigenwaarden van de Hessiaan-matrix in (ˆ α, λ) ˆ dan is (ˆ ˆ inderdaad de maximum likelihooddeze beide eigenwaarden negatief zijn in (ˆ α, λ), α, λ) schatter voor (α, λ).

4.3

Momentenschatters

De methode van de momenten is een alternatief voor de maximum likelihood-methode. Omdat de momentenmethode vaak niet de volledige informatie uit het statistische model gebruikt, zijn momentenschatters vaak minder efficiënt dan maximum likelihood-schatters. Aan de andere kant is de methode soms makkelijker te implementeren. Bovendien vereist de methode alleen de theoretische vorm van de momenten en niet de gehele kansverdeling van de waarnemingen. Aangezien deze momenten vaak gemakkelijker op een realistische manier zijn te modelleren dan de hele kansverdeling, kan dit een groot voordeel zijn. Het gebruik van een verkeerd model om schatters te construeren kan daardoor worden voorkomen. Veronderstel dat X een stochastische variabele is met een verdeling die bekend is op een parameter θ na. Het j e moment van X is gedefinieerd als Eθ (X j ), mits deze verwachting bestaat. Op basis van onderling onafhankelijke en identiek verdeelde variabelen X1 , . . ., Xn uitP dezelfde n verdeling kan het j e moment geschat worden met het j e steekproefmoment: X j = n−1 i=1 Xij . Dat dit een goede schatter is voor Eθ (X j ) volgt uit de Wet van de Grote Aantallen. De momentenschatter voor θ is die waarde θˆ waarvoor het j e moment overeenkomt met het j e steekproefmoment: Eθˆ(X j ) = X j . In de praktijk geven we de voorkeur aan de momentenschatter die gevonden wordt door j zo klein mogelijk te nemen. Voor een eendimensionale parameter θ volstaat j = 1, mits de verwachtingswaarde van de marginale verdeling afhangt van θ. Wanneer het eerste moment niet afhangt van θ, wordt j = 2 gekozen, etc. Indien θ meerdimensionaal is, zijn er meerdere vergelijkingen nodig om een unieke oplossing voor θˆ te krijgen. In dat geval wordt de momentenschatter θˆ opgelost uit de vergelijkingen voor j = 1, . . ., k met k het kleinste gehele getal waarvoor het stelsel vergelijkingen een unieke oplossing bezit. De momentenschatter voor ˆ met θˆ de momentenschatter g(θ) met g: Θ → H een functie met bereik H is gedefinieerd als g(θ) voor θ. Voorbeeld 4.11 (Exponenti¨ ele verdeling). Veronderstel dat X1 , . . ., Xn een steekproef is uit een exponentiële verdeling met onbekende parameter λ. Dan is Eλ Xi = 1/λ. De momentenˆ op te lossen naar λ. ˆ Dit geeft schatter van λ wordt nu gevonden door de vergelijking X = 1/λ ˆ λ = 1/X als momentenschatter voor λ. Deze schatter is ook de maximum likelihood-schatter voor λ (zie Voorbeeld 4.7).

4: Opgaven

23

Voorbeeld 4.12 (Gamma-verdeling). Veronderstel dat X1 , . . ., Xn de gamma-verdeling met onbekende vormen inverse schaalparameter α en λ hebben. Dan is Eα,λ Xi = α/λ en varα,λ Xi = α/λ2 , en dus is het tweede moment gelijk aan Eα,λ Xi2 = varα,λ Xi + (Eα,λ Xi )2 = α(1 + α)/λ2 . De momentenschatters voor α en λ worden gevonden uit het oplossen van de volgende twee vergelijkingen ˆ=X Eα, ˆ /λ ˆ Xi = α ˆλ 2 ˆ2 = X 2 Eα, ˆ (1 + α ˆ )/λ ˆ Xi = α ˆλ

ˆ Dit geeft naar α ˆ en λ. α ˆ=

(X)2 X2

−

(X)2

ˆ= λ

en

X X2

− (X)2

.

Opgaven 1. Bereken de maximum likelihood-schatter voor θ gebaseerd op een steekproef X1 , . . ., Xn uit de Poisson(θ)-verdeling. 2. Zij X1 , . . ., Xn een steekproef uit een Weibull-verdeling, waarvan de kansdichtheid wordt gegeven door a pθ (x) = θaxa−1 e−θx , voor x > 0 en 0 anders. Hierin is a een bekend getal en θ > 0 is een onbekende parameter. (i) Bepaal de maximum likelihood-schatter voor θ. (ii) Bepaal de maximum likelihood-schatter voor 1/θ. 3. Zij X1 , . . ., Xn een steekproef uit een verdeling met kansdichtheid pθ (x) = θxθ−1 ,

voor x ∈ (0, 1)

en 0 anders. Hierin is θ > 0 een onbekende parameter. (i) Bereken µ = g(θ) = Eθ X1 . (ii) Bepaal de maximum likelihood-schatter voor µ. 4. Zij X1 , . . ., Xn een steekproef uit een kansverdeling met dichtheid pθ (x) = θ(1 + x)−(1+θ) ,

voor x ≥ 0

en 0 anders, waarin θ > 0 onbekend. Bepaal de maximum likelihood-schatter voor θ. 5. Zij X1 , . . ., Xn een steekproef uit de binomiale verdeling met parameters n en p, waarbij p ∈ [0, 1] onbekend is. Bepaal de maximum likelihood-schatter en de momentenschatter voor p. 6. Zij X1 , . . ., Xn een steekproef uit een kansverdeling met dichtheid pθ (x) = θ(1 + x)−(1+θ)

voor x > 0,

en 0 anders, met θ > 1 onbekend. Bepaal de momentenschatter voor θ.

5 Toetsen

Bij wetenschappelijk onderzoek, in de industrie en in het dagelijks leven is het vaak gewenst na te gaan of bepaalde vragen al dan niet bevestigend beantwoord kunnen worden. Helpt een bepaalde therapie? Speelt leeftijd of sekse van de patiënt hierbij een rol? Is het ene type auto veiliger dan het andere? Bevat een partij excessief veel defecte artikelen? Gaat het ene type lamp langer mee dan het andere? Komt het DNA-profiel van de verdachte overeen met het DNAprofiel dat gevonden is op de plaats van het misdrijf? Zijn de log returns van beurswaarden op verschillende dagen onafhankelijk? Et cetera. Antwoorden op dergelijke vragen worden gebaseerd op de uitkomsten van experimenten of onderzoeken. In veel gevallen laten de uitkomsten van die experimenten echter geen ondubbelzinnig antwoord toe. Als een nieuwe therapie bij 100 patiënten wordt beproefd, en in 64 gevallen goede resultaten geeft, terwijl dat bij de oude therapie in 50% van de patiënten het geval is, is de nieuwe therapie dan ook werkelijk beter dan de oude of hebben we “geluk” gehad? Als bij 75 van de 100 patiënten verbetering optreedt kan van toeval toch geen sprake meer zijn, . . . of toch wel? Is een waargenomen steekproefcorrelatiecoëfficiënt van 0.17 “significant” verschillend van 0? De toetsingstheorie is erop gericht dit beslissingsproces, waarin gekozen moet worden tussen twee conflicterende hypothesen, te formaliseren.

5.1

Nulhypothese en Alternatieve Hypothese

De beslissing tussen conflicterende hypothesen wordt gebaseerd op een geschikt statistisch model voor de waarneming X. De hypothesen worden gecodeerd in parameterwaarden die de kansverdelingen in het statistische model indiceren. We zullen ons hier beperken tot twee hypothesen. De parameter θ behoort ofwel tot een verzameling Θ0 , corresponderend met de ene hypothese, ofwel tot het complement Θ1 = Θ \ Θ0 , waarbij Θ = Θ0 ∪ Θ1 een disjuncte splitsing van de gehele parameterruimte Θ is. We noemen de hypothese H0 : θ ∈ Θ0 de nulhypothese en de hypothese H1 : θ ∈ Θ1 de alternatieve hypothese. In de standaard toetsingsaanpak (die door de meeste gebruikers van de statistiek wordt gevolgd) worden de nul- en alternatieve hypothese niet symmetrisch behandeld. Het gaat er vooral om te weten te komen of de alternatieve hypothese juist is. Mocht de data hiervoor onvoldoende aanwijzing opleveren, dan besluiten we niet noodzakelijkerwijze dat de alternatieve hypothese onjuist is (en de nulhypothese juist); het is ook mogelijk dat voor geen van beide hypothesen voldoende bewijs is. De statistische analyse kan zo tot twee conclusies leiden: - Verwerp H0 (en accepteer H1 als zijnde correct). - Verwerp H0 niet (maar accepteer H0 niet als correct). De eerst mogelijke conclusie is een sterke conclusie; de tweede is eigenlijk geen conclusie. De tweede conclusie moet worden begrepen als een uitspraak dat meer informatie nodig is om tot een conclusie te komen.

5.2: Toetsingsgrootheid en Kritiek Gebied

25

Doen we op grond van onze waarnemingen uitspraken over de hypothesen, dan kunnen we twee soorten fouten maken, corresponderend met het ten onrechte besluiten tot één van de twee mogelijke conclusies: - Een fout van de eerste soort is H0 verwerpen als H0 correct is. - Een fout van de tweede soort is H0 niet verwerpen als H0 incorrect is. Een fout van de eerste soort correspondeert met het ten onrechte besluiten tot de sterke conclusie. We beschouwen dit als zeer ongewenst. Een fout van de tweede soort correspondeert met het ten onrechte besluiten tot de zwakke conclusie. Dit is ongewenst, maar omdat de zwakke conclusie eigenlijk geen conclusie is, minder erg. Vanwege de asymmetrische behandeling van de hypothesen H0 en H1 bij de keuze van een toets, mag aan het niet-verwerpen van H0 niet veel betekenis worden gehecht. Het is daarom van groot belang de nulhypothese en de alternatieve hypothese geschikt te kiezen. In principe kiezen we de uitspraak die we willen aantonen als alternatieve hypothese. Vervolgens stellen we ons als het ware op het H 0 -standpunt: we verwerpen H0 alleen als er sterke aanwijzingen tegen H0 zijn.

-0.6

-0.4

-0.2

0.0

0.2

0.4

Voorbeeld 5.1 (Twee steekproeven). Figuur 5.1 toont boxplots van de mate van expressie van een gen in twee verschillende typen tumoren. De steekproeven zijn respectievelijk 26 en 15 tumoren groot. De vraag is of het gen in het ene type tumor meer tot expressie is gekomen dan in het andere type tumor. De boxplot geeft niet onmiddellijk uitsluitsel op deze vraag. Weliswaar ligt de box van de tweede steekproef hoger dan die voor de eerste, maar er is een duidelijke overlap en het bereik van de tweede steekproef ligt duidelijk binnen het bereik van de eerste steekproef. Dit laatste zou betekenis kunnen hebben, maar ook een gevolg kunnen zijn van het feit dat de steekproeven ongelijke groottes hebben. Een formele toets kan helpen de vraag te beantwoorden. Een redelijk statistisch model is dat de twee steekproeven X1 , . . ., X26 en Y1 , . . ., Y15 onafhankelijke steekproeven zijn uit de normale verdelingen met, respectievelijk, parameters (µ, σ 2 ) en (ν, τ 2 ). We willen dan de nulhypothese H0 : µ = ν tegen het alternatief H1 : µ 6= ν toetsen. De parameter kunnen we hier gelijk nemen aan θ = (µ, ν, σ 2 , τ 2 ) met parameterruimte Θ = R2 × (0, ∞)2 . De nulhypothese is de deelverzameling Θ0 = {(µ, µ): µ ∈ R} × (0, ∞)2 .

Figuur 5.1. Boxplots van de mate van expressie van een gen gemeten in twee groepen van 26 (links) en 15 (rechts) tumoren.

5.2

Toetsingsgrootheid en Kritiek Gebied

Op basis van de waarneming X moeten we besluiten of er voldoende aanwijzingen tegen de nulhypothese H0 zijn, zodat we H0 willen verwerpen en de bewering onder de alternatieve

26

5: Toetsen

hypothese als correct willen beschouwen. De waarden van X waarvoor de aanwijzingen sterk genoeg zijn vormen het zogenaamde kritieke gebied K. Voor deze waarden van X hebben we voldoende vertrouwen in de alternatieve hypothese om H0 te verwerpen. Definitie 5.2. Bij een gegeven nulhypothese H0 wordt een statistische toets (Engels: test) gegeven door een een verzameling K van mogelijke waarden van de waarneming X, het kritieke gebied. Veronderstel dat x is waargenomen. Als x ∈ K, dan verwerpen we H0 ; als x ∈ / K, dan verwerpen we H0 niet. Met name als X = (X1 , . . ., Xn ) een vector van waarnemingen is, is het vaak lastig om op basis van X in te zien of de uitspraak onder de alternatieve hypothese juist kan zijn. Daarom vatten we de data vaak samen in een toetsingsgrootheid. Een toetsingsgrootheid is een reëelwaardige grootheid T = T (X) die gebaseerd is op de data en informatie geeft over de juistheid van de nul- en de alternatieve hypothese; de toetsingsgrootheid hangt dus niet van de onbekende parameter af. Het kritieke gebied K heeft veelal de vorm {x: T (x) ∈ KT }, of kortweg {T ∈ KT }, voor een toetsingsgrootheid T en een verzameling KT in het bereik van T . In de praktijk wordt de verzameling KT vaak ook wel het kritieke gebied genoemd. Hoe het kritieke gebied K of K T opgesteld kan worden, bespreken we in de volgende paragraaf. Voorbeeld 5.3 (Gauss-toets). Veronderstel dat X1 , . . ., Xn een steekproef vormen uit de normale verdeling met onbekende verwachting µ en bekende variantie σ 2 . We willen graag de nulhypothese H0 : µ ≤ µ0 toetsen tegen het alternatief H1 : µ > µ0 , voor µ0 een vast gekozen getal, bijvoorbeeld voor µ0 = 0. Dit probleem doet zich bijvoorbeeld voor bij de kwaliteitscontrole van producten in een fabriek. Omdat de fabrikant het te duur vindt om alle producten te controleren, wordt een kwaliteitsmaat bij een steekproef van producten gemeten. Uit eerder onderzoek is bekend dat de kwaliteitsmaatstaf normaal verdeeld is. De fabrikant wil bevestigen dat de gemiddelde kwaliteit van de gehele productie groter dan µ0 is. (De aanname van een bekende σ 2 is niet realistisch, maar vergemakkelijkt het voorbeeld. In de praktijk veronderstelt men dat σ 2 onbekend is en gebruikt men bijna altijd de t-toets uit Voorbeeld 5.16.) Het gemiddelde X is de maximum likelihood-schatter voor µ en kan daarom gebruikt worden om een idee te krijgen over de juistheid van de nul- en de alternatieve hypothese. Als het waargenomen gemiddelde x groter is dan µ0 dan wijst dit op het waar zijn van de alternatieve hypothese en hoe groter x is hoe sterker deze aanwijzing is. Het gemiddelde X kunnen we dus gebruiken als toetsingsgrootheid en we verwerpen H0 voor grote waarden van de toetsingsgrootheid. Het kritieke gebied heeft dan de vorm K = {(x1 , . . ., xn ): x ≥ c} voor een zekere waarde c. Maar, hoe groot moet c worden genomen opdat we genoeg vertrouwen hebben in de juistheid van de alternatieve hypothese als x ≥ c en de fout van de eerste soort klein genoeg is? Stel dat voor een statistische toets het kritieke gebied de vorm K = {x: T (x) ∈ K T } heeft waar T een toetsingsgrootheid is en KT een deelverzameling in het bereik van T . De verzameling KT hangt af van de keuze van de toetsingsgrootheid T . Bij een andere toetsingsgrootheid T 0 hoort in het algemeen een andere verzameling KT 0 . Het kritieke gebied K kan echter in beide gevallen hetzelfde zijn; bij twee verschillende toetsingsgrootheden kan hetzelfde kritieke gebied K horen (zie Opgave 5.8). 5.2.1

Onbetrouwbaarheid en Onderscheidend Vermogen

Wanneer bij het toetsen van H0 : θ ∈ Θ0 tegen H1 : θ ∈ Θ1 de ware waarde van θ tot Θ0 behoort, is de nulhypothese waar. Als in dat geval x ∈ K, dan verwerpen we H0 ten onrechte en maken we een fout van de eerste soort. Voor een goede toets moet daarom de kans P θ (X ∈ K) voor θ ∈ Θ0 klein zijn. Aan de andere kant willen we in het geval dat de nulhypothese niet waar is (θ ∈ Θ1 ) juist dat Pθ (X ∈ K) groot is. De kwaliteit van een toets kan daarom worden afgemeten aan de functie θ 7→ Pθ (X ∈ K).


27

Definitie 5.4. Het onderscheidend vermogen (Engels: power function) van een toets K is π(θ; K) = Pθ (X ∈ K).

0.0

0.2

0.4

0.6

0.8

1.0

We zoeken dus een kritiek gebied waarvoor het onderscheidend vermogen “kleine waarden” (dichtbij 0) aanneemt als θ ∈ Θ0 , en “grote waarden” (dichtbij 1) als θ ∈ Θ1 . In Figuur 5.2 zijn van twee toetsen het onderscheidend vermogen weergegeven (als functie van θ op de horizontale as), een “ideale toets” met kans op beide soorten fouten gelijk aan 0 en een reële toets.

Figuur 5.2. Onderscheidend vermogen van een ideale toets (doorgetrokken) en een re¨ ele toets (gestippeld). De parameterruimte onder de nul- en de alternatieve hypothese (Θ0 en Θ1 ) zijn het gedeelte van de horizontale as waar het onderscheidend vermogen van de ideale toets gelijk aan 0, respectievelijk gelijk aan 1 is.

Definitie 5.5. De onbetrouwbaarheid (Engels: size) van een toets K met onderscheidend vermogen π(·; K) is het getal α = sup π(θ; K). θ∈Θ0

Een toets is van niveau α0 (Engels: level α0 ) als α ≤ α0 . De asymmetrie tussen de twee hypothesen wordt nu formeel gemaakt door een afspraak die zeker stelt dat de kans op een fout van de eerste soort hoogstens α0 is. Afspraak 5.6. In iedere praktijksituatie kiezen we eerst een vast getal α 0 , de onbetrouwbaarheidsdrempel. Vervolgens gebruiken we alleen toetsen van niveau α 0 . Met andere woorden, we laten alleen toetsen toe waarvan het onderscheidend vermogen π(·; K) onder de nulhypothese ten hoogste α0 is: sup π(θ; K) ≤ α0 . θ∈Θ0

Het lijkt aantrekkelijk de onbetrouwbaarheidsdrempel α0 extreem klein te kiezen, zodat we vrijwel nooit een fout van de eerste soort maken. Dit kunnen we alleen bereiken door K bijzonder klein te maken. In dat geval zal echter ook het onderscheidend vermogen voor θ ∈ Θ 1 klein worden. De kansen op een fout van de tweede soort Pθ (X ∈ / K) = 1 − π(θ; K),

θ ∈ Θ1 ,

worden nu dus erg groot, hetgeen ook ongewenst is. De eisen om de kansen op fouten van eerste en tweede soort beide klein te maken werken elkaar tegen. We behandelen de twee soorten fouten niet symmetrisch; we streven er bijvoorbeeld niet naar de som van de maximale kansen op fouten van de eerste en tweede soort te minimaliseren. In de praktijk kiest men α0 vaak gelijk aan het magische getal 0.05. Deze keuze leidt ertoe dat, als we vaak toetsen, we ons niet mogen verbazen als we 1 op de 20 keer de nulhypothese ten onrechte verwerpen (en een fout van de eerste soort maken). Eigenlijk zouden we α 0 afhankelijk

28

5: Toetsen

moeten kiezen van de mogelijke consequenties van een fout van de eerste soort. Zijn deze buitengewoon ernstig, dan is α0 = 0.05 wellicht veel te groot. Wat betreft de fouten van de eerste soort beschouwen we Afspraak 5.6 als voldoende garantie dat de kans hierop klein is. Meestal zullen veel toetsen (met evt. verschillende toetsingsgrootheden) aan deze afspraak voldoen. Van deze toetsen geven we de voorkeur aan die toets met de kleinste kansen op een fout van de tweede soort. Hoe klein deze kansen zijn hangt af van de situatie, onder meer van het aantal waarnemingen en de gekozen onbetrouwbaarheidsdrempel α0 . Bij te grote kansen op fouten van de tweede soort is de toets natuurlijk weinig zinvol, omdat we dan bijna altijd H0 niet verwerpen en de tweede, niet-conclusie zullen kiezen. Afspraak 5.7. Gegeven de onbetrouwbaarheidsdrempel α0 , geven we de voorkeur aan een toets van niveau α0 met een zo groot mogelijk onderscheidend vermogen π(θ; K) onder θ ∈ Θ1 . Volgens deze afspraak prefereren we bij een gegeven onbetrouwbaarheidsdrempel α 0 een toets K1 boven een toets K2 , als beide toetsen van niveau α0 zijn en K1 een groter onderscheidend vermogen bezit dan K2 voor alle θ ∈ Θ1 : sup π(θ; Ki ) ≤ α0 , i = 1, 2

θ∈Θ0

en

π(θ; K1 ) ≥ π(θ; K2 ),

∀θ ∈ Θ1 ,

met strikte ongelijkheid voor tenminste één θ ∈ Θ1 . We noemen de toets K1 meer onderscheidend (Engels: more powerful) dan de toets K2 in zekere θ ∈ Θ1 als π(θ; K1 ) > π(θ; K2 ). We noemen K1 uniform meer onderscheidend als de ongelijkheid geldt voor alle θ ∈ Θ1 . In principe zoeken we nu de uniform meest onderscheidende toets van niveau α0 ; dat is een toets waarvan (bij een gegeven onbetrouwbaarheidsdrempel) het onderscheidend vermogen maximaal is voor alle θ ∈ Θ1 . We vergelijken hier twee functies, en het is mogelijk dat de ene toets meer onderscheidend is voor bepaalde θ ∈ Θ1 , en de andere toets juist voor andere θ ∈ Θ1 . Welke toets we dan moeten prefereren is niet onmiddellijk duidelijk. Deze vraag komt in dit boek niet aan de orde. Voorbeeld 5.8 (Binomiale toets). Zij p de kans op succes bij een nieuwe therapie voor een willekeurig gekozen patiënt. Men wil deze therapie vergelijken met een oude therapie, die slechts in de helft van de gevallen succesvol is. Omdat men alleen ge¨ınteresseerd is in de nieuwe therapie wanneer die beter blijkt dan de oude, vergelijken we de onbekende succeskans p van de nieuwe therapie met 0.5; de (bekende) succeskans van de bestaande therapie is. We willen “bewijzen” dat de nieuwe therapie beter is dan de oude. We nemen de uitspraak “p > 0.5” daarom als alternatieve hypothese. De nul- en alternatieve hypothese zijn dan gelijk aan H 0 : p ≤ 0.5 en H1 : p > 0.5. Wanneer we H0 kunnen verwerpen, zullen we overgaan op de nieuwe therapie. De nieuwe therapie wordt bij 100 patiënten toegepast. We noteren het aantal patiënten waarvoor de nieuwe therapie succesvol is als de waarneming X en veronderstellen dat X bin(100, p)verdeeld is. Het ligt voor de hand T (X) = X als toetsingsgrootheid te nemen, en het kritieke gebied van de vorm K = {cα0 , cα0 + 1, . . ., 100}.

Een grote waarde van X geeft immers een aanwijzing dat H0 onjuist is. De waarde cα0 moet zo worden gekozen dat de onbetrouwbaarheid van de toets ten hoogste α 0 is. De onbetrouwbaarheid van de toets wordt gegeven door α = sup Pp (X ≥ cα0 ) = P0.5 (X ≥ cα0 ). p≤0.5

Het supremum wordt aangenomen in p = 0.5, omdat Pp (X ≥ cα0 ) als functie van p monotoon stijgend is. Deze monotonie is met enige moeite analytisch te bewijzen, maar is ook intu¨ıtief duidelijk. De functie p 7→ Pp (X ≥ cα0 ) is voor cα0 = 59 getekend in Figuur 5.3. Veronderstel dat we α0 = 0.05 kiezen. Als we vervolgens c0.05 = 59 nemen, is de onbetrouwbaarheid α = P0.5 (X ≥ 59) = 0.044 kleiner dan α0 = 0.05, terwijl voor c0.05 = 58 de onbetrouwbaarheid P0.5 (X ≥ 58) = 0.067 > 0.05. Voor c0.05 ≤ 58 is de toets dus niet van niveau 0.05 en daarom niet toegestaan bij deze waarde van de onbetrouwbaarheidsdrempel. We

29

0.0

0.2

0.4

0.6

0.8

1.0


0.0

0.2

0.4

0.6

0.8

1.0

p Figuur 5.3. De functie p 7→ Pp (X ≥ 59) voor X bin(100, p) verdeeld.

0.0

0.2

0.4

0.6

0.8

1.0

moeten daarom c0.05 ≥ 59 kiezen. Ter illustratie geeft Figuur 5.4 de functie x 7→ P0.5 (X ≥ x). Volgens Afspraak 5.7 moeten we het kritieke gebied zo kiezen dat het onderscheidend vermogen zo groot mogelijk is. Dit komt er op neer dat we het kritieke gebied zo groot mogelijk moeten kiezen, zodat onder H1 de kans op het (terecht) verwerpen van de nulhypothese, Pp (X ∈ K), zo groot mogelijk is. We kiezen daarom K = {59, 60, . . ., 100}. Onder alle toetsen van de gegeven vorm is dit de toets van niveau 0.05 met het grootste onderscheidend vermogen. De functie p 7→ Pp (X ≥ 59) in Figuur 5.3 is precies het onderscheidend vermogen van deze toets. Vindt men 64 successen bij de nieuwe therapie, dan wordt H0 dus verworpen op niveau 0.05 en luidt de conclusie dat de nieuwe therapie een grotere succeskans heeft dan de oude therapie. Bij 58 successen hadden we deze conclusie niet kunnen trekken: H0 was dan niet verworpen.

0

20

40

60

80

100

x Figuur 5.4. De functie x 7→ P0.5 (X ≥ x) voor X binomiaal verdeeld met parameters 100 (en 0.5). Deze functie is linkscontinu in punten waar x een waarde in N aanneemt. De horizontale lijn is ter hoogte 0.05.

In het geval van een eendimensionale parameter θ spreken we van een eenzijdige hypothese wanneer de nulhypothese van de vorm H0 : θ ≤ θ0 of H0 : θ ≥ θ0 is, waarbij θ0 een vast gegeven getal is. De alternatieve hypothese is dan van de vorm H1 : θ > θ0 respectievelijk H1 : θ < θ0 . De eerste hypothese noemen we rechtseenzijdig en de tweede hypothese linkseenzijdig. Wanneer de nul- en alternatieve hypothese de vorm H0 : θ = θ0 en H1 : θ 6= θ0 hebben, spreken we van een tweezijdige hypothese.

30

5: Toetsen

Voor een toetsingsgrootheid T heeft het kritieke gebied vaak een van de volgende vormen KT = {T ≥ cα0 },

KT = {T ≤ cα0 }, KT = {T ≤ cα0 } ∪ {T ≥ dα0 },

voor getallen cα0 en dα0 met cα0 < dα0 in het laatste kritieke gebied. Welke vorm het kritieke gebied aanneemt is afhankelijk van de opgestelde hypothesen en de keuze voor de toetsingsgrootheid. De eerste twee vormen van KT heten eenzijdig, de laatste tweezijdig. De getallen cα0 en dα0 heten de kritieke waarden. Is de waarde van de toetsingsgrootheid extremer dan de kritieke waarde, dan wordt de nulhypothese verworpen. Merk op dat “extremer” zowel “groter dan” als “kleiner dan” kan betekenen, afhankelijk van de context en de toetsingsgrootheid. De Gauss-toets in Voorbeeld 5.9 is een voorbeeld van een toets waar een eenzijdige nulhypothese leidt tot een eenzijdig kritiek gebied KT en een tweezijdige nulhypothese tot een tweezijdig kritiek gebied KT . Dit is echter niet in het algemeen het geval; de vorm van het kritieke gebied hangt af van de opgestelde hypothesen en de keuze van de toetsingsgrootheid. Voorbeeld 5.9 (Gauss-toets, vervolg). Stel dat X1 , . . ., Xn een steekproef is uit de N (µ, σ 2 )verdeling, waarbij σ 2 een bekende constante is. Beschouw het toetsingsprobleem H0 : µ ≤ µ0 tegen H1 : µ > µ0 , waarbij µ0 een vast gegeven getal is (bijvoorbeeld µ0 = 0). In Voorbeeld 5.3 hebben we gezien dat het gemiddelde X een geschikte toetsingsgrootheid zou kunnen zijn. Het blijkt echter handiger te zijn om deze grootheid te standaardiseren tot T =

√ X − µ0 n , σ

zodat T onder de aanname µ = µ0 een N (0, 1)-verdeling bezit. Zowel µ0 als σ 2 is bekend, waardoor T ook daadwerkelijk een toetsingsgrootheid is. De toetsingsgrootheid X leidt tot hetzelfde kritieke gebied K (zie Opgave 5.8). Grote waarden van X, groter dan µ0 , en dus van T , zijn waarschijnlijker onder H1 dan onder H0 . Immers X is normaal verdeeld met verwachting µ en variantie σ 2 /n, en deze verdeling schuift naar rechts als µ toeneemt. We kiezen daarom een kritiek gebied, gebaseerd op de toetsingsgrootheid T , van de vorm K = (x1 , . . ., xn ): T ≥ cα0 . In de volgende twee alinea’s beargumenteren we dat de juiste keuze voor cα0 wordt gegeven door het (1−α0 )-kwantiel ξ1−α0 van de standaard normale verdeling. (We noteren met ξα het getal zodanig dat Φ(ξα ) = α, waarbij Φ de standaard normale verdelingsfunctie is.) Volgens Afspraak 5.6 zoeken we een toets met een onbetrouwbaarheid die hoogstens α 0 is, dat wil zeggen: sup Pµ ((X1 , . . ., Xn ) ∈ K) = sup Pµ (T ≥ cα0 ) ≤ α0 .

(5.1)

µ≤µ0

µ≤µ0

√ Omdat n(X − µ)/σ, als µ de ware parameterwaarde is, de standaard normale verdeling volgt, is de kans Pµ (T ≥ cα0 ) gelijk aan Pµ

√ X − µ √ X − µ √ µ0 − µ 0 ≥ c α0 = P µ n ≥ c α0 + n n σ σ σ √ µ0 − µ = 1 − Φ cα0 + n . σ

Deze kans is een stijgende functie van µ (hetgeen ook intu¨ıtief duidelijk is uit het feit dat de normale verdeling met verwachting µ naar rechts schuift als µ toeneemt), zodat het supremum supµ≤µ0 Pµ (T ≥ cα0 ) wordt aangenomen voor de grootst mogelijke waarde van µ, µ = µ0 . De eis (5.1) dat de onbetrouwbaarheid hoogstens α0 is, reduceert dus tot Pµ0 (T ≥ cα0 ) ≤ α0 . Aangezien T standaard normaal verdeeld is onder de aanname dat µ = µ0 , volgt hieruit dat cα0 ≥ ξ1−α0 .


31

Onder de toetsen van niveau α0 (van bovenstaande vorm) zoeken we nu de meest onderscheidende toets, volgens Afspraak 5.7. Dit is natuurlijk de toets met het grootste kritieke gebied, dat wil zeggen met een zo klein mogelijke kritieke waarde cα0 . In combinatie met de ongelijkheid in de vorige alinea nemen we cα0 = ξ1−α0 . Merk op dat de onbetrouwbaarheid nu precies gelijk is aan de onbetrouwbaarheidsdrempel α0 . Samengevat, de toets verwerpt de nulhypothese H0 : µ ≤ µ0 voor waarden van X zodanig √ dat T = n(X − µ0 )/σ ≥ ξ1−α0 . Dit is de gebruikelijke toets voor dit probleem, de Gauss-toets (genoemd naar de wiskundige die als een van de eersten de normale verdeling hanteerde). Het bijbehorende kritieke gebied is gelijk aan o n √ x − µ0 ≥ ξ1−α0 . K = {(x1 , . . ., xn ): T ∈ KT } = (x1 , . . ., xn ): n σ

0.0

0.2

0.4

0.6

0.8

1.0

De verzameling KT is dus gelijk aan [ξ1−α0 , ∞). Merk op dat de gevonden kritieke waarde cα0 = ξ1−α0 niet afhangt van de waarden van µ0 en σ 2 . Voor alle waarden van µ0 en σ 2 wordt hetzelfde kritieke gebied KT = [ξ1−α0 , ∞) gevonden. Dit is het voordeel van de gestandaardiseerde toetsingsgrootheid T boven X als toetsingsgrootheid. Daarom is het bij de Gauss-toets gebruikelijk de gestandaardiseerde toetsingsgrootheid te gebruiken. De verzameling KT = [ξ1−α0 , ∞) wordt dan ook vaak het kritieke gebied genoemd van de rechtseenzijdige Gauss-toets. Op analoge wijze kan de nulhypothese H0 : µ ≥ µ0 worden getoetst tegen de alternatieve hypothese H1 : µ < µ0 . Voor deze toets hanteert men dezelfde toetsingsgrootheid T . De nul√ hypothese H0 wordt bij onbetrouwbaarheidsdrempel α0 verworpen als T = n(X − µ0 )/σ ≤ ξα0 = −ξ1−α0 . Het kritieke gebied voor het toetsen van de nulhypothese H0 : µ = µ0 tegen het tweezijdige alternatief H1 : µ 6= µ0 bij onbetrouwbaarheidsdrempel α0 wordt gevonden door samenvoeging van de kritieke gebieden van de beide eenzijdige α 0 /2. Dit √ √ toetsen met elk onbetrouwbaarheid leidt tot verwerping van de nulhypothese als n(X −µ0 )/σ ≤ ξα0 /2 of n(X −µ0 )/σ ≥ ξ1−α0 /2 , √ of, equivalent, als n|X − µ0 |/σ ≥ ξ1−α0 /2 .

−3

−2

−1

0

1

2

3

mu

Figuur 5.5. Onderscheidend vermogens als functie van µ van de twee eenzijdige Gauss-toetsen (gestreept en gestippeld) en de tweezijdige Gauss-toets (doorgetrokken) voor µ0 = 0 bij α0 = 0.05 en n = 5.

Uiteraard bezit de tweezijdige toets een kleiner onderscheidend vermogen dan de linkseenzijdige toets voor waarden µ < µ0 en dan de rechtseenzijdige toets voor waarden µ > µ0 , zie Figuur 5.5. Is men alleen in één van deze typen alternatieven ge¨ınteresseerd, dan verdient een geschikte eenzijdige toets dus de voorkeur boven een tweezijdige toets. Dit kan bijvoorbeeld het geval zijn als men overweegt een nieuwe productiemethode in te voeren of een nieuw apparaat aan te schaffen. Men is dan niet zozeer ge¨ınteresseerd in de vraag of deze innovatie tot achteruitgang kan leiden, maar men wil weten of een verbetering te verwachten is. De keuze tussen eenzijdig en tweezijdig toetsen hangt dus af van de praktische vraagstelling. Als men het idee achter “onbetrouwbaarheid” serieus wil nemen, dan mag men zich bij die keuze niet

32

5: Toetsen

laten leiden door de uitkomsten van de experimenten! In het bijzonder zou het verkeerd zijn bijvoorbeeld voor de rechtseenzijdige toets te kiezen nadat is vastgesteld dat X > µ 0 . In het bovenstaande hebben we de Gauss-toets ge¨ıntroduceerd middels een ad hoc argument. Behalve intu¨ıtief redelijk zijn deze toetsen ook de best mogelijke. Men kan namelijk bewijzen dat de eenzijdige Gauss-toetsen uniform meest onderscheidend zijn; dat wil zeggen dat bij deze toetsen het onderscheidend vermogen in alle mogelijke waarden onder de alternatieve hypothese maximaal zijn. De tweezijdige Gauss-toets is uniform meest onderscheidend onder de zuivere toetsen. Zuivere toetsen zijn toetsen met π(θ0 ) ≤ α0 ≤ π(θ1 ) voor alle θ0 ∈ Θ0 en θ1 ∈ Θ1 en voor α0 de onbetrouwbaarheidsdrempel. Voorbeeld 5.10 (Binomiale toets, vervolg). Voorbeeld 5.8 betreft een speciaal geval van de volgende binomiale toets. Veronderstel dat voor een vast gekozen getal p 0 ∈ (0, 1) we de nulhypothese H0 : p ≤ p0 willen toetsen tegen H1 : p > p0 op grond van een bin(n, p)-verdeelde waarneming X. We kiezen X zelf als toetsingsgrootheid en verwerpen H0 voor grote waarden van X. Het kritieke gebied heeft derhalve de vorm {x ∈ {0, 1, . . ., n}: x ≥ cα0 } = {cα0 , . . ., n}. We kiezen de kritieke waarde cα0 ∈ {0, . . ., n} zodanig dat de onbetrouwbaarheid van de toets kleiner dan of gelijk is aan α0 en, onder deze nevenvoorwaarde, het onderscheidend vermogen maximaal is (vergelijk Voorbeeld 5.8). De onbetrouwbaarheid van deze toets is gelijk aan α = sup Pp (X ≥ cα0 ) = Pp0 (X ≥ cα0 ), p≤p0

aangezien de kans Pp (X ≥ x) stijgend is in p bij vaste x. Om het onderscheidend vermogen zo groot mogelijk te maken onder de alternatieve hypothese, nemen we kriteke gebied zo groot mogelijk, ofwel de kritieke waarde zo klein mogelijk: o n cα0 = min t ∈ {0, . . ., n}: Pp0 (X ≥ t) ≤ α0 . Uiteraard geldt dan dat α ≤ α0 . Vanwege het sprongkarakter van de binomiale verdelingsfunctie zal deze ongelijkheid strikt zijn voor de meeste waarden van α0 . Voor niet te kleine waarden van n kunnen we de kans Pp0 (X ≥ t) normaal benaderen en zijn de sprongen in de verdelingsfunctie van X te verwaarlozen. Voor de onbetrouwbaarheid van de binomiale toets levert dit α0 ≥ Pp0 X ≥ cα0 = Pp0 X ≥ cα0 − 12 X − np c − np − 1 cα − np0 − 21 0 α 0 2 = P p0 p ≈ 1 − Φ p0 ≥ p0 np0 (1 − p0 ) np0 (1 − p0 ) np0 (1 − p0 ) waarbij het ≈-teken volgt uit de benadering van de binomiale verdelingsfunctie door de normale verdelingsfunctie en de term 1/2 in de teller de continu¨ıteitscorrectie is (zie Appendix 8). Bij gegeven α0 is de waarde van cα0 het kleinste gehele getal waarvoor geldt dat (5.2)

cα − np0 − 21 ξ1−α0 ≤ p0 . np0 (1 − p0 )

Aanpassingen van deze eenzijdige toets voor het geval van het andere eenzijdige probleem, H1 : p < p0 , of het tweezijdige probleem, H1 : p 6= p0 , liggen voor de hand.

5.3

Statistische Significantie

De algemene opzet van de toetsingstheorie zoals hiervoor is beschreven is zowel tamelijk ingewikkeld, als verbluffend eenvoudig, omdat er slechts twee beslissingen mogelijk zijn. In veel praktijksituaties is de eenvoud misleidend. Een effect wordt statistisch significant genoemd als de relevante nulhypothese wordt verworpen, bij de gegeven onbetrouwbaarheidsdrempel. Dit

5.4: Overschrijdingskansen

33

moet als volgt worden ge¨ınterpreteerd: het effect dat we in de data hebben waargenomen is waarschijnlijk niet aan toevalsvariatie te wijten; zouden we het hele experiment herhalen, dan zouden we waarschijnlijk hetzelfde effect opnieuw vinden. Dit hoeft geenszins te betekenen dat het “effect” praktisch significant is. Het is heel denkbaar dat de toetsingsprocedure terecht heeft aangetoond dat de nieuwe therapie beter is, maar dat de verbetering verwaarloosbaar klein is. Als de oude therapie kans p = 0.5 op succes heeft, en de nieuwe kans p = 0.500001, dan zullen we dit effect vinden en H0 verwerpen mits we voldoende waarnemingen doen, maar praktisch gesproken zal het ons waarschijnlijk weinig uitmaken welke therapie we volgen. Om deze reden is het wenselijk een toetsingsprocedure die leidt tot verwerping van H 0 altijd aan te vullen met een schattingsprocedure die een indicatie geeft van de grootte van een mogelijk effect. De context bepaalt dan of dit effect van praktisch belang is. Een andere mogelijkheid om de discrepantie tussen statistische en praktische significantie te overbruggen zou zijn om de nulhypothese anders te formuleren. We zouden bijvoorbeeld de nulhypothese kunnen toetsen dat het verschil p2 −p1 in kans op succes bij de nieuwe therapie en de oude therapie minstens 0.2 is, in plaats van de hypothese dat p2 − p1 > 0. De waarde 0.2 zou dan de praktische significantie kunnen uitdrukken. In de praktijk is men echter meestal tevreden met het vaststellen van een kwalitatief verschil en toetst men de hypothese H1 : p2 − p1 > 0.

5.4

Overschrijdingskansen

In het voorgaande hebben we toetsen middels een toetsingsgrootheid T en een kritiek gebied K beschreven. Veronderstel dat het kritieke gebied de vorm K = {x: T (x) ≥ dα0 } bezit waarbij de constante dα0 het kleinste getal is zodanig dat een toets van deze vorm niveau α0 heeft. Dat wil zeggen o n (5.3) dα0 = min t: sup Pθ (T ≥ t) ≤ α0 . θ∈Θ0

Veelal correspondeert het minimaal nemen van dα0 met het maximaliseren van het onderscheidend vermogen in Θ1 . De formule is daarom een gevolg van Afspraak 5.7. De gelijkheid (5.3) impliceert dat, voor iedere t ∈ R, sup Pθ (T ≥ t) ≤ α0

θ∈Θ0

⇐⇒

t ≥ d α0 .

We kunnen de toets daarom op de volgende wijze uitvoeren: gegeven dat de waarde t is waargenomen voor de toetsingsgrootheid T , bereken de overschrijdingskans of p-waarde (Engels: observed significance level, of, p-value) sup Pθ (T ≥ t).

θ∈Θ0

Is de overschrijdingskans kleiner dan of gelijk aan α0 , dan verwerpen we H0 ; anders verwerpen we H0 niet. Dit voorschrift geeft precies de toets met kritiek gebied K = {x: T (x) ≥ d α0 }, want de overschrijdingskans is kleiner dan of gelijk aan α0 dan en slechts dan als t ≥ dα0 . Bovenstaande wordt met behulp van de Gauss-toets in Figuur 5.6 ge¨ıllustreerd. In de figuur is duidelijk te zien dat voor waarden t in het kritieke gebied geldt dat supµ≤µ0 Pµ0 (T ≥ t) ≤ α0 , en anders om. In woorden is de overschrijdingskans het maximum over alle mogelijkheden onder de nulhypothese van de kans dat bij een identiek experiment een extremere waarde van de toetsingsgrootheid wordt gevonden dan de waarde t van het uitgevoerde experiment. Het supremum over alle mogelijkheden onder de nulhypothese maakt de overschrijdingskans enigszins gecompliceerd. In veel gevallen is het supremum overbodig omdat één van de parameters θ 0 ∈ Θ (vaak een randpunt van Θ0 ) altijd de maximumkans geeft. In dat geval is de overschrijdingskans gelijk aan Pθ0 (T ≥ t).

5: Toetsen

0.0

0.2

0.4

0.6

0.8

1.0

34

−2

0

2

4

Figuur 5.6. Rechter overschrijdingskans t 7→ supµ≤µ0 Pµ (T ≥ t) = Pµ0 (T ≥ t) (doorgetrokken curve) voor de Gauss-toets met µ0 = 0. Op de hoogte van α0 = 0.05 is een stippellijn getekend. De dik gedrukte lijn is het bijbehorende kritieke gebied.

De overschrijdingskans zoals we hem zojuist hebben gedefinieerd is specifiek voor kritieke gebieden van de vorm {x: T (x) ≥ dα0 }. Een uitbreiding naar kritieke gebieden van de vorm {x: T (x) ≤ cα0 } ligt voor de hand, waarbij nu de aanname is dat o n cα0 = max t: sup Pθ (T ≤ t) ≤ α0 . θ∈Θ0

Gegeven de waargenomen waarde t berekenen we de overschrijdingskans sup θ∈Θ0 Pθ (T ≤ t). Is dit getal kleiner dan of gelijk aan α0 , dan verwerpen we H0 . Tweezijdige kritieke gebieden van de vorm {x: T (x) ≤ c} ∪ {x: T (x) ≥ d} bestaan vaak uit een combinatie van twee eenzijdige gebieden in de zin dat c = cα0 /2 en d = dα0 /2 voor cα0 en dα0 als eerder gedefinieerd. De onbetrouwbaarheidsdrempel α0 wordt dus gesplitst in twee gelijke delen van α0 /2 in de linker- en rechterstaart. In dit geval wordt de overschrijdingskans bij waargenomen waarde t gedefinieerd als 2 min sup Pθ (T ≤ t), sup Pθ (T ≥ t) . θ∈Θ0

θ∈Θ0

Is dat getal kleiner dan of gelijk aan α0 , dan verwerpen we H0 ; anders verwerpen we H0 niet. Dit komt neer op het kijken of één van de twee “eenzijdige overschrijdingskansen” kleiner is dan of gelijk is aan α0 /2: 2 min(a, b) ≤ α0 dan en slechts dan als a ≤ α0 /2 of b ≤ α0 /2. Toetsen middels overschrijdingskansen verdient in de meeste gevallen de voorkeur boven toetsen middels een kritiek gebied, omdat de resulterende uitspraak informatiever is. Bij rapportering van de overschrijdingskans is het immers mogelijk alsnog (en op heel eenvoudige wijze) de hypothese bij ieder gewenste onbetrouwbaarheidsdrempel α0 te toetsen, terwijl bij rapportering van het kritieke gebied en de waarde van de toetsingsgrootheid bij een vaste α 0 dit niet mogelijk is. Bovendien geeft, bijvoorbeeld, een heel kleine overschrijdingskans onmiddellijk aan dat H0 overduidelijk wordt verworpen. Voorbeeld 5.11 (Binomiale toets, vervolg). In Voorbeeld 5.8 werd geconcludeerd dat bij 64 successen de nulhypothese wordt verworpen bij α0 = 0.05, terwijl bij 58 successen de nulhypothese niet wordt verworpen. De overschrijdingskansen bij 64 en 58 successen zijn respectievelijk sup Pp (X ≥ 64) = P0.5 (X ≥ 64) = 0.0033 p≤0.5

sup Pp (X ≥ 58) = P0.5 (X ≥ 58) = 0.0666.

p≤0.5

De eerste kans is heel klein en inderdaad kleiner dan 0.05 en de tweede is groter dan 0.05. We zien bovendien dat de nulhypothese bij 64 successen verworpen wordt voor alle onbetrouwbaarheidsdrempels α0 ≥ 0.0033. De overschrijdingskans geeft dus meer informatie dan alleen

5.5: Enkele Standaard Toetsen

35

de vaststelling dat de nulhypothese wordt verworpen bij α0 = 0.05, hetgeen de conclusie was in Voorbeeld 5.8. Voorbeeld 5.12 (Gauss-toets,√vervolg). De Gauss-toets verwerpt de nulhypothese H 0 : µ ≤ µ0 voor grote waarden van T = n(X − µ0 )/σ. De kritieke waarde ξ1−α0 van de toets voldoet aan (5.3). De overschrijdingskans van de toets is daarom gelijk aan, bij waargenomen waarde x, √ x − µ √ x − µ0 √ x − µ0 0 sup Pµ T ≥ n = P µ0 T ≥ n =1−Φ n . σ σ σ µ≤µ0 Is deze kans kleiner dan of gelijk aan α0 , dan wordt H0 verworpen op niveau α0 . De overschrijdingskans voor het toetsen van de andere eenzijdige nulhypothese H 0√ : µ ≥ µ0 tegen de alternatieve hypothese H1 : µ < µ0 wordt gegeven door de kans Pµ0 (T ≤ n(x − µ0 )/σ). We verwerpen de nulhypothese als deze kans kleiner dan of gelijk is aan α0 . De tweezijdige Gauss-toets is niets anders dan de combinatie van de twee eenzijdige toetsen, ieder met onbetrouwbaarheidsdrempel α0 /2. We kunnen deze toets daarom uitvoeren door het berekenen van zowel de linker- als de rechteroverschrijdingskans. De overschrijdingskans van de tweezijdige toets is dan gelijk aan twee maal het minimum van de linker- en rechteroverschrijdingskans. Is één van de twee overschrijdingskansen kleiner dan of gelijk aan α 0 /2, dan is de overschrijdingskans kleiner dan of gelijk aan α0 en verwerpen we de nulhypothese H0 : µ 6= µ0 .

5.5

Enkele Standaard Toetsen

In deze paragraaf bespreken we enkele toetsen die, naast de Gauss-toets en de binomiale toets, veel toegepast worden. De meeste van deze toetsen kunnen op intu¨ıtieve gronden worden begrepen. Het algemene idee is om voor een gegeven toetsingsprobleem een toetsingsgrootheid te vinden die “redelijk” is (vaak gebaseerd op een goede schatter van de parameter) en waarvoor gemakkelijk een kritieke waarde of overschrijdingskans kan worden berekend. Voor dat laatste is het nodig dat de kansverdeling van de toetsingsgrootheid onder de (“rand” van de) nulhypothese getabelleerd of berekenbaar is. Vaak behoort de kansverdeling onder de nulhypothese echter niet tot het gebruikelijke rijtje bekende verdelingen uit de kansrekening. We kunnen dan een nieuwe standaard kansverdeling introduceren en tabelleren. Een alternatief is om de kansverdeling “on-the-fly” te benaderen door stochastische simulatie. We bespreken voorbeelden van beide manieren van aanpak. We beginnen deze paragraaf met een bespreking van de twee belangrijkste statistische kansverdelingen, de chikwadraat- en t-verdelingen. Deze families van kansverdelingen zijn allebei gerelateerd aan de normale verdeling, en treden zowel op bij het toetsen van de parameters van de normale verdeling als bij benaderingen bij grote steekproeven. 5.5.1

Chikwadraat- en t-Verdeling

De chikwadraat- en t-verdelingen zijn continue kansverdelingen, waarvan de dichtheden worden gegeven door relatief eenvoudige uitdrukkingen. Voor ons doel zijn de volgende structurele definities van deze kansverdelingen echter aantrekkelijker. Definitie 5.13. Een stochastische grootheid W bezit dePchikwadraat-verdeling met n vrijheidsgraden, notatie χ2n , als W dezelfde verdeling bezit als ni=1 Zi2 voor Z1 , . . ., Zn een steekproef uit de N (0, 1)-verdeling.

5: Toetsen

0.0

0.05

0.10

0.15

0.20

0.25

36

0

5

10

15

20

2

Figuur 5.7. Dichtheden van de χ -verdelingen met 4 (doorgetrokken) en 10 (gestippeld) vrijheidsgraden.

Definitie 5.14. Een stochastische grootheid T bezit de t-verdeling of Student-verdeling met n vrijheidsgraden, notatie tn , als T dezelfde verdeling bezit als p

Z W/n

,

0.0

0.1

0.2

0.3

0.4

voor Z en W twee onafhankelijke stochastische grootheden uit respectievelijk de N (0, 1)verdeling en de χ2n -verdeling.

−4

−2

0

2

4

Figuur 5.8. Dichtheden van de t-verdelingen met 1 (gestreept) en 5 (gestippeld) en ∞ (doorgetrokken) vrijheidsgraden.

Met standaard technieken uit de kansrekening is het mogelijk om formules voor de dichtheden van chikwadraat- en t-verdeling af te leiden. Deze uitdrukkingen zijn in “klassieke” tijden gebruikt om tabellen van de verdelingsfuncties te maken. Meer recentelijk zijn zij de basis voor standaard algoritmes in statistische software. We zullen de tabellen en software als gegeven beschouwen, en bespreken de precieze vorm van de dichtheden niet. Figuren 5.7 en 5.8 geven een kwalitatief idee van de dichtheden. De volgende stelling laat zien waarom de chikwadraat- en t-verdeling belangrijk zijn. Stelling 5.15. Als X1 , . . ., Xn een steekproef uit de N (µ, σ 2 )-verdeling is, dan geldt (i) X is N (µ, σ 2 /n)-verdeeld. 2 (ii) (n − 1)SX /σ 2 is χ2n−1 -verdeeld. 2 (iii) X en SX zijn onderling onafhankelijk.


(iv)

√

37

n(X − µ)/SX bezit de tn−1 -verdeling.

Bewijs. Bewering (i) is bekend uit de kansrekening: de som van onafhankelijke normaal verdeelde stochastische grootheden is weer normaal verdeeld. Voor het bewijs van beweringen (ii) en (iii) kunnen we zonder verlies van de algemeenheid aannemen dat µ = 0 en σ 2 = 1. De simultane dichtheid van de stochastische vector X = (X1 , . . ., Xn )T is dan gelijk aan (x1 , . . ., xn ) 7→

n Y 1 1 2 2 1 1 √ e− 2 xi = e− 2 kxk , n/2 (2π) 2π i=1

Pn 2 2 waarin kxk √ xi het kwadraat van de Euclidische lengte van x is. Definieer de vector √ = i=1 f1 = (1/ n, . . ., 1/ n) ∈ Rn met kf1 k2 = 1 en vul f1 op een willekeurige wijze aan tot een orthonormale basis {f1 , . . ., fn } van Rn . Zij O de (n × n)-matrix met rijen f1 , . . ., fn . Uit de definitie volgt onmiddellijk dat OO T = I (met I de eenheidsmatrix), zodat O T = O−1 en O een orthogonale matrix is: OO T = OT O = I en dus kOxk2 = xT OT Ox = kxk2 voor alle x. Definieer de stochastische vector Y = OX. Dan is √ Y1 = f1 X = n X, n n X X 2 Yi2 = kY k2 − Y12 = kXk2 − nX = (Xi − X)2 . i=2

i=1

Beweringen (ii) en (iii) volgen daarom als we kunnen bewijzen dat Y1 , . . ., Yn onderling onafhankelijk en N (0, 1)-verdeeld zijn. De verdelingsfunctie van Y wordt gegeven door Z Z 1 1 − 2 kxk2 P(Y ≤ y) = · · · dx1 · · · dxn e n/2 x:Ox≤y (2π) Z Z 1 1 − 2 kuk2 = ··· e du1 · · · dun , n/2 u:u≤y (2π) waar we gebruik maken van de substitutie Ox = u. Dan is kxk = kOxk = kuk en de Jacobiaan van de transformatie Ox = u is gelijk aan det O = 1. Uit de laatste uitdrukking volgt onmiddellijk dat Y dezelfde simultane dichtheid bezit als X. Daaruit volgt dat Y een multivariaat-normale verdeling heeft met verwachtingsvector (0, . . ., 0)T en variantiematrix gelijk aan de eenheidsmatrix I. Het kan bewezen worden dat Y1 , . . ., Yn onderling onafhankelijk en normaal verdeeld met verwachting 0 en variantie 1. Voor het bewijs van (iv) schrijven we √ √ X −µ n(X − µ)/σ = q . n 2 /σ 2 (n−1)SX SX (n−1)

Volgens (iii) zijn de teller en de noemer onafhankelijk, en volgens (i) en (ii) bezitten zij, respectievelijk, een N (0, 1)-verdeling en de wortel van een χ2n−1 -verdeling gedeeld door n − 1. Volgens Definitie 5.14 bezit het quotiënt dan de tn−1 -verdeling. Ieder statistisch pakket bevat functies om de verdelingsfunctie en kwantielen van de Student- en chikwadraat-verdeling numeriek te berekenen. 5.5.2

Eensteekproeftoetsen

Gegeven een steekproef X1 , . . ., Xn wil men vaak toetsen of de locatie van de verdeling van de steekproef zich links of rechts van een bepaalde waarde bevindt. Hierbij wordt “locatie” bijvoorbeeld gepreciseerd als “verwachting”, of als “mediaan”. Nemen we ook aan dat de steekproef uit een normale verdeling afkomstig is, dan is de zeer bekende t-toets de correcte toets voor het probleem wanneer de variantie onbekend is. Als de variantie bekend zou zijn, zouden we gebruikmaken van de Gauss-toets (zie Voorbeeld 5.9).

5: Toetsen

0

0

1

2

2

4

3

4

6

5

38

-1.0

-0.5

0.0

0.5

1.0

0.5

1.0

1.5

2.0

2.5

3.0

Figuur 5.9. Scatterplot van het steekproefgemiddelde (x-as) tegen de steekproefvariantie (y-as) voor 1000 steekproeven ter grootte 5 uit de standaard normale verdeling (links) en 1000 steekproeven ter grootte 10 uit de exponenti¨ ele verdeling (rechts). Links zijn de twee co¨ ordinaten onafhankelijk, rechts bestaat een positief verband.

Voorbeeld 5.16 (t-Toets). Veronderstel dat X1 , . . ., Xn een steekproef is uit de N (µ, σ 2 )verdeling met µ en σ 2 onbekend. Beschouw het toetsingsprobleem H0 : µ ≤ µ0 tegen H1 : µ > µ0 , waarbij µ0 een vast gegeven getal is (bijvoorbeeld µ0 = 0). Formeel gesproken wordt de 2 parameter in dit geval gegeven door het paar θ = (µ, σ ) en is de nulhypothese gelijk aan Θ0 = (µ, σ 2 ): µ ≤ µ0 , σ 2 > 0 . Aangezien de toetsingsgrootheid en het kritieke gebied K van de Gauss-toets uit Voorbeeld 5.9 afhankelijk is van σ en deze parameter nu onbekend is, is die toets hier niet bruikbaar. Een logische uitbreiding van de Gauss-toets is om σ in de definitie van de toetsingsgrootheid te vervangen door een schatter. We gebruiken hiervoor de steekproef standaarddeviatie S X . Dit geeft de toetsingsgrootheid √ X − µ0 T = n . SX We verwerpen de nulhypothese voor grote waarden van deze grootheid. Aangezien de substitutie van SX voor σ ook de verdeling van deze grootheid verandert, is deze niet meer normaal verdeeld als µ = µ0 . Het is daarom niet onmiddellijk duidelijk welke kritieke waarde we dienen te nemen. In de volgende alinea beargumenteren we, met behulp van Stelling 5.15, dat dit het (1 − α0 )-kwantiel van de t-verdeling met n − 1 vrijheidsgraden moet zijn, welke we noteren met tn−1,1−α0 . Volgens Afspraak 5.6 dient de kritieke waarde cα0 , voor het verkrijgen van een toets van niveau α0 , te voldoen aan sup µ≤µ0 ,σ 2 >0

Pµ,σ2

√ X − µ 0 ≥ c α0 ≤ α 0 , n SX

voor α0 de onbetrouwbaarheidsdrempel van de toets. Merk op dat het supremum zowel over µ ≤ µ0 als over alle mogelijke waarden van σ 2 moet worden berekend; over de gehele parameterruimte onder de nulhypothese. Het supremum over µ (voor iedere σ) wordt echter aangenomen in het randpunt µ = µ0 , zoals intu¨ıtief wel duidelijk is (maar niet triviaal om te bewijzen), zodat de ongelijkheid reduceert tot √ X − µ 0 n ≥ c α0 ≤ α 0 . SX σ 2 >0 √ Nu hangt volgens Stelling 5.15(iv) de verdeling van n(X − µ0 )/SX onder (µ0 , σ 2 ) helemaal niet af van (µ0 , σ 2 ) en is gelijk aan de tn−1 -verdeling. Uit bovenstaande ongelijkheid volgt nu dat cα0 ≥ tn−1,1−α0 . Om een zo groot mogelijk onderscheidend vermogen te krijgen, in overeenstemming met Afspraak 5.7, kiezen we het kritieke gebied zo groot mogelijk en nemen we cα0 = tn−1,1−α0 . De onbetrouwbaarheid α van de toets is dan precies gelijk aan de onbetrouwbaarheidsdrempel: α = α0 . sup Pµ0 ,σ2


39

√De resulterende toets, die de t-toets of Student-toets wordt genoemd, luidt: “Verwerp H 0 als n(X − µ0 )/SX ≥ tn−1,1−α0 ”. De corresponderende overschrijdingskans, bij waargenomen waarden x en sx , is gelijk aan √ x − µ0 √ x − µ0 = P Tn−1 ≥ n , Pµ0 ,σ2 T ≥ n sx sx

waarin Tn−1 een stochastische grootheid met de tn−1 -verdeling voorstelt. Aanpassingen van de t-toets voor de toetsingsproblemen H0 : µ ≥ µ0 en H0 : µ = µ0 gaan analoog aan de aanpassingen van de Gauss-toets. Het is belangrijk hierbij te weten dat de t-verdeling symmetrisch is rond 0, net als de normale verdeling, zodat tn,α = −tn,1−α . Voor kleine waarden van n (zoals n ≤ 10) verschilt de tn -verdeling aanzienlijk van de normale verdeling. Het gebruik van normale kwantielen in plaats van tn−1 -kwantielen (dat wil zeggen de Gauss-toets met σ gelijk genomen aan SX ) leidt dan tot een toets met een onbetrouwbaarheid die veel groter is dan de bedoelde onbetrouwbaarheid. Dit is in strijd met Afspraak 5.6. Voor toenemende waarden van n geldt dat de tn -verdeling steeds meer op de standaard normale verdeling gaat lijken, met convergentie naar de normale verdeling voor n → ∞. Voor n ≥ 20 is de gelijkenis al zo goed dat de Student- en Gauss-toetsen in dat geval praktisch identieke resultaten geven. We hebben de t-toets ge¨ıntroduceerd met ad hoc argumenten. Men kan echter laten zien dat de toets uniform meest onderscheidend is binnen de klasse van alle zuivere toetsen. De t-toets is de correcte toets voor het toetsen van locatie in het geval de waarnemingen X1 , . . ., Xn een steekproef uit de normale verdeling vormen. Is aan de laatste aanname niet voldaan, dan is het wellicht mogelijk en wenselijk de waarnemingen te transformeren (bijvoorbeeld met behulp van de logaritmische functie) tot waarnemingen waarvoor de normaliteitsaanname wel redelijk is. Een alternatief is het gebruik van een toets die de normaliteitsaanname niet vereist. Hiervan bestaan veel voorbeelden, waarvan we er slechts één bespreken. Voorbeeld 5.17 (Tekentoets). De tekentoets is toepasbaar onder minimale aannames en is daarom ook geschikt als de verdeling waar de waarnemingen uit afkomstig zijn niet de normale verdeling is. Het is een toets voor de mediaan en niet de verwachtingswaarde, zoals bij de Gaussen de t-toets. Veronderstel dat we willen toetsen of de mediaan µ van de verdeling waar de onafhankelijke waarnemingen X1 , . . ., Xn uit afkomstig zijn, groter is dan een gegeven waarde µ0 : H0 : µ ≤ µ0 tegen H1 : µ > µ0 . De toetsingsgrootheid wordt gegeven door T = #(1 ≤ i ≤ n: Xi > µ0 ); er wordt geteld hoeveel waarnemingen groter zijn dan µ0 , ofwel hoeveel verschillen Xi − µ0 positief zijn. We krijgen nu in feite de binomiale toets toegepast op de tekens (postief of negatief) van de verschillen X1 − µ0 , . . ., Xn − µ0 . De toetsingsgrootheid is binomiaal verdeeld met parameters n en pµ = Pµ (Xi > µ0 ). Als de mediaan van de verdeling van de waarnemingen gelijk is aan µ0 , dan zijn de parameters n en 12 . Bezit de verdeling van de waarnemingen echter een mediaan µ ≤ µ0 , dan is de kans pµ kleiner dan of gelijk aan 1/2. De nulhypothese H0 : µ ≤ µ0 kan daarom worden getoetst door de equivalente nulhypothese H0 : pµ ≤ 1/2 te toetsen op grond van T . We verwerpen voor grote waarden van T , waarbij de kritieke waarde wordt bepaald als in Voorbeeld 5.10. Voorbeeld 5.18 (Toetsen voor σ 2 ). Veronderstel dat X1 , . . ., Xn een steekproef is uit de N (µ, σ 2 )-verdeling met µ en σ 2 onbekende parameters. Beschouw het toetsingsprobleem H0 : σ 2 ≤ σ02 tegen H1 : σ 2 > σ02 , waarbij σ02 een vast gegeven getal is. Formeel gesproken wordt 2 de parameter gegeven door het paar θ = (µ, de nulhypothese σ ) en is de parameterruimte onder 2 2 2 gelijk aan Θ0 = (µ, σ ): µ ∈ R, σ ≤ σ0 . Een redelijke schatter voor σ 2 is de steekproefvari2 2 antie SX . Grote waarden van SX geven een indicatie dat de alternatieve hypothese juist zou 2 kunnen zijn. We verwerpen daarom de nulhypothese voor grote waarden van S X . 2 2 2 De kansverdeling van (n − 1)SX /σ onder (µ, σ ) hangt niet van de parameter (µ, σ 2 ) af, en is precies de chikwadraat-verdeling met n − 1 vrijheidsgraden (zie Stelling 5.15). Noteren we het α-kwantiel van de chi-kwadraat verdeling met n − 1 vrijheidsgraden als χ 2n−1,α , dan is 2 de voor de hand liggende toets: “Verwerp H0 als (n − 1)SX /σ02 ≥ χ2n−1,1−α0 ” (met α0 de onbetrouwbaarheidsdrempel van de toets). Men kan laten zien dat deze toets onbetrouwbaarheid α0 bezit, op dezelfde wijze als in de voorgaande voorbeelden.

40

5: Toetsen

De toetsen voor het andere eenzijdige toetsingsprobleem H0 : σ 2 ≥ σ02 en het tweezijdige toetsingsprobleem H0 : σ 2 = σ02 liggen voor de hand. Merk echter op dat de chikwadraatverdeling niet symmetrisch is en alle kansmassa op (0, ∞) legt. Er is daarom geen directe relatie tussen de kwantielen χ2n−1,α en χ2n−1,1−α , en we kunnen niet de absolute waarde van de toetsingsgrootheid gebruiken om de tweezijdige toets te beschrijven. De tweezijdige toets luidt: 2 2 “Verwerp H0 als (n − 1)SX /σ02 ≤ χ2n−1,α0 /2 of als (n − 1)SX /σ02 ≥ χ2n−1,1−α0 /2 .” We kunnen deze toetsen ook met overschrijdingskansen uitvoeren. 5.5.3

Tweesteekproeventoetsen

Bij het tweesteekproevenprobleem beschikken we over twee steekproeven X 1 , . . ., Xm en Y1 , . . ., Yn uit mogelijk verschillende kansverdelingen, en zijn we ge¨ınteresseerd in het vergelijken van deze kansverdelingen, bijvoorbeeld wat betreft hun locatie. Afhankelijk van de aannames bestaan verschillende typen tweesteekproeventoetsen. We kunnen een belangrijk onderscheid maken tussen toetsen voor gepaarde en ongepaarde waarnemingen. In het eerste geval ontstaan de twee steekproeven uit een steekproef (X1 , Y1 ), . . ., (Xn , Yn ) van paren waarnemingen, waarbij de X- en de Y -variabelen binnen ieder paar gerelateerd kunnen zijn, maar de paren onderling onafhankelijk worden geacht. De Xmeting geeft bijvoorbeeld de toestand van een patiënt voor de behandeling weer en de Y -meting de toestand na de behandeling, bij een onderzoek naar de effectiviteit van een behandeling. Omdat Xi en Yi metingen aan dezelfde patiënt zijn, ligt het voor de hand dat ze stochastisch afhankelijk zijn. Een lage waarde bij de eerste meting is immers een indicatie dat de patiënt in slechte gezondheid verkeert, hetgeen het waarschijnlijk maakt dat de tweede meting eveneens laag zal zijn (ten opzichte van de rest van de populatie, alhoewel misschien wel hoger dan de eerste meting als de behandeling succes heeft). Bij herhaalde metingen aan eenzelfde object of persoon (Engels: repeated measures; longitudinal data) is afhankelijkheid van de metingen niet te vermijden. In andere toepassingen worden de X- en Y -component zelfs bewust stochastisch afhankelijk gemaakt door de opzet van het experiment. Een groep proefpersonen wordt bijvoorbeeld van tevoren gepaard volgens achtergrondvariabelen als leeftijd, sekse, voorgaande behandeling, of ziektegeschiedenis, zodat de twee personen in ieder paar wat deze variabelen betreft vergelijkbaar zijn. Vervolgens krijgt in elk paar één (willekeurig gekozen) persoon het medicijn en de ander een placebo. Een verschil in toestand na deze behandeling geeft een indicatie voor de werkzaamheid van het medicijn. Het doel van het paren van de proefpersonen in deze opzet is om het effect van de behandeling duidelijker naar voren te brengen. Een waargenomen verschil binnen een paar kan immers niet worden verklaard door fluctuaties in de achtergrondvariabelen, maar moet te wijten zijn aan de behandeling (of een nog onbekende achtergrondsvariabele). Paart men de waarnemingen niet, dan kan de extra toevalsfluctuatie, die veroorzaakt is door de achtergrondvariabelen, het behandelingseffect maskeren. Voorbeeld 5.19 (t-Toets voor gepaarde waarnemingen). Het ligt voor de hand om een toets voor het vergelijken van de locaties van twee gepaarde steekproeven (X 1 , Y1 ), . . ., (Xn , Yn ) te baseren op de verschillen Zi = Xi −Yi . De t-toets voor paren is dan de gewone t-toets toegepast op de verschillen Z1 , . . ., Zn . Voor de toepassing van de t-toets veronderstellen we dat de verschillen Z1 , . . ., Zn onafhankelijk en N (∆, σ 2 )-verdeeld zijn, waarbij de parameter ∆ gelijk is aan het verschil EXi − EYi van de verwachtingen. Stel dat we de nulhypothese H0 : ∆ = 0 dat de behandeling geen effect heeft, willen toetsen, dan wel één van de hypothesen H0 : ∆ ≥ 0 of H0 : ∆ ≤ 0. Het onderscheidend vermogen van de t-toets is sterk afhankelijk van de variantie σ 2 . Is de variantie groot, dan is een verschil in verwachtingswaarde moeilijk te detecteren, en is het onderscheidend vermogen van de t-toets gering. Een kleine variantie is gunstig en zorgt voor een groot onderscheidend vermogen. Deze vaststelling maakt duidelijk dat het verstandig kan zijn de steekproeven in het tweesteekproevenprobleem met opzet afhankelijk te maken. Volgens de rekenregel voor varianties geldt immers dat var Zi = var Yi + var Xi − 2 cov(Xi , Yi ), hetgeen kleiner is dan var Yi + var Xi , indien Xi en Yi positief gecorreleerd zijn. Een intu¨ıtieve verklaring is dat het nemen van verschillen toevalsvariatie elimineert die in zowel de X- als


41

de Y -component aanwezig is en waarin we niet ge¨ınteresseerd zijn. Na eliminatie van deze variatie is het gemakkelijker een eventueel verschil te ontdekken dat veroorzaakt wordt door de behandeling. Correcte toepassing van de t-toets vereist wel dat de verschillen Z1 , . . ., Zn kunnen worden opgevat als een steekproef uit een normale verdeling. Voorbeeld 5.20 (Twee steekproeven t-toets). Veronderstel dat de waarnemingen X 1 , . . ., Xm en Y1 , . . ., Yn twee onafhankelijke steekproeven zijn uit, respectievelijk, de N (µ, σ 2 ) en de N (ν, σ 2 )-verdelingen. We willen eenzijdig toetsen of µ−ν > 0: H0 : µ−ν ≤ 0 tegen H1 : µ−ν > 0. Een voor de hand liggende schatter voor µ − ν is het verschil X − Y van de gemiddelden van de steekproeven. Grote waarden van dit verschil zijn een aanwijzing dat H1 juist is. De verdeling van X − Y is normaal met verwachting µ − ν en variantie var(X − Y ) = var X + var Y =

1 σ2 σ2 1 , + = σ2 + m n m n

door onafhankelijkheid van de twee steekproeven. Omdat deze verdeling afhangt van de onbekende parameter σ 2 , kiezen we als toetsingsgrootheid niet X − Y , maar de grootheid T =

X −Y q

SX,Y waarin 2 SX,Y =

1 m

+

1 n

,

m n X X 1 (Xi − X)2 + (Yj − Y )2 m + n − 2 i=1 j=1

een zuivere schatter is voor σ 2 (de maximum likelihood-schatter voor σ 2 is gelijk aan (m + n − 2 2)/(m+n)SX,Y , ga na). Als µ = ν dan bezit T een t-verdeling met m+n−2 vrijheidsgraden (zie de volgende alinea voor de afleiding). Net als bij de t-toets voor een steekproef in Voorbeeld 5.16 geldt ook hier dat het voldoende is om de verdeling van T in het randpunt µ = ν te beschouwen en dat de verdeling van T dan onafhankelijk is van σ 2 . We kiezen de kritieke waarde daarom gelijk aan tm+n−2,1−α0 , en de toets is: “Verwerp H0 als T ≥ tm+n−2,1−α0 ”. Om in te zien dat T een t-verdeling bezit, schrijven we de toetsingsgrootheid T als q 2 2 (X − Y )/ σm + σn . T = q 2 2 (m+n−2)SX,Y /σ m+n−2

De teller van deze uitdrukking bezit onder µ = ν de N (0, 1)-verdeling. Om de verdeling van de noemer te bepalen, merken we op dat de som van twee onafhankelijke chikwadraat-verdeelde stochastische grootheden weer een chikwadraat-verdeling bezit, met het aantal vrijheidsgraden gelijk aan de som van de aantallen vrijheidsgraden. Gebruikmakend van Stelling 5.15 zien we 2 dan dat (m + n − 2)SX,Y /σ 2 een χ2m+n−2 -verdeling bezit en onafhankelijk is van X − Y . De teller en de noemer van de toetsingsgrootheid zijn dus onafhankelijk. Dat T onder µ = ν de tm+n−2 -verdeling bezit, volgt nu uit de definitie van de t-verdeling. Toetsen voor het andere eenzijdige en het tweezijdige hypothesen, en de bijbehorende overschrijdingskansen kunnen worden afgeleid op analoge wijze als in het eensteekproefprobleem.

Men noemt de toets in Voorbeeld 5.20 de t-toets (of Student-toets) voor twee steekproeven. Deze toets wijkt essentieel af van de eensteekproeftoets voor verschillen uit Voorbeeld 5.19, omdat daar op natuurlijke wijze paren (Xi , Yi ) gedefinieerd waren, hetgeen hier niet het geval is. Als de co¨ ordinaten Xi en Yi binnen een paar (Xi , Yi ) afhankelijk zijn, mogen we de tweesteekproeven t-toets niet gebruiken, of in ieder geval niet met de t m+n−2,1−α0 kritieke waarde. Er is dan immers geen garantie dat de onbetrouwbaarheid kleiner dan of gelijk is aan α0 . Zijn echter zowel de paren (X1 , Y1 ), . . ., (Xn , Yn ) als de co¨ ordinaten Xi en Yi binnen ieder paar onderling onafhankelijk en normaal verdeeld, dan bezitten zowel de t-toets voor

42

5: Toetsen

paren (Voorbeeld 5.19) als de tweesteekproeven t-toets (Voorbeeld 5.20) onbetrouwbaarheid α0 en zijn beide toegestaan. De tweesteekproeven t-toets verdient dan echter de voorkeur vanwege zijn grotere onderscheidend vermogen. De intu¨ıtieve reden is dat de onbekende parameter σ 2 bij de t-toets voor paren geschat wordt met behulp van n onafhankelijke waarnemingen (de verschillen Zi = Xi − Yi , waarvan er n − 1 “vrij” zijn, ofwel met n − 1 vrijheidsgraden), 2 terwijl SX,Y is gebaseerd op 2n onafhankelijke waarnemingen (met 2n − 2 vrijheidsgraden). Dat tweede is natuurlijk beter. In Voorbeeld 5.20 hebben we aangenomen dat de variantie σ 2 voor beide steekproeven gelijk is, maar bij veel praktische problemen is dit onzeker of niet waar. Een algemenere probleemstelling wordt verkregen door aan te nemen dat de twee steekproeven uit de N (µ, σ 2 )- en N (ν, τ 2 )-verdeling afkomstig zijn. Gevraagd wordt dezelfde nulhypothese H0 : µ ≤ ν te toetsen, maar nu bij onbekende σ 2 en τ 2 . Dit is het befaamde Behrens-Fisher -probleem. Anders dan voor het probleem waarin σ 2 = τ 2 , waarin de zojuist besproken toets uniform meest onderscheidend is onder de zuivere toetsen, bestaat in de situatie van het Behrens-Fisher-probleem geen absoluut beste toets (vandaar: “probleem”). Wel zijn er verschillende redelijke toetsen (waarvoor we verwijzen naar de handboeken). Trekt men zich van de mogelijke ongelijkheid van σ 2 en τ 2 niets aan, maar past men gewoon de tweesteekproeven t-toets uit Voorbeeld 5.20 toe, dan kan de ware onbetrouwbaarheid van de toets flink afwijken van de gewenste onbetrouwbaarheid (die in dit verband de nominale onbetrouwbaarheid wordt genoemd). Tabel 5.1 geeft hiervan een indruk. Het effect van ongelijke varianties is relatief klein als m en n ongeveer gelijk zijn en niet te klein. (Men kan bewijzen dat de onbetrouwbaarheid naar α0 convergeert als m = n → ∞ voor iedere σ 2 /τ 2 !) Dit leidt tot het advies om zo mogelijk gelijke steekproefomvangen te kiezen. Dit is overigens ook verstandig als σ 2 = τ 2 , omdat het onderscheidend vermogen van de tweesteekproeven t-toets maximaal is als m = n (bij vaste m + n).

m 5 15 7

σ 2 /τ 2 n 3 5 7

0.2

0.5

1

2

3

0.100 0.180 0.063

0.072 0.098 0.058

0.050 0.050 0.050

0.038 0.025 0.058

0.031 0.008 0.063

Tabel 5.1. Ware onbetrouwbaarheid van tweezijdige tweesteekproeven t-toets bij ongelijke varianties en nominaal niveau 0.05.

Voorbeeld 5.21 (Asymptotische t-toets). Correcte toepassing van de tweesteekproeven ttoets uit Voorbeeld 5.20 veronderstelt dat de twee steekproeven normaal verdeeld zijn met gelijke varianties. Als de twee steekproeven allebei voldoend groot zijn, dan is noch de normaliteit, noch de veronderstelling van gelijke varianties essentieel, mits de toets als volgt wordt aangepast. Als toetsingsgrootheid kiezen we T =q

X −Y 2 SX m

+

2 SY n

.

Deze grootheid verschilt van de toetsingsgrootheid in Voorbeeld 5.20 door het gebruik van een andere schatter voor de standaarddeviatie in de noemer. Met behulp van de Centrale Limietstelling 8.21 kan men laten zien dat onder de hypothese µ = ν van gelijke verwachtingen voor de twee steekproeven, de grootheid T = T m,n in verdeling naar een standaard normale verdeling convergeert, als m, n → ∞, mits de varianties van de twee steekproeven bestaan en eindig zijn. Voor grote waarden van m en n kunnen we de nulhypothese H0 : µ ≤ ν daarom toetsen met behulp van de toets: “Verwerp H0 als T ≥ ξ1−α0 ”. De onbetrouwbaarheid van deze toets convergeert naar de onbetrouwbaarheidsdrempel α 0 als m, n → ∞, voor ieder paar van onderliggende verdelingen met eindige varianties. Voor niet te asymmetrische verdelingen is dit resultaat al bruikbaar voor m = n = 20. Het is lang niet altijd redelijk om aan te nemen dat de data afkomstig zijn uit normale verdelingen. Als goede redenen bestaan voor een ander parametrisch model, bijvoorbeeld exponentiële verdelingen, dan zal dit in het algemeen tot een andere toets aanleiding geven, omdat


43

de t-toets in dat geval niet het juiste niveau en mogelijk ook een onnodig klein onderscheidend vermogen bezit. De algemene methoden voor toetsconstructie, zoals de likelihood-ratiotoets suggereren welke toets redelijk is. Het is ook mogelijk correcte toetsen te vinden die heel weinig aannames over de verdeling van de data vereisen. Zogenaamde verdelingsvrije toetsen werken voor heel ruime klassen van verdelingen. De tekentoets uit Voorbeeld 5.17 behoort tot deze groep. We bespreken hieronder een voorbeeld van een verdelingsvrije tweesteekproeventoets. Voorbeeld 5.22 (Wilcoxon). Bij gegeven steekproeven X1 , . . ., Xm , Y1 , . . ., Yn definiëren we de rangnummers R1 , . . ., Rm van de eerste steekproef in de totale steekproef als de positienummers van X1 , . . ., Xm na ordening op grootte van X1 , . . ., Xm , Y1 , . . ., Yn . (Bijvoorbeeld, als X1 de op drie na kleinste is onder alle waarnemingen, dan definiëren we R1 = 4; als X2 de P grootste m is, dan is R2 = m + n, etc.) De toetsingsgrootheid van de Wilcoxon -toets is W = i=1 Ri . Grote waarden van W geven een indicatie dat X1 , . . ., Xm relatief groot zijn ten opzichte van Y1 , . . ., Yn . Dit leidt tot het verwerpen van de nulhypothese H0 dat de twee steekproeven identiek verdeeld zijn tegen het alternatief dat de eerste steekproef uit een “stochastisch grotere verdeling” komt voor grote waarden van W . Uiteraard kunnen we ook eenzijdig naar de andere kant en tweezijdig toetsen. Onder de nulhypothese zijn X1 , . . ., Xm , Y1 , . . ., Yn te beschouwen als een steekproef ter grootte m + n uit een vaste (onbekende) verdeling. De rangnummers R1 , . . ., Rm kunnen dan worden beschouwd als een willekeurige greep van m getallen uit de getallen {1, 2, . . ., m + n}. (We veronderstellen gemakshalve dat de waarnemingen continu verdeeld zijn, zodat de rangnummers altijd eenduidig zijn bepaald.) De verdeling van de Wilcoxon-grootheid onder de nulhypothese is daarom onafhankelijk van deze verdeling, en kan worden bepaald op grond van combinatorische argumenten. Deze verdeling is getabelleerd en via statistische pakketten op de computer beschikbaar.

5.5.4

Aanpassingstoetsen

Een toets om vast te stellen of de ware verdeling van de waarneming tot een bepaalde klasse verdelingen behoort wordt een toets voor aanpassing (Engels: goodness of fit) genoemd. Eigenlijk past deze categorie toetsen slecht in de algemene filosofie van het toetsen, omdat bij aanpassingstoetsen men de nulhypothese meestal liever niet verwerpt. De nulhypothese zegt bijvoorbeeld dat de data opgevat kunnen worden als een steekproef uit een normale verdeling, en het zou ons de meeste informatie verschaffen als we deze nulhypothese zouden kunnen aantonen. De algemene opzet van de toetsingstheorie geeft ons deze mogelijkheid echter niet: de enig mogelijke sterke conclusie is dat de nulhypothese onjuist is; in het andere geval houden we ons op de vlakte. Men zou kunnen denken dat het omdraaien van de nul- en de alternatieve hypothese het probleem oplost. Immers, wanneer we dan de nulhypothese verwerpen, zouden we de sterke conclusie hebben dat de data uit een normale verdeling afkomstig zijn. Echter, deze nulhypothese zal in de praktijk nooit verworpen kunnen worden. De nulhypothese bevat in dat geval alle niet-normale verdelingen. Iedere normale verdeling in de alternatieve hypothese kan willekeurig dicht benaderd worden door een niet-normale verdeling in de nulhypothese. Het is daardoor onmogelijk een duidelijk onderscheid te maken tussen de nul- en de alternatieve hypothese en de sterke conclusie te trekken. Daarom kiezen we voor de eerder genoemde nulhypothese dat de verdeling waaruit de waarnemingen afkomstig zijn een normale verdeling is. In overeenstemming met deze handelswijze is het verstandig de resultaten van aanpassingstoetsen pragmatisch te interpreteren. Wordt bijvoorbeeld de nulhypothese van normaliteit niet verworpen, dan nemen we dit als indicatie dat het gebruik van de normale verdeling niet onredelijk is, zonder dat we het als afdoend bewijs van normaliteit opvatten. Het is simpelweg onmogelijk de juistheid van een bepaalde verdeling aan te tonen. Voorbeeld 5.23 (Kolmogorov-Smirnov). Veronderstel dat de waarnemingen X 1 , . . ., Xn een steekproef zijn uit een onbekende verdeling F , en dat we de nulhypothese H0 : F = F0 dat deze

44

5: Toetsen

gelijk is aan een gegeven verdeling F0 willen toetsen tegen het alternatief H1 : F 6= F0 dat dit niet zo is. De verdeling F0 zou bijvoorbeeld de standaard normale verdeling kunnen zijn. De Kolmogorov-Smirnov-toets is gebaseerd op de empirische verdelingsfunctie F n van X1 , . . ., Xn , welke is gedefinieerd als n

Fn (x) =

1X 1 #(Xi ≤ x) = 1{Xi ≤x} . n n i=1

Fn (x) is gelijk aan het aantal waarnemingen dat kleiner dan of gelijk is aan x, gedeeld door n. P Wegens de Wet van de Grote Aantallen geldt dat Fn (x) → E1{X≤x} = F (x) als n → ∞. Voor niet te kleine waarden van n moet Fn daarom dicht bij de echte verdelingsfunctie liggen, dus bij F0 als H0 juist is. De Kolmogorov-Smirnov-statistiek is de maximale afstand tussen Fn en F0 , T = sup Fn (x) − F0 (x) . x∈R

We verwerpen H0 : F = F0 voor grote waarde van T . De kritieke waarde voor de toets kunnen we afleiden uit de kansverdeling van T onder H0 . Deze heeft geen bijzondere naam, maar is wel getabelleerd en in statistische pakketten op de computer beschikbaar. Handig hierbij is dat de verdeling hetzelfde is voor iedere continue verdelingsfunctie F0 , zodat één tabel volstaat. Voor grote waarden van n kunnen we ook gebruik maken van het limietresultaat ∞ X √ 2 2 (−1)j+1 e−j z . lim PF0 sup Fn (x) − F0 (x) > z/ n = 2

n→∞

x∈R

j=1

0.0

0.2

0.4

0.6

0.8

1.0

De reeks aan de rechterkant kan op eenvoudige wijze numeriek berekend worden voor gegeven z. Daardoor is bovenstaande gelijkheid met name handig om overschrijdingskansen te bepalen.

-3

-2

-1

0

1

2

3

Figuur 5.10. De empirische verdelingsfunctie van een steekproef ter grootte 25 uit de N (0, 1)-verdeling en de ware verdelingsfunctie.

In veel gevallen in de praktijk is het zojuist besproken probleem te eenvoudig. We willen vaak nieteen enkelvoudige nulhypothese H0 : F = F0 toetsen, maar een hypothese van de vorm H0 : F ∈ Fθ : θ ∈ Θ voor een gegeven statistisch model {Fθ : θ ∈ Θ . Bijvoorbeeld, toetsen of de waarnemingen “normaal verdeeld” zijn, correspondeert met de keuzes θ = (µ, σ 2 ) ∈ R × (0, ∞) en Fµ,σ2 = N (µ, σ 2 ). Een uitbreiding van de Kolmogorov-Smirnov-toetsingsgrootheid is T ∗ = sup Fn (x) − Fθˆ(x) , x∈R

voor een schatter θˆ van θ. We verwerpen weer voor grote waarden van T ∗ . Door de substitutie van θˆ is de verdeling van T ∗ echter niet gelijk aan die van T . In het algemeen hangt deze verdeling af van het te toetsen model, van welke schatter θˆ wordt gebruikt, en zelfs van de ware parameter θ. Voor sommige speciale gevallen is de verdeling getabelleerd. In andere gevallen gebruikt men benaderingen of bepaalt men kritieke waarden met behulp van computersimulatie.

5: Opgaven

45

Beschouw bijvoorbeeld de toepassing op het toetsen van normaliteit. Het ligt voor de hand de onbekende parameter θ = (µ, σ 2 ) te schatten met het steekproefgemiddelde en de steekproefvariantie. We verwerpen de nulhypothese van normaliteit voor grote waarden van de statistiek x − X (5.4) T ∗ = sup Fn (x) − Φ . S X x∈R

Om de kritieke waarde voor de toets, of een overschrijdingskans te bepalen, is het nodig de verdeling van deze statistiek te kennen onder de aanname dat de nulhypothese correct is. Hoewel de nulhypothese samengesteld is, kan men laten zien dat de verdeling van T ∗ hetzelfde is onder ieder element van de nulhypothese. Het is niet eenvoudig hier een analytische uitdrukking voor af te leiden, maar het is heel eenvoudig deze verdeling te benaderen door middel van simulatie. We simuleren een groot aantal keer, bijvoorbeeld 1000 keer, een steekproef uit de normale verdeling van dezelfde omvang als de data, en berekenen de waarde van de Kolmogorov-Smirnov-statistiek T ∗ voor ieder van de 1000 steekproeven. Een benadering voor de overschrijdingskans is dan de fractie van de 1000 waarden die groter zijn dan de waarde van de statistiek op de echte data. Voorbeeld 5.24 (Chikwadraat-toets). Veronderstel dat de waarnemingen X1 , . . ., Xn een steekproef zijn uit een onbekende verdeling F . Een alternatief voor de Kolmogorov-Smirnovtoets bij een enkelvoudige nulhypothese, H0 : F = F0 tegen H1 : F 6= F0 , is de chikwadraat-toets voor samenhang. Daarbij verdeelt men het bereik van X1 in een aantal aansluitende intervallen I1 , I2 , . . ., Ik . Het aantal waarnemingen in de steekproef in ieder interval, genoteerd als Nj voor j = 1, . . ., k, is de stochastische grootheid Nj = #(1 ≤ i ≤ n: xi ∈ Ij ). Onder de nulhypothese volgt de kans dat een waarneming Xi valt in het interval Ij , pj : = PH0 (X1 ∈ Ij ), uit de verdeling F0 voor j = 1, . . ., k en is de verwachting van het aantal waarnemingen in interval Ij gelijk aan npj . In de toetsingsgrootheid bij de chikwadraat-toets wordt de afwijking tussen het gerealiseerde aantal waarnemingen en het onder de nulhypothese verwachte aantal waarnemingen in de intervallen op een genormaliseerde manier gemeten: X2 =

k X (Nj − npj )2 . npj j=1

Onder de nulhypothese heeft X 2 voor vaste k bij benadering een chikwadraat-verdeling met k − 1 vrijheidsgraden. Deze benadering is betrouwbaar voor niet te kleine waarden van n. Als vuistregel hanteert men dat het verwachte aantal waarnemingen onder de nulhypothese in ieder interval, npj voor j = 1, . . ., k, minstens 5 moet zijn. Chikwadraat-toetsen komen ook in andere situaties voor. Het wezenlijke van een chikwadraat-toets is dat de toetsingsgrootheid X 2 een som is van termen van de vorm (Yi − EYi )2 /EYi waar Yi stochastische grootheden zijn.

Opgaven 1. Bij McRonald adverteert men met kwartpond hamburgers. De consumentenbond wil onderzoeken of het hier inderdaad kwartponders betreft. Men meet het gewicht van 100 als kwartpond hamburgers geafficheerde producten. Formuleer een statistisch model en beschrijf het toetsingsprobleem. 2. Een koffiebar trekt ’s ochtends voor 10 uur weinig klanten. Om meer klanten te trekken wordt overwogen de prijs van een kopje koffie voor 10 uur 50 eurocent te verlagen. Beschrijf een experiment om te meten of een dergelijke maatregel effect heeft. Geef het statistische model aan en beschrijf het toetsingsprobleem.

46

5: Toetsen

3. Traditioneel veronderstelt men een lineair verband y = α + βx1 + γx2 tussen de opbrengst y van een industrieel proces, de temperatuur x1 , en de hoeveelheid toegevoegde catalysator x2 . Een onderzoeker meent echter dat (binnen zekere grenzen) de temperatuur niet van invloed is op de opbrengst. Zijn collega gelooft daar niets van en wil met behulp van een statistische toets bewijzen dat de temperatuur wel degelijk een rol speelt. Beschrijf hoe deze vraagstelling past binnen het statistisch toetsen (o.a. statistisch model, hypothesen). 4. Volgens de verpakking behoort een pakje shag 50 gram tabak te bevatten. Om te kijken of de fabrikant voldoende shag in een pakje doet, is van 100 pakjes de inhoud gewogen. De gemiddelde inhoud blijkt 49,82 gram te zijn. Bekend is dat bij vulling de variantie gelijk is aan 1. Formuleer een statistisch model en beschrijf het toetsingsprobleem. Ga d.m.v. een geschikte toets na of de fabrikant aan de eis voldoet. Neem α0 = 0.05. 5. Zij X1 , . . ., X25 een steekproef uit de N (µ, 4)-verdeling. Men wenst de nulhypothese H0 : µ ≤ 0 te toetsen tegen H1 : µ > 0 bij onbetrouwbaarheidsdrempel α0 = 0.05. Het waargenomen steekproefgemiddelde is 0.63. (i) Bepaal het kritieke gebied van een geschikte toets. (ii) Dient H0 verworpen te worden? (iii) Bepaal het onderscheidend vermogen van de toets in µ = 1/2. (iv) Bepaal de overschrijdingskans bij deze toets. 6. Zij X1 , . . ., X100 onderling onafhankelijk N (µ, 25)-verdeelde stochastische grootheden. Men wenst de nulhypothese H0 : µ = 0 te toetsen tegen H1 : µ 6= 0 bij onbetrouwbaarheidsdrempel α0 = 0.05. Men vindt x = −1.67. (i) Ga door middel van een geschikte toets na of H0 verworpen dient te worden. (ii) Bepaal de overschrijdingskans. 7. Zij X1 , . . ., Xn een steekproef uit de N (µ, 4)-verdeling. Men wenst de nulhypothese H0 : µ ≤ 1 te toetsen tegen H1 : µ > 1 bij onbetrouwbaarheidsdrempel α0 = 0.05. Daar het van groot belang is in het onderhavige geval dat H0 daadwerkelijk wordt verworpen als µ = 2, wenst men n zo te kiezen dat bij de Gauss-toets de kans op een fout van de tweede soort bij µ = 2 ten hoogste 0.1 is. Hoe groot moet n tenminste zijn? 8. Stel dat X1 , . . ., Xn een steekproef is uit de N (µ, σ 2 )-verdeling met µ onbekend en σ 2 > 0 bekend. Beschouw het toetsingsprobleem H0 : µ ≤ µ0 tegen H1 : µ > µ0 waarbij µ0 een vast gegeven getal is. Veronderstel dat, in tegenstelling tot Voorbeeld 5.9, toetsingsgrootheid √ X wordt genomen. (i) Laat zien dat het kritieke gebied K = {(x1 , . . ., xn ): x ≥ ξ1−α0 σ/ n + µ0 } een toets geeft met onbetrouwbaarheid α0 . (ii) Laat zien dat het kritieke gebied√K uit het vorige onderdeel gelijk is aan het kritieke gebied op basis van toetsingsgrootheid n(X − µ0 )/σ dat wordt gegeven in Voorbeeld 5.9. 9. Zij X verdeeld volgens de bin(25, p)-verdeling. Men wenst H0 : p ≥ 0.4 te toetsen tegen H1 : p < 0.4. Als men in p = 0.3 een onderscheidend vermogen van tenminste 0.6 wil hebben, hoe groot moet men dan de onbetrouwbaarheid van de in aanmerking komende toets tenminste kiezen? Is dit bevredigend? 10. Op grond van twee onafhankelijke steekproeven X1 , . . ., X25 en Y1 , . . ., Y16 uit de N (µ, σ 2 ) respectievelijk N (ν, τ 2 )-verdeling wensen we te toetsen H0 : σ 2 ≥ 2τ 2 tegen H1 : σ 2 < 2τ 2 bij onbekende µ en ν en α0 = 0.01. (i) Wat is de conclusie als we als kwadraatsommen vinden: s2x = 46.7 en s2y = 45.1? (ii) Bepaal de bijbehorende overschrijdingskans. 11. Zij X1 , . . ., Xn een steekproef uit de N (µ, σ 2 )-verdeling, waarbij µ ∈ R en σ 2 > 0 onbekend zijn. 2 (i) Bewijs dat de toets “Verwerp H0 : σ 2 ≤ σ02 als (n − 1)SX /σ02 ≥ χ2n−1,1−α ” (beschreven in Voorbeeld 5.18) onbetrouwbaarheid α bezit. (ii) Het onderscheidend vermogen van deze toets is een functie van (µ, σ). Druk deze functie uit in de verdelingsfunctie van de chikwadraat-verdeling. (iii) Maak een schets van deze functie. 12. Zij X1 , . . ., Xn een steekproef uit de N (µ, σ 2 )-verdeling, waarbij µ bekend is. Hoe zou je de bekende waarde van µ kunnen gebruiken voor het construeren van een toets voor het toetsen van H 0 : σ 2 = σ02 tegen H0 : σ 2 6= σ02 ? Verwacht je dat deze toets een groter onderscheidend vermogen bezit dan de toets uit Voorbeeld 5.18?

5: Opgaven

47

13. Een chemisch proces behoort tenminste 800 ton chemicaliën te produceren per dag. De dagelijkse opbrengsten van een bepaalde week zijn 785, 805, 790, 793 en 802 ton. Geven deze gegevens aanleiding om te concluderen dat er iets mis is met het proces? Neem α0 = 0.05. Welke veronderstellingen zijn gemaakt? 14. In een experiment is de bloeddruk van 32 hypertensiepatiënten gemeten na inname van het bloeddrukverlagend medicijn Cozaar. In een tweede experiment is de bloeddruk van 20 hypertensiepatiënten gemeten na inname van Diovan, een ander bloeddrukverlagend medicijn. Noteer de bloeddrukwaardes in de twee experimenten als X1 , . . ., X32 en Y1 , . . ., Y20 . De gemeten uitkomsten zijn x = 163, y = 158, sX = 7.8 en sY = 9.0. (i) Ga door middel van een geschikte toets na of een van de beide medicijnen beter werkt. Gebruik een onbetrouwbaarheidsdrempel van 5 %. (ii) Bepaal de overschrijdingskans (bij benadering). 15. De bepaling van isolerende eigenschappen van olie kan geschieden door een glazen buisje, waarin zich 2 polen bevinden, met olie te vullen en vervolgens op de polen een spanningsverschil aan te brengen, dat men laat stijgen tot een vonk de isolatie doorbreekt. Men kan deze bepaling van de doorslagspanning zo vaak herhalen als men wil. In een door Youden en Cameron beschreven experiment voerde men steeds 2 bepalingen uit (“duplobepalingen”). Noemen we de doorslagspanningen bij de eerste bepaling X en bij de tweede bepaling Y , dan ligt het voor de hand te onderstellen dat X en Y dezelfde verdeling hebben (ook al is deze in de regel verschillend voor verschillende oliesoorten). Dit is echter geenszins zeker, daar een door de olie slaande vonk ionen kan achterlaten, die bij de tweede bepaling de uitkomst kan be¨ınvloeden. Men wenst na te gaan of een dergelijk effect inderdaad aanwezig is. In het experiment werden 10 oliemonsters (ieder van een ander soort olie) betrokken; van ieder monster werden twee vullingen onderzocht en bij iedere vulling twee bepalingen verricht. De uitkomsten zijn hieronder aangegeven: oliemonster 1 2 3 4 5 6 7 8 9 10

1e vulling 16 11 14 19 23 13 16 20 15 14

12 10 14 17 20 15 15 19 11 12

2e vulling 17 12 15 18 21 14 16 19 16 13

14 10 14 19 19 14 14 20 13 15

Toets de nulhypothese dat er geen systematisch verschil is tussen de duplobepalingen tegen de alternatieve hypothese, dat dit wel het geval is, bij onbetrouwbaarheidsdrempel α0 = 0.01 in de veronderstelling, dat alle doorslagspanningen onafhankelijk en normaal verdeeld zijn met dezelfde (onbekende) variantie. Geef aan hoe groot de overschrijdingskans ongeveer is.

6 Betrouwbaarheidsgebieden

6.1

Introductie

In Hoofdstuk 4 hebben we gezien hoe een parameter θ geschat kan worden met de waargenomen waarde t = T (x) van een schatter T . Binnen de context van het huidige hoofdstuk zullen we dergelijke schattingen ook wel aanduiden met puntschattingen. Een schatting t verschilt als regel van de te schatten θ. Met behulp van de in dit hoofdstuk te bespreken betrouwbaarheidsgebieden kan de mogelijke afwijking van de schatter T tot θ worden gekwantificeerd. Dit leidt in veel gevallen tot een intervalschatting L(x), R(x) , met de interpretatie dat θ met grote kans in dit interval ligt. Dergelijke intervalschattingen zijn bijvoorbeeld de basis van de onzekerheidsmarges die worden vermeld bij de resultaten van opiniepeilingen. De precieze definitie van een betrouwbaarheidsgebied is als volgt. Definitie 6.1. Zij X een variabele met een kansverdeling die van een parameter θ ∈ Θ afhangt. Een afbeelding X 7→ GX die als bereik de collectie deelverzamelingen van Θ heeft, is een betrouwbaarheidsgebied voor θ met onbetrouwbaarheid α als Pθ GX 3 θ ≥ 1 − α, voor alle θ ∈ Θ.

Met andere woorden, een betrouwbaarheidsgebied is een “stochastische deelverzameling” GX van Θ die met “grote kans” de ware parameter θ zal bevatten. Omdat van tevoren niet bekend is welke waarde van θ de ware waarde is, geldt de eis in de definitie voor alle mogelijke waarden van θ: onder aanname dat θ de ware waarde is, moet deze ware waarde met kans minstens 1 − α bevat zijn in GX . Nadat X = x waargenomen is, verandert de stochastische verzameling GX in een gewone, niet-stochastische, deelverzameling Gx van Θ. Gewoonlijk kiest men α klein, bijvoorbeeld α = 0.05, zodat de kans dat θ in het betrouwbaarheidsgebied ligt groot is. Naarmate we α kleiner kiezen, zal het betrouwbaarheidsgebied natuurlijk groter moeten zijn en dus minder informatie geven over θ, die dan echter wel “zekerder” is. Hier is weer sprake van een trade-off tussen twee gewenste doelen, zoals we die al zijn tegen gekomen bij het toetsen. Vaak wordt gezegd dat de kans dat de realisatie Gx de ware waarde θ zal bevatten minstens 1 − α is. Deze kansuitspraak is gemakkelijk verkeerd te interpreteren. In onze interpretatie is de ware waarde van θ vast; het gerealiseerde betrouwbaarheidsgebied G x is evenmin een toevalsvariabele. Derhalve is de ware θ bevat in het betrouwbaarheidsgebied G x , of niet. (Helaas weten we niet welke van de twee gevallen zich voordoet.) De kansuitspraak kan worden ge¨ınterpreteerd in de zin dat als we bijvoorbeeld 100 keer onafhankelijk het experiment dat aanleiding geeft tot X zouden uitvoeren, en 100 keer het betrouwbaarheidsgebied G x berekenen, dan mogen we verwachten dat (minstens) ongeveer 100(1 − α) van de gebieden de ware θ bevatten. Dit wordt ge¨ıllustreerd in Figuur 6.1, waarin 100 onafhankelijke realisaties van een

6.1: Introductie

49

0

20

40

60

80

100

90 % betrouwbaarheidsinterval voor de verwachtingsparameter van de normale verdeling zijn weergegeven. De ware waarde van de parameter is 0 en is in 89 van de intervallen bevat. In de praktijk kunnen we natuurlijk niet herhalen, en kunnen we slechts één betrouwbaarheidsgebied bepalen. Dit kan één van de 100α gebieden zijn waar de ware parameter niet in zit, zonder dat we dit kunnen weten!

-2

-1

0

1

2

3

Figuur 6.1. 100 realisaties van het betrouwbaarheidsinterval voor de verwachtingsparameter van de normale verdeling (als in Voorbeeld 6.4) gebaseerd op 100 onafhankelijke steekproeven ter grootte 5.

Omdat GX stochastisch is en θ deterministisch hebben we GX 3 θ geschreven in plaats van θ ∈ GX . In onze notatie voor kansen staat de stochast immers meestal links. Sommigen vinden om dezelfde reden ook een uitspraak als “θ ligt met grote kans in GX ” uit den boze. We volgen de laatste conventie niet, maar benadrukken nogmaals dat betrouwbaarheidsgebieden een subtiele interpretatie bezitten. Als θ een numerieke parameter is (dwz. Θ ⊂ R), dan gebruiken we meestal betrouwbaar heidsintervallen. Dit zijn betrouwbaarheidsgebieden van de vorm GX = L(X), R(X) voor twee functies L en R van X. We spreken dan ook wel van het betrouwbaarheidsinterval [L, R] voor de parameter θ. Soms is het midden van het betrouwbaarheidsinterval precies de gebruikte puntschatter T = T (X) voor θ. Dan noteren we het interval ook wel in de vorm θ = T ± η, met η = 12 R(X) − L(X) de helft van de lengte van het interval. In andere gevallen is het interval echter bewust asymmetrisch rond de gebruikte puntschatting, hetgeen een uitdrukking kan zijn van een “grotere precisie” naar boven of beneden. Voorbeeld 6.2 (Normale verdeling). Veronderstel dat X = (X1 , . . ., Xn ) een steekproef is uit de normale verdeling N (µ, σ 2 ) met onbekende µ ∈ R en bekende variantie σ 2 . De natuurlijke schatter voor µ is het steekproefgemiddelde X. Deze bezit een N (µ, σ 2 /n)-verdeling en dus is √ n(X − µ)/σ standaard normaal verdeeld. Dan geldt √ X −µ ≤ ξ1−α/2 = 1 − α, Pµ ξα/2 ≤ n σ

waar ξα het α-kwantiel van de standaard normale verdeling is. We kunnen dit herschrijven in de vorm σ σ Pµ X − √ ξ1−α/2 ≤ µ ≤ X + √ ξ1−α/2 = 1 − α n n

waar we gebruik hebben gemaakt van ξα/2 = −ξ1−α/2 . Hieruit volgt dat Pµ (GX 3 µ) = 1 − α voor σ σ GX = X − √ ξ1−α/2 , X + √ ξ1−α/2 n n

50

6: Betrouwbaarheidsgebieden

en dat GX een betrouwbaarheidsinterval is voor µ met onbetrouwbaarheid α. Dit interval ligt symmetrisch rond de schatter X en wordt vaak geschreven als σ µ = X ± √ ξ1−α/2 . n Hoe kleiner σ en hoe groter n, des te smaller (en dus informatiever) het interval. Merk op dat voor het halveren van de intervallengte 4 maal zoveel waarnemingen nodig zijn. Ook bij grotere α is het interval smaller, maar dit gaat ten koste van de betrouwbaarheid.

6.2

Pivots en Bijna-Pivots

Veel betrouwbaarheidsgebieden zijn geconstrueerd met behulp van een pivot. Definitie 6.3. Een pivot is een functie (X, θ) 7→ T (X, θ), zodanig dat de kansverdeling van T (X, θ) onder de aanname dat θ de ware parameter is een vaste verdeling bezit, die niet afhangt van θ of andere onbekende parameters. Een pivot is dus geen statistiek omdat de pivot af mag hangen van zowel de waarneming X als de parameter θ. Wanneer T (X, θ) een pivot is, is de kans Pθ T (X, θ) ∈ B in principe bekend voor iedere verzameling B.“Bekend” betekent hier “onafhankelijk van θ”; de twee θ’s in de uitdrukking Pθ T (X, θ) ∈ B moeten elkaar dus opheffen. In Voorbeeld 6.2 hebben we √ al een voorbeeld van een pivot gezien: n(X − µ)/σ die de standaard normale verdeling bezit. Voor iedere verzameling B zodanig dat Pθ T (X, θ) ∈ B ≥ 1 − α is de verzameling n o θ ∈ Θ: T (X, θ) ∈ B een betrouwbaarheidsgebied voor θ met onbetrouwbaarheid α. Meestal bestaan vele verzamelingen B die hieraan voldoen, en we willen daar nu een “geschikte” kandidaat uit kiezen. Hoewel het voor de hand ligt te zoeken naar verzamelingen waarvoor het volume van het betrouwbaarheidsgebied klein is, ligt de keuze niet eenduidig vast. We illustreren dit in de volgende voorbeelden.

Voorbeeld 6.4 (Normale verdeling). Veronderstel dat X1 , . . ., Xn een steekproef is uit de N (µ, σ 2 )-verdeling met µ ∈ R en σ 2 > 0 onbekend. Volgens Stelling 5.15 bezit √ X −µ n SX een tn−1 -verdeling, welke niet afhangt van de parameter (µ, σ 2 ). Deze grootheid is dus een pivot en er geldt √ X −µ Pµ tn−1,α/2 ≤ n ≤ tn−1,1−α/2 = 1 − α. SX Uit berekeningen analoog aan die in Voorbeeld 6.2 volgt direct dat SX SX X − √ tn−1,1−α/2 , X + √ tn−1,1−α/2 n n

een betrouwbaarheidsinterval voor µ is met onbetrouwbaarheid α. Omdat het interval symmetrisch rond X ligt, wordt het ook wel geschreven als SX µ = X ± √ tn−1,1−α/2 . n Dit interval lijkt sterk op het interval in het vorige voorbeeld, met σ vervangen door S X en ξα door tn−1,α . Omdat de t-verdeling dikkere staarten heeft dan de standaard normale verdeling,

6.2: Pivots en Bijna-Pivots

51

liggen de t-kwantielen verder van nul vandaan dan de kwantielen van de standaard normale verdeling, en is het hier gevonden interval gewoonlijk iets wijder dan in het geval dat σ bekend is (hoewel dat ook van de waarde van SX afhangt). Dit is de prijs die we voor het onbekend zijn van σ moeten betalen. Voor n → ∞ lijkt de tn -verdeling in toenemende mate op de normale verdeling en geldt dat SX in kans naar σ convergeert. Daarom verdwijnt het verschil tussen de twee intervallen voor n → ∞. Door de keuze van de kwantielen is bovenstaand interval symmetrisch rond de maximum likelihood-schatter voor µ. Niet-symmetrische intervallen met onbetrouwbaarheid α kunnen worden geconstrueerd door andere kwantielen van de t-verdeling te kiezen: √ X −µ Pµ tn−1,β ≤ n ≤ tn−1,1−γ = 1 − α, SX voor β + γ = α. Het betrouwbaarheidsinterval voor µ op basis van deze kwantielen is gelijk aan SX SX X − √ tn−1,1−γ , X − √ tn−1,β . n n

Het smalste betrouwbaarheidsinterval met onbetrouwbaarheid α wordt verkregen door β = γ = α/2 te nemen; dit resulteert in het eerder gegeven interval. Exacte betrouwbaarheidsgebieden afleiden uit pivots is slechts incidenteel mogelijk, eenvoudig weg omdat niet altijd een pivot beschikbaar is. Het lukt bijvoorbeeld niet voor de parameter p in de binomiale verdeling, of de parameter µ in de Poisson-verdeling. In zo’n geval wordt vaak genoegen genomen met een benaderend betrouwbaarheidsgebied, dat kan worden afgeleid uit een “bijna-pivot”. Hebben we te maken met grote steekproeven, dan zijn zulke bijna-pivots meestal ruim voorradig. Voorbeeld 6.5 (Binomiale verdeling). Als X binomiaal verdeeld is met parameters n en p, dan is X − np p np(1 − p)

voor grote n bij benadering N (0, 1)-verdeeld, vanwege de Centrale Limietstelling, zie Paragraaf 8.6. Bij benadering is deze functie van X en p dus een pivot. De verzameling n

o X − np p: ξα/2 ≤ p ≤ ξ1−α/2 np(1 − p)

0

2

4

6

8

10

12

is derhalve bij benadering een betrouwbaarheidsgebied voor p met onbetrouwbaarheid α. Deze verzameling is een interval, dat gevonden kan worden door de kwadratische vergelijking (X − 2 np)2 ≤ ξ1−α/2 np(1 − p) op te lossen. In Figuur 6.2 wordt dit interval grafisch weergegeven voor bepaalde waarden van α, n en p.

0.0

0.2

0.4

0.6

0.8

1.0

Figuur 6.2. Betrouwbaarheidsinterval voor de parameter p van een binomiale verdeling. De grafiek toont de p functies p 7→ |x − np| en p 7→ 1.96 np(1 − p) voor een geval dat 0 < x < n (namelijk x = 13 en n = 20). Het betrouwbaarheidsinterval is het interval op de horizontale as tussen de twee snijpunten.

52


Als we toch aan het benaderen zijn, kunnen we ook een stap verder gaan. Vanwege de Wet van de Grote Aantallen convergeert X/n in kans naar p voor n naar oneindig. Daarom is ook de stochastische grootheid X − np p n(X/n)(1 − X/n)

bij benadering N (0, 1)-verdeeld. (Dit volgt met behulp van Slutzky’s lemma. We geven geen precies bewijs, maar verwijzen naar Hoofdstuk 2 van [vdV].) Het benaderende betrouwbaarheidsinterval op basis van deze bijna-pivot heeft de eenvoudige vorm r X X X 1 p= 1− ξ1−α/2 . ±√ n n n n

Dit interval wordt veel gebruikt als betrouwbaarheidsindicatie bij het schatten van een fractie elementen in een populatie met een bepaalde eigenschap, bijvoorbeeld een opiniepeiling. Het is opmerkelijk dat de grootte van de populatie geen rol speelt in de lengte van het interval. Alleen de steekproefgrootte telt, en in mindere mate de ware fractie. Als p = 12 en n = 1500, dan is het 95% betrouwbaarheidsinterval ongeveer (X/n) ± 2%. Deze 2% is vermoedelijk de waarde die wordt bedoeld als in de krant een afwijking van hoogstens 2% wordt beloofd in een gegeven opiniepeiling. De correcte interpretatie van deze marge is dat in 95% van de opiniepeilingen de afwijking van de steekproeffractie tot de ware fractie niet groter zal zijn dan 2%. Helaas vertaalt de Nederlandse pers deze ingewikkelde bewering vaak in een te stellige foutenmarge. p Voor p → 0 of p → 1 nadert de functie p 7→ p(1 − p) de waarde 0. Het betrouwbaarheidsinterval is daarom korter voor extreme waarden van p. Voor p = 1/2 is de lengte van het betrouwbaarheidsinterval het minst gunstig. De bijna-pivot in Voorbeeld 6.5 komt van een asymptotische benadering voor de verdeling van de schatter. Veel schatters Tn voor een parameter g(θ) zijn asymptotisch normaal verdeeld in de zin dat voor zekere getallen σn,θ (vaak de standaardafwijking van Tn ), Tn − g(θ) σn,θ

θ

N (0, 1)

als n → ∞. Hierin is het pijltje θ onze notatie voor “convergentie in verdeling” aangenomen dat θ de ware parameter is. Preciezer gezegd betekent de bewering dat T − g(θ) n ≤ x = Φ(x), voor alle x. lim Pθ n→∞ σn,θ Een informele interpretatie is dat Tn −g(θ) /σn,θ voor grote n bij benadering N (0, 1)-verdeeld is, als θ de ware parameter is. Derhalve is Tn − g(θ) σn,θ een bijna-pivot. Dit noemen we ook wel de grote steekproefmethode. Dit leidt tot het benaderende betrouwbaarheidsgebied voor g(θ) n o g(θ): Tn − σn,θ ξ1−α/2 ≤ g(θ) ≤ Tn + σn,θ ξ1−α/2

met onbetrouwbaarheid α. Voor het gemak wordt hierin de uitdrukking σn,θ vaak vervangen door een schatter σ ˆn , hetgeen leidt tot het symmetrische interval g(θ) = Tn ± σ ˆn ξ1−α/2 . De uitdrukking σ ˆn is meestal een schatting voor de standaardafwijking van Tn , de standaardfout (Engels: standard error of s.e.) van de schatter (of schatting). In veel wetenschappelijke rapportages wordt volstaan met het vermelden van de schatting met een bijbehorende standaardfout. Mits de gebruikte schatter bij benadering normaal verdeeld is, kunnen we deze

6.3: Maximum Likelihood-Schatters als Bijna-Pivots

53

informatie ruwweg interpreteren in de zin van een 95 % betrouwbaarheidsinterval van de vorm g(θ) = Tn ± ξ0.975 s.e. = Tn ± 1.96 s.e.. Goede statistische software zal naast een parameterschatting ook de standaardfout van de schatter vermelden. Voor een schatting van een vectorwaardige parameter levert dit voor iedere co¨ ordinaat een standaardfout, en worden bovendien meestal de geschatte covarianties tussen de schatters vermeld, in de vorm van een matrix met de geschatte varianties van de schatters (de kwadraten van de standaardfouten) op de diagonaal.

6.3

Maximum Likelihood-Schatters als Bijna-Pivots

Een belangrijk speciaal geval van de in de voorgaande paragraaf besproken bijna-pivots is dat waarin Tn de maximum likelihood-schatter is. Onder bepaalde voorwaarden is de maximum likelihood-schatter asymptotisch normaal verdeeld. We bespreken het eenvoudigste geval, dat van een steekproef van onderling onafhankelijke stochasten, en beperken ons eerst tot parameters θ ∈ R. Zij θˆ de maximum likelihood-schatter gebaseerd op een steekproef X1 , . . ., Xn uit de verdeling met (marginale) kansdichtheid pθ . Veronderstel dat de functie θ 7→ `θ (x): = log pθ (x) (partieel) differentieerbaar is voor alle x, met afgeleide ∂ `˙θ (x) = log pθ (x). ∂θ Pn De maximum likelihood-schatter voldoet dan (meestal) aan de likelihood-vergelijking i=1 `˙θˆ(Xi ) = 0. De functie `˙θ heet de score-functie van het model en de Fisher-informatie is gedefinieerd als het getal iθ = varθ `˙θ (X1 ). Stelling 6.6. Veronderstel dat de afbeelding θ 7→ pθ (x) differentieerbaar is voor alle x en dat iθ eindig is. Zij θˆn de maximum likelihood-schatter gebaseerd op een steekproef ter grootte n uit de √verdeling met kansdichtheid pθ . Dan geldt onder bepaalde voorwaarden dat, onder θ, de rij n(θˆn − θ) in verdeling naar een normale verdeling met verwachting 0 en variantie i −1 θ convergeert. Dus √ θ n(θˆn − θ) N (0, i−1 ). θ

Voor een precieze formulering en een bewijs van deze stelling verwijzen we naar [vdV]. Indien de uitspraak van de stelling van toepassing is, dan is voor grote n de stochastische grootheid √ niθ (θˆ − θ) onder θ bij benadering standaard normaal verdeeld, en derhalve een bijna-pivot. Voor het gemak kunnen we iθ vervangen door een schatter ibθ en krijgen we voor θ het benaderende betrouwbaarheidsinterval 1 θ = θˆ ± p ξ1−α/2 nibθ

met onbetrouwbaarheid α. Dit interval wordt het Wald-interval genoemd. De uitspraak van ˆ en de stelling wordt vaak zo begrepen dat 1/(niθ ) een benadering is voor de variantie van θ, de wortel hieruit een benadering voor de standaardfout. Voor α = 0.05 is het Wald-interval daarom in feite een interval van de algemene vorm θ = θˆ±ξ0.975 s.e. ≈ θˆ±2s.e.. (De stelling zegt overigens niets over convergentie van de variantie van de MLE schatters, maar de voorgaande interpretatie is meestal wel te verdedigen.) Gebruikelijke schatters voor de Fisher-informatie iθ zijn de plug-in schatter en de waargenomen informatie (Engels: observed information). De plug-in schatter is ibθ = iθˆ, ofwel

54


de parameter θ in de uitdrukking voor iθ wordt vervangen door de maximum likelihood-schatter ˆ De waargenomen informatie is gedefinieerd als θ. n 1 X¨ ibθ = − ` ˆ(Xi ), n i=1 θ

met

∂2 `¨θ (x) = 2 log pθ (x). ∂θ

De plug-in schatter vereist de (analytische) berekening van de Fisher-informatie i θ , terwijl de waargenomen informatie eenvoudiger uit de data volgt. De waargenomen informatie is (−1/n) Pn ˆ Zo keer de tweede afgeleide van de log likelihood-functie θ 7→ i=1 `θ (Xi ) geëvalueerd in θ = θ. nodig kan men een numerieke afgeleide (differentiequotiënt) gebruiken in plaats van een analytische afgeleide. Grafisch geeft de waargenomen informatie de kromming van de log likelihoodfunctie in het punt θ = θˆ waar de log likelihood maximaal is. Als de likelihood-functie een platte top bezit, dan is de waargenomen informatie klein, en is het betrouwbaarheidsinterval voor θ breed: de maximum likelihood-schatter is dan weinig nauwkeurig. (Dit reflecteert geen zwakte van deze schattingsmethode, maar is te wijten aan een intrinsiek moeilijk te schatten parameter.) De zinvolheid van de waargenomen informatie als schatter voor i θ is niet onmiddellijk duidelijk, maar volgt (voor Pn grote n) uit het volgende lemma, en de Wet van de Grote Aantallen, volgens welke n−1 i=1 `¨θ (Xi ) → Eθ `¨θ (Xi ) als n → ∞ met kans 1 als θ de ware parameter is. Lemma 6.7. Veronderstel dat θ 7→ `θ (x) = log pθ (x) twee maal differentieerbaar is voor alle x. Dan geldt onder regulariteitsvoorwaarden dat Eθ `˙θ (X1 ) = 0 en Eθ `¨θ (X1 ) = −iθ . Bewijs. We schrijven de formules onder de aanname dat X1 continu verdeeld is. (Voor een discrete kansdichtheid vervangen we de integralen door sommen.) Omdat pθ een kansdichtheid R is, geldt dat 1 = pθ (x) dx voor alle θ. Derhalve geldt Z Z Z ∂ ∂ 0= pθ (x) dx = pθ (x) dx = p˙ θ (x) dx, ∂θ ∂θ met p˙ θ (x) = ∂/∂θ pθ (x). De verwisseling van differentiatie (naar θ) en integraal (naar x) is toegestaan onder “regulariteitsvoorwaarden”. Aangezien `˙θ (x) = ∂/∂θ log pθ (x) = p˙ θ (x)/pθ (x) kunnen we de rechterkant herschrijven als Z Z p˙θ (x) pθ (x) dx = `˙θ (x) pθ (x) dx = Eθ `˙θ (X1 ). pθ (x) Dit voltooit het bewijs vanRde eerste bewering: Eθ `˙θ (X1 ) = 0. Voor het bewijs van de tweede bewering differentiëren we pθ (x)dx tweemaal naar θ en vinden we ∂2 0= 2 ∂θ

Z

pθ (x) dx =

Z

p¨θ (x) dx,

met p¨θ (x) = ∂ 2 /∂θ2 pθ (x). Differentiatie van de gelijkheid `˙θ (x) = p˙ θ (x)/pθ (x) naar θ geeft p¨θ (x) p˙ θ (x) 2 p¨θ (x) `¨θ (x) = = − − `˙θ (x)2 . pθ (x) pθ (x) pθ (x) We vermenigvuldigen dit met pθ (x) en nemen de integraal met betrekking tot x om te vinden dat Z Z Eθ `¨θ (X1 ) = p¨θ (x) dx − `˙θ (x)2 pθ (x) dx

= 0 − Eθ `˙θ (X1 )2 = − varθ `˙θ (X1 ) = −iθ , aangezien varθ `˙θ (X1 ) = Eθ `˙θ (X1 )2 − (Eθ `˙θ (X1 ))2 = Eθ `˙θ (X1 )2 , vanwege de eerste bewering. Dit bewijst de tweede bewering.

6.3: Maximum Likelihood-Schatters als Bijna-Pivots

55

Voorbeeld 6.8 (Poisson-verdeling). Zij X1 , . . ., Xn een steekproef uit de Poisson(θ)verdeling, waarbij θ > 0 onbekend is. De maximum likelihood-schatter voor θ is θˆ = X mits X > 0. De score-functie is gelijk aan e−θ θx x ∂ log = − 1. `˙θ (x) = ∂θ x! θ De Fisher-informatie is dan

1 −1 = . θ θ Volgens het vorige lemma hadden we dezelfde uitdrukking gekregen met de vergelijking iθ = varθ

X

1

iθ = −Eθ `¨θ (X1 ) =

E θ X1 1 = . θ2 θ

Indien we θ schatten met X, dan is de plug-in schatter voor iθ gelijk aan 1/X. De waargenomen informatie geeft dezelfde schatter, aangezien n

−

n

X 1 X Xi 1 1 X¨ = `θˆ(Xi ) = = . 2 n i=1 n i=1 θˆ2 (X) X

Het symmetrische benaderende betrouwbaarheidsinterval met onbetrouwbaarheid α wordt dan √ X θ = X ± √ ξ1−α/2 . n Dit interval hadden we ook via een√meer directe √ weg kunnen vinden, door toepassing van de Centrale Limietstelling op X. De rij n(X −θ)/ θ is immers bij benadering standaard normaal verdeeld (zie Voorbeeld 8.23). Voorbeeld 6.9 (Exponenti¨ ele verdeling). Zij X1 , . . ., Xn een steekproef uit de exponenti¨ ele ˆ = 1/X (zie verdeling met onbekende parameter λ. De maximum likelihood-schatter voor λ is λ Voorbeeld 4.7). De score-functie is gelijk aan ∂ 1 log λe−λx = − x `˙λ (x) = ∂λ λ en de Fisher-informatie is

1

1 − X1 = 2 . λ λ Volgens Lemma 6.7 is de Fisher-informatie ook te vinden middels de vergelijking iλ = −Eλ `¨λ (X1 ) = 1/λ2 . Als λ wordt geschat met de maximum likelihood-schatter, dan is de plug-in schatter voor iλ gelijk aan (X)2 . De waargenomen informatie geeft dezelfde schatter voor iλ : n n 1 X¨ 1X 1 2 − `λˆ (Xi ) = =X . 2 ˆ n n λ iλ = varλ

i=1

i=1

Voor beide schatters van iλ vinden we dan het symmetrische benaderende betrouwbaarheidsinterval voor λ 1 1 λ= ±√ ξ1−α/2 X nX

met onbetrouwbaarheid α.

56


Opgaven 1. In een laboratorium tracht men een bepaalde grootheid θ te meten. Bij deze metingen treden normaal verdeelde meetfouten op met bekende standaardafwijking 2.3 en met verwachting 0. Men voert 25 onafhankelijke metingen uit en vindt de gemiddelde waarde 18,61. Bepaal een (numeriek) betrouwbaarheidsinterval voor θ met een onbetrouwbaarheid 0.01. 2. Als in de voorgaande opgave de standaardafwijking niet bekend mag worden verondersteld, en de steekproefstandaardafwijking S bedraagt 2.3, hoe is dan het (numeriek) betrouwbaarheidsinterval voor θ met onbetrouwbaarheid 0.01? 3. Veronderstel dat X1 , . . ., Xm en Y1 , . . ., Yn onafhankelijke aselecte steekproeven zijn uit een normale N (µ, σ 2 ) respectievelijk N (ν, σ 2 )-verdeling. Bepaal een betrouwbaarheidsinterval voor µ − ν met onbetrouwbaarheid α (i) als σ 2 bekend is; (ii) als σ 2 onbekend is. 4. Veronderstel dat X1 , . . ., Xn een steekproef is uit de exponentiële verdeling met parameter λ. (i) Bepaal een exact betrouwbaarheidsinterval voor λ gebaseerd op een geschikte pivot; (ii) Bepaal een benaderend betrouwbaarheidsinterval voor λ gebaseerd op de maximum likelihoodschatter en de grote steekproevenmethode. 5. Men voert 25 onafhankelijke Bernoulli-proeven uit, elk met onbekende kans p op succes. Men vindt 18 successen. Neem onbetrouwbaarheid 0.05. (i) Bereken een exact betrouwbaarheidsinterval voor p. (ii) Bereken een benaderend betrouwbaarheidsinterval voor p gebaseerd op de grote steekproefmethode. Is 25 in dit verband te beschouwen als “groot”? 6. Een fabrikant van weegschalen beweert dat door hem gefabriceerde weegschalen een nauwkeurigheid hebben van 2 promille. Dat betekent dat als X het gewicht voorstelt van een object van 1000 mg gemeten met een willekeurige weegschaal van de fabrikant, de variantie van X gelijk is aan 22 = 4 mg2 . We willen onderzoeken of de fabrikant gelijk heeft. Daartoe nemen we een object van 1000 mg en bepalen de massa van dit object met behulp van 100 willekeurige weegschalen van de fabrikant. De verschillende metingen worden genoteerd met X1 , . . ., X100 . We veronderstellen dat de waarnemingen X1 , . . ., X100 onderling onafhankelijk en normaal verdeeld zijn met verwachting µ en onbekende variantie σ 2 . De waargenomen steekproefvariantie is gelijk aan 4.8. (i) Construeer een betrouwbaarheidsinterval voor σ 2 met onbetrouwbaarheid 0.05 onder de aanname dat µ = 1000 mg bekend is. (ii) Construeer een betrouwbaarheidsinterval voor σ 2 met onbetrouwbaarheid 0.05 onder de aanname dat µ onbekend is. (iii) Beschrijf een toets om te toetsen of σ 2 significant afwijkt van de door de fabrikant opgegeven variantie. Doe dit zowel voor het geval dat µ bekend als onbekend is. Geef de nulhypothese. Deze toets mag je uitvoeren met behulp van het betrouwbaarheidsinterval in het vorige onderdeel, een kritiek gebied of een overschrijdingskans. Neem de onbetrouwbaarheidsdrempel gelijk aan 0.05.

7 Enkele Regressiemodellen

Het woord regressie heeft in het hedendaagse taalgebruik een negatieve lading, hoewel het in de statistiek de standaard aanduiding is voor het verklaren van een variabele Y met behulp van een variabele X. Een afhankelijke variabele Y wordt “teruggebracht” op een verklarende variabele X, die in deze context ook wel de onafhankelijke variabele wordt genoemd. Een greep uit de vele voorbeelden van toepassingen is: - het voorspellen van de opbrengst van een biochemisch proces als functie van temperatuur, hoeveelheid katalysator, etc. - het voorspellen van de prijs van onroerend goed als functie van grootte, ligging, aanwezige infrastructuur, etc. - het voorspellen van de respons op een mailing als functie van postcode, opleidingsniveau, inkomen, etc. - het voorspellen van het nationaal product uit macro-economische variabelen als beroepsbevolking, rentestand, begrotingstekort, inflatie, etc. - het voorspellen van studieduur als functie van eindexamencijfer, studieprofiel, etc. - het voorspellen van de eindlengte van een kind op basis van de lengte van de ouders en het geslacht van het kind Omdat de verklarende variabele X op vele wijzen invloed kan hebben op de afhankelijke variabele Y en de kansverdeling van Y niet voor elke toepassing dezelfde zal zijn, bestaan verschillende soorten regressiemodellen. Paragrafen 7.1, en 7.2 behandelen enkele veel gebruikte regressiemodellen, waarin de variabele Y een reëelwaardige stochastische grootheid is. Alle genoemde toepassingen hebben gemeen dat geen perfecte relatie tussen de variabelen X en Y bestaat, hoewel we wel een zeker verband verwachten. We verwachten bijvoorbeeld dat de grootte en de leeftijd van een pand, de ligging en mogelijk nog andere factoren, de prijs van het pand zullen be¨ınvloeden, maar het zal in het algemeen niet mogelijk zijn de prijs perfect te voorspellen uit een aantal van zulke indicatoren. Dit zou te wijten kunnen zijn aan een gebrek aan informatie (sommige relevante variabelen zijn nog onbekend), of aan toevallige factoren. In beide gevallen is het niet onredelijk om (x, y) te beschouwen als realisatie van een stochastische vector (X, Y ). Het verband tussen x en y kunnen we vervolgens onderzoeken middels de kansverdeling van de vector (X, Y ). Daarbij zijn we vooral ge¨ınteresseerd in de voorwaardelijke verdeling van Y gegeven X = x, en minder in de marginale verdeling van X. In sommige gevallen is de waarde van de verklarende variabele X onder controle. Bijvoorbeeld, bij het bepalen van optimale productieomstandigheden worden de verschillende instellingen x systematisch gevarieerd, waarna de opbrengsten y worden geanalyseerd. In een dergelijk geval is het niet redelijk de verklarende variabele als realisatie van een stochastische grootheid te beschouwen, maar beschouwen we alleen de afhankelijke variabele Y als een stochastische variabele. In het geval dat X wel stochastisch is, wordt meestal de voorwaardelijke verdeling van Y gegeven X = x gemodelleerd. Daarom maakt het voor de verschillende regressiemodellen praktisch niet veel uit of we X nu wel of niet stochastisch veronderstellen. In de volgende paragrafen zal worden aangegeven of X wel of niet stochastisch wordt verondersteld. We gaan steeds uit van beschikbare data (x1 , y1 ), . . ., (xn , yn ) als realisaties van ofwel de stochastische

58

7: Enkele Regressiemodellen

vector (X, Y ) ofwel als realisaties van de stochastische variabele Y in combinatie met gemeten niet-stochastische waarnemingen (x1 , . . ., xn ). Omdat causaliteit een belangrijk onderwerp is binnen de context van regressie, bespreken we twee korte voorbeelden waarin dit begrip wordt ge¨ıllustreerd. Het eerste voorbeeld gaat over het fabeltje dat baby’s gebracht worden door de ooievaar. In sommige gebieden is er inderdaad een positieve correlatie waargenomen in de schommelingen van de ooievaarspopulatie en het geboortecijfer; in periodes waarin het slecht ging met de ooievaarspopulatie daalde het geboortecijfer en op momenten dat de grootte van de ooievaarspopulatie steeg, nam ook het geboortecijfer toe. Dit is een opmerkelijk resultaat aangezien we er al lange tijd van overtuigd zijn dat baby’s niet door de ooievaar gebracht worden. Ondanks de correlatie is het niet te verwachten dat als in deze gebieden de grootte van de ooievaarspopulatie kunstmatig opgehoogd wordt, door bijvoorbeeld extra ooievaars uit te zetten, het geboortecijfer zal stijgen. Het tweede voorbeeld gaat over de positieve correlatie tussen inkomen en uitgaven: mensen die meer verdienen geven in het algemeen ook meer uit. In dit geval is het veelal wel zo dat als een persoon meer gaat verdienen zijn uitgaven ook omhoog gaan. Wat is het verschil tussen deze twee voorbeelden? In beide gevallen kunnen we zeggen dat de verklarende variabele (aantal ooievaars en inkomen) voorspellende waarde heeft voor de afhankelijke variabele (geboortecijfer en uitgaven). Echter, in het ooievaarsvoorbeeld kunnen we niet zeggen dat er een causaal verband is, terwijl dat in het inkomensvoorbeeld wel het geval is: het kunstmatig ophogen van het aantal ooievaars zal geen effect hebben op het geboortecijfer, terwijl zo’n effect waarschijnlijk wel waar te nemen is in het inkomensvoorbeeld. Waarom in sommige gebieden het geboortecijfer en de grootte van de ooievaarspopulatie positief gecorreleerd zijn, is niet geheel duidelijk. Mogelijk zijn beide afhankelijk van de industriële ontwikkeling: meer industrie betekent enerzijds meer rijkdom en dat veroorzaakt traditiegetrouw een afname in het geboortecijfer, en anderzijds toenemende luchtvervuiling en onrust in het gebied, waardoor ooievaars wegtrekken.

7.1

Lineaire Regressie

Het lineaire model, de basis voor lineaire regressie en variantie-analyse, is het werkpaard van de “klassieke” statistiek, enerzijds omdat het, met enig verstand van zaken, breed toepasbaar is, en anderzijds omdat de benodigde berekeningen op eenvoudige matrixalgebra zijn gebaseerd. Hoewel moderne rekenmiddelen de toepassing van veel flexibelere modellen mogelijk hebben gemaakt, is het lineaire model nog steeds van veel waarde. Het standaard lineaire regressiemodel veronderstelt dat, gegeven X = x = (x 1 , . . ., xp ), de variabele Y normaal verdeeld is, met conditionele verwachting en variantie E(Y |X = x) =

p X

βi xi ,

i=1

var(Y |X = x) = σ 2 ,

waarbij deze laatste onafhankelijk van x is. Wanneer X niet stochastisch is, heeft Y bij verklarende variabelen (x1 , . . ., xp ) (onvoorwaardelijk) bovenstaande verwachting en variantie. In de rest van deze paragraaf gaan we ervan uit dat de verklarende variabele X niet stochastisch is, zodat we de conditionering achterwege kunnen laten. Het model bevat p + 1 reëelwaardige parameters, die we in de parametervector θ = (β1 , . . ., βp , σ 2 ) samen kunnen nemen. In het lineaire regressiemodel beschikken we naast de meting Y ook over de verklarende variabele x, die we gebruiken om de verwachte waarde van Pp de afhankelijke variabele Y te modelleren. Als we een “meetfout” definiëren door e = Y − i=1 βi xi , dan kunnen we schrijven Y =

p X

βi xi + e.

i=1

In het standaard regressiemodel zijn de meetfouten onderling onafhankelijk en normaal verdeeld met verwachting 0 en variantie σ 2 . De variabele Y heeft dan eveneens een normale verdeling,

7.1: Lineaire Regressie

59

P maar nu met verwachting pi=1 βi xi en variantie σ 2 . We kunnen het lineaire regressiemodel beschouwen als een uitbreiding van het meetfoutenmodel van Voorbeeld 2.3. Het lineaire regressiemodel maakt een aantal specifieke keuzes: - de verwachting van Y is afhankelijk van x, maar de variantie niet. - de verwachting van Y is een lineaire functie van x. - de meetfout is normaal verdeeld. In een verrassend groot aantal toepassingen is aan deze keuzes voldaan, maar natuurlijk niet altijd. Hierbij moet worden bedacht dat de variabelen vaak door “preprocessing” in een vorm worden gebracht, die aansluit bij het lineaire regressiemodel. Zo kan de regressie worden uitgevoerd na transformatie van de afhankelijke variabele Y (bijvoorbeeld door het nemen van de logaritme, log Y ), en kunnen in het bijzonder de verklarende variabelen op vele manieren worden bewerkt. Wanneer bijvoorbeeld in het geval van een eendimensionale variabele x een polynomiaal verband tussen x en Y wordt verwacht, dan kan de regressie worden uitgevoerd met de verklarende variabele (1, x, x2 , . . ., xk ) in plaats van x. Wanneer men in het geval van een tweedimensionale variabele x = (x1 , x2 ) een simultaan kwadratisch verband verwacht tussen de variabelen (x1 , x2 ) en Y , dan kan de vector (1, x1 , x2 , x21 , x22 , x1 x2 ) in het model worden opgenomen, etc. Hier zien we dat de lineariteit van het lineaire regressiemodel betrekking heeft op de lineariteit in de regressieparameters, en niet zozeer in de onafhankelijke variabelen. In veel gevallen wordt er een intercept aan het model toegevoegd. In de praktijk komt dit er op neer dat de eerste verklarende variabele x1 gelijk aan 1 wordt gesteld, en de gemeten verklarende variabelen worden opgenomen in x2 , . . ., xp . Een lineair regressiemodel met intercept ziet er dus als volgt uit: Y = β1 + β2 x2 + . . . + βp xp + e. De regressieparameter β1 wordt het intercept genoemd en is de verwachting van Y in het geval dat de regressieparameters β2 , . . ., βp gelijk aan nul zijn. Het kan zinvol zijn de aannames van het regressiemodel te verzachten. In plaats van normaliteit van de meetfout zouden we alleen kunnen veronderstellen dat de meetfouten verwachting 0 bezitten, en ook de variantie van de meetfouten kan worden gemodelleerd als een functie van x. Deze modellen worden in dit boek niet besproken. 7.1.1

Enkelvoudige lineaire regressie

Stel dat de variabele Y afhangt van een eendimensionale reëelwaardige verklarende variabele x en dat we n waarnemingen (x1 , y1 ), . . ., (xn , yn ) hebben verkregen. Met behulp van een scatterplot van de waarnemingen (x1 , y1 ), . . ., (xn , yn ) kan inzicht verkregen worden in de relatie tussen x en Y . Indien deze relatie lineair lijkt te zijn, dan kunnen de waarnemingen gemodelleerd worden met een zogenaamd enkelvoudig lineair regressiemodel. Het enkelvoudige lineaire regressiemodel met intercept, wordt dan beschreven door (7.1)

Yi = α + βxi + ei ,

i = 1, . . ., n,

waarin de “meetfouten” e1 , . . ., en onderling onafhankelijke N (0, σ 2 )-verdeelde, niet-waarneembare stochastische grootheden zijn. Onder deze aanname zijn de variabelen Y1 , . . ., Yn eveneens onderling onafhankelijk en normaal verdeeld, waarbij de variabele Yi verwachtingswaarde α + βxi en variantie σ 2 heeft. De waarnemingen zijn dus niet identiek verdeeld. Wanneer er geen meetfouten zouden worden gemaakt, dan zou een exact lineair verband tussen Y en x bestaan. We nemen de parameterverzameling voor de parameter θ = (α, β, σ 2 ) zo groot mogelijk: (α, β) ∈ R2 en σ 2 > 0. 7.1.1.1

Schatten

Door de parameters α en β te schatten kunnen we het (lineaire) verband tussen Y en x achterhalen. Stelling 7.1. De maximum likelihood-schatters voor α, β en σ 2 in het enkelvoudige lineaire regressiemodel (7.1) zijn gelijk aan ˆ α ˆ = Y − βx,

sY βˆ = rx,Y , sx

n

σ ˆ2 =

1X ˆ i )2 , (Yi − α ˆ − βx n i=1

60


waar sx , sY en rx,Y de steekproefstandaardafwijkingen en de steekproefcorrelatie zijn (zie Notatie 3.4). Bewijs. De log likelihood-functie voor het model in (7.1) wordt gegeven door (α, β, σ 2 ) 7→ log

n Y

i=1

√

1

1

2πσ 2

e− 2 (Yi −α−βxi )

2

/σ 2

= − 21 n log 2π − 21 n log σ 2 −

n 1 X (Yi − α − βxi )2 . 2σ 2 i=1

Net als in Voorbeeld 4.8 gaat het maximaliseren van de log likelihood-functie in twee stappen. Maximalisatie van deze functie naar (α, β) is equivalent met minimalisatie van de kwadratische P vorm ni=1 (Yi − α − βxi )2 naar (α, β). Nulstellen van de partiële afgeleiden van de som naar α en β geeft het stelsel vergelijkingen n X i=1

(7.2)

n X i=1

ˆ i ) = 0, (Yi − α ˆ − βx

ˆ i )xi = 0. (Yi − α ˆ − βx

Met enig rekenwerk zijn de schatters voor α en β hieruit op te lossen: ˆ α ˆ = Y − βx, Pn Pn xi (Yi − Y ) (x − x)(Yi − Y ) sY i=1 Pn i βˆ = Pi=1 = = rx,Y , n 2 sx x) x) x (x − (x − i i i i=1 i=1

Dat deze oplossing een minimum geeft voor de kwadratische vorm voor iedere waarde van ˆ Substitueren σ 2 > 0, is na te gaan door bijvoorbeeld de Hessiaan-matrix te berekenen in (ˆ α, β). ˆ we de gevonden waarden van α ˆ en β in de log likelihood-functie en maximaliseren we deze functie vervolgens naar σ 2 dan vinden we n

σ ˆ2 =

1X ˆ i )2 (Yi − α ˆ − βx n i=1

als maximum likelihood-schatter voor σ 2 . De gevonden schatters voor α enP β zijn zogenaamde kleinste kwadratenschatters (KKn schatters), omdat ze de kwadraatsom i=1 (Yi − α − βxi )2 minimaliseren. Meetkundig komt dit neer op het minimaliseren van de som van de kwadratische verticale afstanden van de meetpunten (xi , Yi ) tot de beoogde regressielijn y = α + βx, zie Figuur 7.1; vandaar de naam. De KK-schatters α ˆ en βˆ zijn zuivere schatters voor het schatten van α en β (ook wanneer de aanname van normaal verdeelde meetfouten niet gemaakt wordt). P De KK-schatters α ˆ en βˆ zijn gevonden door de kwadraatsom ni=1 (Yi − α − βxi )2 te minimaliseren en voldoen derhalve aan de likelihood-vergelijkingen in (7.2). Meer algemeen kunnen schatters voor α en β gevonden worden door de vergelijkingen n X ˆ i )w(xi ) = 0, ψ(Yi − α ˆ − βx i=1

n X i=1

ˆ i )xi w(xi ) = 0, ψ(Yi − α ˆ − βx

ˆ voor geschikte functies ψ en w. Dit leidt in het algemeen tot andere op te lossen naar (ˆ α, β) schatters. De rol van de functie ψ en de gewichten w is vaak om de invloed van mogelijke ˆ i of de variabelen xi te verkleinen, dan wel om extreme waarden van de residuen Yi − α ˆ − βx de efficiëntie van de schatters te vergroten. Dit heet robuuste regressie.


61

• • •

15

•

•

•

•

•

•

10

• • • • • •

• •

•

5

•

0

•

0

2

4

6

8

10

Figuur 7.1. Een collectie meetpunten (xi , yi ) met de kleinste kwadratenlijn.

De maximum likelihood-schatter voor σ 2 kan worden geschreven als σ ˆ2 =

n n n X 1X 1 X ˆ i )2 = 1 (Yi − α ˆ − βx (Yi − Y )2 − βˆ2 (xi − x)2 . n i=1 n i=1 n i=1

Als a priori bekend zou zijn dat β = 0, dan hangt Yi niet van xi af, en zou de maximum likelihood-schatter van σ 2 gegeven worden door de eerste term aan de rechterkant. In het huidige model weten we niet of β gelijk is aan 0 en is de maximum likelihood-schatter voor σ 2 kleiner (tenzij βˆ = 0). Dit is intu¨ıtief redelijk: een deel van de variatie in Y is nu het gevolg van de variatie in x, en derhalve is de steekproefvariantie van de Yi een overschatting van σ 2 . Definitie 7.2. De uitdrukkingen SStot =

n X i=1

(Yi − Y )

2

en

SSres =

n X i=1

ˆ i )2 (Yi − α ˆ − βx

heten de totale kwadraatsom (Engels: total sum of squares) en de residuele kwadraatsom (Engels: residual sum of squares) of, meer volledig, de “residuele kwadraatsom na lineaire regressie op x”. De determinatiecoëfficiënt is gelijk aan 1−

SSres . SStot

P De totale kwadraatsom SStot is het minimum van ni=1 (Yi −α)2 over α, terwijl de residuele Pn kwadraatsom SSres het minimum is van i=1 (Yi − α − βxi )2 over (α, β) samen. Het tweede minimum isPnatuurlijk kleiner. Als SSres ongeveer even groot is als SStot , dan ligt SStot − SSres = βˆ2 ni=1 (xi −x)2 vlakbij 0 (dat wil zeggen, βˆ is ongeveer gelijk aan 0 bij genormaliseerde xi ) en heeft xi weinig voorspellende waarde voor Yi . De determinatiecoëfficiënt geeft de fractie door regressie op x verklaarde variantie (Engels: explained variance) weer en is te schrijven als Pn (xi − x)2 SSres 2 = βˆ2 Pni=1 1− = rx,Y . 2 SStot (Y − Y ) i i=1

Als de determinatiecoëfficiënt bijna gelijk aan 1 is, betekent dit dat de punten (xi , yi ) vrij behoorlijk op een rechte lijn liggen. Bij een determinatiecoëfficiënt van bijvoorbeeld 0.2, liggen de punten ofwel sterk verspreid rond een rechte lijn, of het lineaire regressiemodel is niet zinvol, omdat het verband tussen x en Y sterk niet-lineair is. De interpretatie van een determinatiecoëfficiënt is daarom niet eenvoudig. Merk op dat de schaal van de coëfficiënt kwadratisch is, hetgeen moeilijk te verdedigen is. De determinatiecoëfficiënt kan wel worden beschouwd als een zinvolle samenvatting van de data, en wordt standaard bij de rapportage van een regressieanalyse vermeld.

62


-1

0

y

1

2

3

ˆ Deze lijn kan worden gebruikt om de y-waarde De geschatte regressielijn is y = α ˆ + βx. bij een bepaalde x te voorspellen. Na invullen van de formules voor α ˆ en βˆ kunnen we de lijn herschrijven in de prettige vorm x−x y−Y = rx,Y . SY sx Aangezien |rx,Y | ≤ 1 vanwege de ongelijkheid van Cauchy-Schwarz, met strikte ongelijkheid tenzij de meetfouten exact 0 zijn, betekent dit dat, gemeten in standaarddeviaties, de voorspelde y dichter bij Y ligt dan de instelwaarde x bij x. Dit noemt men regressie naar het gemiddelde, in het bijzonder in het geval dat de standaarddeviaties van de twee variabelen ongeveer gelijk zijn. Is x de intelligentie van een vader en y de intelligentie van een zoon, dan leidde men hier ooit uit af dat de mensheid steeds middelmatiger wordt. Een verklaring voor “regressie naar het gemiddelde” is als volgt. We kunnen een x-waarde van een willekeurig geselecteerd individu uit een populatie (zoals de intelligentie van de vader) beschouwen als opgebouwd uit een toevallige en een systematische component. Als de x-waarde extreem groot is, dan is het een redelijke veronderstelling dat de toevallige component tot die relatieve grootte heeft bijgedragen. Bij het voorspellen van de afgeleide y-waarde is het verstandig hiermee rekening te houden: we voorspellen dat de y-waarde relatief minder extreem in de populatie van y-waarden zal liggen dan de x-waarde in de populatie van x-waarden. Deze interpretatie van “regressie naar het gemiddelde” binnen het kader van voorspellen wordt ondersteund door Figuur 7.2. Deze figuur laat twee regressielijnen zien. De gestippelde lijn lijkt de puntenwolk het best te volgen, maar de doorgetrokken lijn is de kleinste kwadratenlijn. De richtingscoëfficiënt van de gestippelde lijn is sy /sx , terwijl de kleinste kwadratenlijn de kleinere richtingscoëfficiënt rx,y sy /sx bezit. Dat de kleinste kwadratenlijn een betere voorspeller is, is te zien in het gebied tussen de twee verticale lijnen. De kleinste kwadratenlijn deelt de punten in deze strip (net als in iedere andere verticale strip) in ruwweg gelijke aantallen, terwijl de gestippelde lijn te hoog ligt.

-3

-2

-1

0

1

2

3

x

Figuur 7.2. De kleinste kwadratenlijn en “regressie naar het gemiddelde”.

7.1.1.2

Toetsen en betrouwbaarheidsintervallen

Omdat β in het algemeen de meest interessante parameter is in een enkelvoudig lineair regressiemodel worden in deze paragraaf toetsen en betrouwbaarheidsintervallen voor β afgeleid. Toetsen en betrouwbaarheidsintervallen voor de regressieparameter α kunnen op analoge wijze worden bepaald. In deze paragraaf beschrijven we de veel gebruikte t-toets en de likelihoodratiotoets. Om te bepalen of de onafhankelijke variabele x van lineaire invloed is op Y , toetsen we de nulhypothese H0 : β = 0. De gebruikelijke toetsingsgrootheid voor de meer algemene nulhypothese H0 : β = β0 is βˆ − β0 , T =q var c βˆ


63

waar βˆ de maximum likelihood-schatter voor β is en waar de variantie van βˆ geschat wordt met Pn 1 ˆ i )2 ˆ − βx i=1 (Yi − α n−2 P var c βˆ = . n 2 i=1 (xi − x)

Deze laatste schatter wordt verkregen door βˆ te schrijven als Pn Pn (xi − x)Yi (xi − x)(Yi − Y ) i=1 ˆ P = Pi=1 , β= n n 2 2 (x − x) i=1 i i=1 (xi − x) Pn ˆ i )2 /(n−2). (Merk te gebruiken dat var Yi = σ 2 en vervolgens σ 2 te schatten met i=1 (Yi −α ˆ −βx op dat deze schatter een factor n/(n − 2) afwijkt van de maximum likelihood-schatter voor σ 2 .) In Paragraaf 7.1.2.3 wordt in een algemener model bewezen dat de toetsingsgrootheid T onder H0 : β = β0 een tn−2 -verdeling heeft; het aantal vrijheidsgraden is gelijk aan het aantal waarnemingen minus het aantal geschatte regressiecoëfficienten. De nulhypothese wordt dan verworpen als |T | ≥ tn−2,1−α0 /2 waar α0 de onbetrouwbaarheidsdrempel van de toets is. Gebaseerd op deze t-toets kan een betrouwbaarheidsinterval voor β met onbetrouwbaarheid α0 worden bepaald: q ˆ c β. β = βˆ ± tn−2,1−α /2 var 0

7.1.2

Meervoudige Lineaire Regressie

In het meervoudige lineaire regressiemodel is de onafhankelijke variabele meerdimensionaal in plaats van eendimensionaal zoals in het enkelvoudige lineaire regressiemodel. Het meervoudige lineaire regressiemodel voor n afhankelijke variabelen Y1 , . . ., Yn met corresponderende pdimensionale verklarende variabelen (x1,1 , . . ., x1,p ), . . ., (xn,1 , . . ., xn,p ) wordt beschreven door Yi =

p X

βj xi,j + ei ,

i = 1, . . ., n,

j=1

waarin e1 , . . ., en onderling onafhankelijke normaal verdeelde stochastische grootheden zijn met verwachting 0 en eindige variantie σ 2 . De verklarende variabelen worden wederom niet stochastisch verondersteld, zodat de waarden xi,1 , . . ., xi,p als bekende constanten kunnen worden beschouwd. Het is handig om dit model in matrixnotatie weer te geven. De waarneming is een vector Y = (Y1 , . . ., Yn )T in Rn en de regressiecoëfficiënten vormen een vector β = (β1 , . . ., βp )T in Rp . Definiëren we de (n × p)-matrix X als de matrix met (i, j)-element xi,j , dan kunnen we het model schrijven als (7.3)

Y = Xβ + e,

waarin e = (e1 , . . ., en )T in Rn de foutenvector is. De matrix X wordt de design matrix genoemd. Merk op dat de notatie X hier gebruikt wordt voor een niet-stochastische matrix. In modellen met een intercept worden de elementen in de eerste kolom van de design matrix gelijk aan 1 verondersteld. De onbekende modelparameters zijn de regressiecoëfficiënten β = (β 1 , . . ., βp ) en de variantie σ 2 . 7.1.2.1

Keuze van de variabelen

In de praktijk zijn vaak veel mogelijke verklarende variabelen voorhanden om in de regressie te betrekken. De vraag is dan welke variabelen wel en welke niet opgenomen dienen te worden om een “best passend” model te krijgen. Zo zijn voor de prijs van een huis de ligging en inhoud informatief, net als het aantal vierkante meters en het bouwjaar. Met welke variabelen kan de verkoopprijs goed worden ingeschat? Of, in het geval van een ziekte, geven twee specifieke genen een goede beschrijving van de genetische component van een ziekte of is het verstandig alle ca. 20000 genen in de regressievergelijking te betrekken? Deze vragen zijn complexer dan ze lijken. Behalve dat de vragen op zeer uiteenlopende situaties betrekking hebben, lijkt ook het doel van de vragen uiteen te lopen. De vraag over de prijs van een huis is waarschijnlijk

64


afkomstig van de belastingdienst of een makelaar, welke een eenvoudige formule zoekt om op een enigszins objectieve manier een prijs vast te stellen. Best passend is dan het model dat de beste voorspellingen levert. De vraag uit de genetica is vermoedelijk causaal bedoeld: bepaalde genen zullen via biologisch-chemische processen invloed hebben op de vatbaarheid voor een bepaalde ziekte, andere genen niet. Als het alleen om het voorspellen van een ziekte gaat, is het niet zo erg ook enkele genen uit die tweede groep in de vergelijking op te nemen (ze kunnen samenhangen met genen die wel directe invloed hebben en toch voorspellende waarde hebben), maar voor het begrijpen van het ontstaan van de ziekte (causale verklaring) is een scherp onderscheid tussen de twee groepen genen essentieel. De gekozen verklarende variabelen kunnen zowel reëelwaardig als categorisch zijn. Een categorische verklarende variabele, ook wel nominale variabele genoemd, is een variabele waarvan de waarden een klassenindeling weergeven in plaats van een relevante numerieke grootte. Zo kunnen de waarden 0 en 1 voor man en vrouw gebruikt worden of kan een code staan voor een bepaalde regio. Om naast reëelwaardige variabelen ook categorische variabelen in een lineair regressiemodel op te nemen, is het gebruik van dummy-variabelen een standaard techniek. Een dummy-variabele is een indicatorvariabele en kan slechts de waarden 0 en 1 aannemen. Wanneer de categorische verklarende variabele k mogelijke klassen heeft, dan worden er k verklarende dummy-variabelen x1 , . . ., xk aan het regressiemodel (zonder intercept) toegevoegd. Wanneer de categorische variabele behoort tot klasse i, dan wordt de dummy-variabele x i gelijk gesteld aan 1, en de overige dummy-variabelen aan 0. Met deze k verklarende dummy-variabelen corresponderen in het lineaire regressiemodel k regressieparameters β1 , . . ., βk . Voor een waarneming Y behorende bij een verklarende variabele die tot de ie klasse behoort komt alleen de parameter βi in het regressiemodel, zie Tabel 7.1. x

x 1 x2 · · · x k

“1” “2” .. .

1 0 ··· 0 0 1 ··· 0

“k”

0 0 ··· 1

Pk

i=1

βi xi

β1 β2 .. . βk

Tabel 7.1. Definitie van dummy-variabelen x1 , . . ., xk voor regressie op een categorische variabele x met k klassen met labels “1”, . . ., “k”.

Op deze manier is de parameter βi in feite het intercept voor de klasse i. Er wordt dan ook geen extra intercept meer toegevoegd aan het model. Wanneer men toch een intercept wil opnemen in het model, dient het aantal dummy-variabelen een minder te zijn dan het aantal klassen. Dan wordt bijvoorbeeld de dummy-variabele voor de eerste klasse eruit gelaten. De parameter β1 is dan het gebruikelijke intercept en de parameter βi (i = 2, . . ., k) geeft het effect van klasse i op de afhankelijke variabele Y weer ten opzichte van klasse 1. Ook in het geval van twee categorische verklarende variabelen (bijvoorbeeld wanneer zowel regio als geslacht in het model wordt opgenomen) wordt voor de tweede variabele een dummy-variabele minder dan het aantal bijbehorende klassen opgenomen in het model. Dit is nodig om de design matrix volle rang te laten behouden, zie Paragraaf 7.1.2.2. Hoewel het gebruikelijk is om voor een dummy-variabele alleen de waarden 0 en 1 toe te laten, ligt het in sommige situaties voor de hand aangepaste waarden te hanteren. Zo is in Voorbeeld 2.5 een dummy-variabele in het model opgenomen die de waarden -1 en 1 aanneemt. De keuze voor de waarden van de dummy-variabele evenals het wel of niet opnemen van een intercept hangt af van de gewenste interpretatie van de parameters βi . Wanneer er alleen categorische variabelen zijn, is het model voor variantie-analyse van toepassing. 7.1.2.2

Schatten

De volgende stelling geeft de maximum likelihood-schatters voor de parameters in het meervoudige lineaire regressiemodel.


65

Stelling 7.3. Wanneer de design matrix X in het meervoudige lineaire regressiemodel (7.3) volle rang heeft, worden de maximum likelihood-schatters voor β en σ 2 gegeven door βˆ = (X T X)−1 X T Y,

σ ˆ2 =

ˆ 2 kY − X βk . n

Bewijs. De log likelihood-functie voor het meervoudige lineaire regressiemodel wordt gegeven door Pp n 1 Y 1 βj xi,j )2 /σ 2 − (Y − j=1 √ e 2 i (β, σ 2 ) 7→ log 2 i=1 2πσ 1 = − 21 n log 2π − 21 n log σ 2 − 2 kY − Xβk2 , 2σ waar k·k de notatie voor de Euclidische norm is. Geheel analoog aan het geval van enkelvoudige lineaire regressie wordt eerst de schatter voor β afgeleid voor willekeurige σ 2 en daarna de schatter voor σ 2 . De maximum likelihood-schatter voor β is de KK-schatter βˆ die de functie β 7→ kY − Xβk2 ,

β ∈ Rp

minimaliseert. Voor iedere β behoort de vector Xβ tot het bereik (de kolommenruimte) van de matrix X, gezien als afbeelding X: Rp → Rn . De KK-schatter βˆ die kY − Xβk2 minimaliseert, is derhalve de vector zodanig dat X βˆ het element van het bereik van X is dat zo dicht mogelijk bij de vector Y ligt. In de lineaire algebra heet X βˆ de projectie van Y op het bereik van X. De projectie ten opzichte van de Euclidische norm voldoet aan de orthogonaliteitsrelatie ˆ Xγi = γ T X T (Y − X β) ˆ = 0, hY − X β,

∀γ ∈ Rp .

Met andere woorden, het “residu” Y − X βˆ staat loodrecht op ieder willekeurig element in de kolommenruimte van X, dat in algemene vorm geschreven kan worden als Xγ voor een γ ∈ R p . ˆ = De eis dat dit gelijk aan nul moet zijn voor willekeurige γ ∈ Rp betekent dat X T (Y − X β) 0. Dit is de zogenaamde normaalvergelijking. Uit de aanname dat X volle rang heeft, volgt dat X T X inverteerbaar is, en daarmee volgt dat βˆ = (X T X)−1 X T Y . Dan is X βˆ gelijk aan X(X T X)−1 X T Y , hetgeen inderdaad de projectie van Y op de kolommenruimte van X is, omdat X(X T X)−1 X T de projectiematrix is die projecteert op deze ruimte. In de tweede stap wordt βˆ voor β in de log likelihood-functie gesubstitueerd. Het is eenˆ 2 /n als maximum likelihood-schatter voor σ 2 oplevvoudig na te gaan dat dit σ ˆ 2 = kY − X βk ert. De maximum likelihood-schatter βˆ is zuiver: Eβˆ = (X T X)−1 X T EY = (X T X)−1 X T Xβ = β. Definitie 7.4. De residuen van de regressie van Y op X zijn de co¨ ordinaten van de vector ˆ De uitdrukkingen Y − X β. SStot = kY − Y 1k2

en

ˆ 2, SSres = kY − X βk

met Y 1 = (Y , . . ., Y ), heten de totale kwadraatsom en de residuele kwadraatsom. De determinatiecoëfficiënt is gelijk aan SSres 1− . SStot De determinatiecoëfficiënt neemt, net als in het geval van enkelvoudige lineaire regressie, waarden aan tussen 0 en 1. Dit kan als volgt worden ingezien. De vector Y 1 = (Y , . . ., Y ) is de beste voorspelling van Y in een model bestaande uit alleen een geschat intercept. Het is de projectie van Y op de eendimensionale lineaire ruimte opgespannen door de vector 1: = (1, 1, . . ., 1). Omdat we voor een model met intercept hebben gekozen, is deze ruimte bevat in de kolommenruimte van X. Daarom staat de residuvector Y − X βˆ loodrecht op de vector 1,

66


ˆ 1i = 0. Daarmee volgt ook dat hY − X β, ˆ Y 1i = 0. Verder hadden we voor iedere γ in hY − X β, p ˆ ˆ X βi ˆ = 0. Hieruit volgt dat Y − X βˆ R dat hY − X β, Xγi = 0, dus in het bijzonder hY − X β, ˆ loodrecht staat op zowel de vector Y 1 als de vector X β. Daarmee staan de vectoren Y − X βˆ en X βˆ − Y 1 loodrecht op elkaar. Met de regel van Pythagoras volgt nu dat het kwadraat van ˆ + (X βˆ − Y 1), gelijk is aan de lengte van de som van deze twee vectoren, Y − Y 1 = (Y − X β) ˆ 2 + kX βˆ − Y 1k2 . kY − Y 1k2 = kY − X βk Het linkerlid in deze vergelijking is gelijk aan SStot en de eerste term in het rechterlid is gelijk aan SSres . We zien dat 0 ≤ SSres ≤ SStot en dat de determinatiecoëfficiënt tussen 0 en 1 2 ligt. Men kan laten zien dat, analoog aan 1 − SSres /SStot = rx,Y in het enkelvoudige lineaire 2 regressiemodel, voor het meervoudige model geldt dat 1 − SSres /SStot = rX ˆ . β,Y Wanneer men geen intercept in het model opneemt, moet ook in de berekening van de totale kwadraatsom niet uitgegaan worden van een intercept. In dat geval hanteert men SS tot = kY k2 . 7.1.2.3

Toetsen

Net als in het enkelvoudige lineaire regressiemodel, zijn er twee belangrijke soorten toetsen om de invloed van één of meerdere onafhankelijke variabelen op de afhankelijke variabele Y te bepalen in het geval van een meervoudig lineair regressiemodel. Een veel voorkomende nulhypothese bij het meervoudige lineaire regressiemodel is H 0 : βj = 0 voor een zekere j ∈ {1, . . ., p}. Als βj = 0 heeft de j e verklarende variabele geen invloed op de afhankelijke variabele Y . Onder de nulhypothese kan het regressiemodel daarom worden vereenvoudigd door deze verklarende variabele uit het model weg te laten. Dit betekent concreet dat de j e kolom uit de design matrix verwijderd wordt, evenals de j e co¨ ordinaat van β. De nieuwe design matrix noteren we met X−j en de verkorte vector met regressieparameters met β−j . De maximum likelihood-schatter voor β−j wordt op analoge wijze als de schatter voor β T T afgeleid en is gelijk aan βˆ−j = (X−j X−j )−1 X−j Y . de toetsingsgrootheid T =

βˆj , σ ˆ 2 (X T X)−1 j,j

waarvoor men kan laten zien dat T onder de nulhypothese een t-verdeling met n − p vrijheidsgraden volgt. Voorbeeld 7.5 (Lichaamslengte). In Voorbeeld 2.5 wordt een meervoudig lineair regressiemodel beschreven voor het schatten van de eindlengte van een kind op basis van de lichaamslengte van de (biologische) ouders en het geslacht van het kind. De regressieparameters voor de verklarende variabelen “lichaamslengte van de vader” en “lichaamslengte van de moeder” zijn in dat voorbeeld gelijk aan 1/2 genomen, zodat het geschatte model gemakkelijk te interpreteren is. In deze paragraaf zullen we deze aanname achterwege laten en een meervoudig lineair regressiemodel schatten op basis van onze eigen data. Voor Y de eindlengte van een kind, x2 de lengte van de vader, x3 de lengte van de moeder en x4 een variabele voor het geslacht van het kind, ziet het meervoudige lineaire regressiemodel er als volgt uit Y = β1 + β2 x2 + β3 x3 + β4 x4 + e, met e een normaal verdeelde stochastische variabele met verwachting 0 en variantie σ 2 . De eerste verklarende variabele, x1 , hebben we gelijk aan 1 genomen, zodat het model een intercept heeft. De verklarende variabele x4 geeft aan of het kind een jongen of een meisje is. Omdat er een intercept in het model is, hebben we aan een dummy-variabele genoeg. We willen dat β4 gelijk is aan de helft van het gemiddelde verschil in lengte tussen mannen en vrouwen en kiezen daarom x4 voor een meisje gelijk aan -1 en voor een jongen aan 1. Aangezien gemiddeld genomen jongens langer zullen worden dan meisjes, zal β4 positief zijn. Onze data bestaan uit de eindlengtes van 111 jong-volwassenen (44 mannen en 67 vrouwen), zijn of haar geslacht en de lengte van hun ouders. Voor elk van de regressieparameters toetsen we of de waarde significant van nul afwijkt met de t-toets zoals beschreven staat


67

in de vorige paragraaf. De onbetrouwbaarheid van de toetsen wordt gelijk aan 0.05 genomen. Alle toetsen worden verworpen en het uiteindelijke model met geschatte regressieparameters wordt gegeven door: Y = 2.52 + 0.46x2 + 0.55x3 + 6.27x4 + e, waarbij e normaal verdeeld verondersteld wordt met verwachting nul en (geschatte) variantie 25.78. De determinatiecoëfficient van het model is 0.69. Om te onderzoeken of de normaliteitsaanname plausibel is, kan men scatterplots tekenen en eventueel aanvullende toetsen uitvoeren. Bovendien is het verstanding om te onderzoeken of het redelijk is om een lineair verband te veronderstellen. Ondanks dat de geschatte regressieparameters in bovenstaand model niet overeenkomen met de schattingen in de Vierde Landelijke Groeistudie (zie (2.1)) liggen de verwachte eindlengtes niet ver van elkaar. Zo zijn de streeflengtes van kinderen van een man van 180 cm en een vrouw van 172 cm in bovenstaand model gelijk aan 186.2 cm (zoon) en 173.7 cm (dochter), terwijl deze lengtes volgens het regressiemodel in (2.1) 187 cm en 174 cm zijn. In bovenstaand model is de invloed van de lengte van de moeder op de eindlengte van het kind groter dan die van de lengte van de vader. Aangezien het meer voor de hand ligt dat de invloeden gelijk zijn, kunnen we eveneens het volgende regressiemodel schatten: Y = β1 + β2 (x2 + x3 ) + β4 x4 + e. De geschatte regressieparameters voor dit model zijn βˆ1 = 3.47, βˆ2 = 0.50 en βˆ4 = 6.30. Ook voor dit tweede geschatte model is de determinatiecoëfficiënt gelijk aan 0.69. Dit model komt meer overeen met het model in (2.1). De geschatte waarden voor β1 en β4 zijn evenwel iets lager in ons geschatte model waardoor de voorspelde eindlengtes op basis van ons model wat lager zijn dan op basis van het model in (2.1). Een mogelijke oorzaak voor het verschil tussen de geschatte modellen is dat wij slechts gegevens hebben van 111 jongvolwassenen, terwijl het model dat in de Vierde Landelijke Groeistudie is geschat gebaseerd is op veel meer gegevens. De geschatte modellen in dit voorbeeld zijn dus minder betrouwbaar.

7.1.3

Het ontwerp van een experiment

Bij veel onderzoeken en experimenten kunnen de waarden van de onafhankelijke variabelen gekozen worden. Denk maar eens aan een natuurkundig experiment waarbij de remweg wordt gemeten bij verschillende snelheden en massa’s van een auto, de instelwaarden van een machine bij een produktie proces, of een onderzoek waarbij men is ge¨ınteresseerd in de voorspellende waarde van lichaamslengte op het lichaamsgewicht. Immers, bij het onderzoek naar de remweg kan voor aanvang worden bepaald bij welke snelheden en massa’s de remweg zal worden gemeten, de knoppen van een machine kunnen handmatig worden ingesteld, en bij het onderzoek naar lichaamsgewicht kunnen deelnemers worden geselecteerd op hun lichaamslengte. Een geschikte keuze van de waarden van de onafhankelijke variabelen, indien mogelijk, is niet alleen van belang bij het onderzoek naar de aard van de relatie tussen de onafhankelijke en afhankelijke variabele (bijvoorbeeld lineair of kwadratisch), maar ook om de regressieparameters zo nauwkeurig mogelijk te schatten (kleine variantie). De vraag is nu hoe de waarden van de onafhankelijke variabelen het beste gekozen kunnen worden. Alle verklarende variabelen gelijk kiezen is een mogelijkheid, maar zal het onmogelijk maken de relatie tussen de onafhankelijke en afhankelijke variabele te onderzoeken en de bijbehorende regressieparameter te schatten. De waarden gelijk verdelen over het domein van de variabele, is misschien een goede manier om een eventueel verband te onderzoeken, maar blijkt in het (enkelvoudige) lineaire regressiemodel niet optimaal te zijn als we de variantie van de schatters van de regressieparameters zo klein mogelijk willen krijgen. In deze paragraaf zullen we onderzoeken hoe de waarden van de onafhankelijke variabelen gekozen moeten worden om de regressieparameters zo nauwkeurig mogelijk te schatten onder de aanname dat een enkelvoudig lineair model bij de data past. Veronderstel dat Y1 , . . ., Yn onderling onafhankelijk zijn en dat Yi = α + βxi + ei

68


met x1 , . . ., xn de waargenomen verklarende variabelen en e1 , . . ., en onderling onafhankelijke fouten met verwachting 0 en variantie σ 2 . De kleinste kwadratenschatters voor de regressieparameters α en β zijn zuiver en hun varianties zijn gelijk aan Pn σ2 x2 σ2 var βˆ = Pn . var α ˆ = Pn i=1 i 2 2 n i=1 (xi − x) i=1 (xi − x)

Voor zowel α als β hangt de variantie af van de waarden van de onafhankelijke variabelen x1 , . . ., xn . Dat betekent dat de nauwkeurigheid waarmee de regressieparameters geschat kunnen worden afhankelijk is van de keuze van deze variabelen. Met name de nauwkeurigheid van βˆ is interessant, omdat deze parameter het verband tussen de afhankelijke en onafhankelijke Pn variabele weergeeft. Om de variantie van βˆ te minimaliseren, moet i=1 (xi − x)2 zo groot mogelijk worden; ofwel de variatie in de onafhanelijke variabele moet zo groot mogelijk genomen worden. Dit kan bereikt worden door de helft van de variabelen zo klein en de andere helft zo groot mogelijk te nemen. Deze keuze van de waarden minimaliseert de variantie, maar is niet geschikt wanneer modelaannamen als lineairiteit in de onafhankelijke variabele gevalideerd moeten worden. Echter, soms is uit de natuurkunde bekend wat de relatie tussen de afhankelijke en onafhankelijke variabelen is, omdat dat wordt beschreven door een natuur- of scheikundige wet en wil men slechtsP de waarde van een natuur- of scheikundige constante verifiëren. In dat n geval volstaat het om i=1 (xi − x)2 te maximaliseren. Als bekend is dat het verband tussen de afhankelijke en onafhankelijke variabele wordt beschreven door een enkelvoudig lineair regressiemodel zonder intercept, dan verandert de kleinste kwadratenschatter voor β, net als zijn variantie en de optimale waarden voor de onafhankelijke variabele. Voorbeeld 7.6 (Relatie lichaamslengte en -gewicht). Veronderstel dat we op zoek zijn naar een relatie tussen lichaamslengte en gewicht, zodat op basis van de lichaamslengte een voorspelling van het lichaamsgewicht gemaakt kan worden. Het ligt voor de hand dat de variabelen een positieve relatie hebben, maar welke relatie is niet bekend. Het is daarom niet verstandig om alleen korte en lange mensen te vragen naar hun lengte en gewicht, maar zeker ook mensen met een meer gemiddelde lichaamslengte.

In de natuur- en scheikunde zijn relaties tussen sommige variabelen vast. Zo is uit te rekenen hoelang het duurt totdat een bol op de grond valt als je het op een hoogte van 2 meter loslaat, of hoe groot de slingertijd van een bol is als de lengte van het koord bekend is. Als we meerdere keren hetzelfde experiment uitvoeren en de gevonden data uitzetten in een grafiek is een duidelijk verband te zien. De waarnemingen zullen echter niet precies overeenkomen met de theoretische waarde en dus zullen de waarnemingen in de scatterplot lichtjes afwijken van de lijn of kromme die bekend is uit de theorie. De kleine afwijkingen worden niet veroorzaakt doordat de relatie tussen afhankelijke en de onafhankelijke variabele niet exact klopt, maar doordat er kleine meetfouten gemaakt worden bij het meten en aflezen van zowel de waarden van de afhankelijke als onafhankelijke variabele. Investeert men in betere meetapparatuur, dan kun je verwachten dat de meetfouten kleiner worden en de waarnemingen dichter bij de theoretische lijn of kromme komen te liggen. Lengte en gewicht hebben duidelijk een positief verband, maar dit verband is niet exact; de waarnemingen zullen mooi verspreid om een denkbeeldige lijn of kromme liggen. Het nauwkeuriger meten van lengte en gewicht, zal misschien iets uitmaken in de scatterplot, maar de mate van correlatie tussen de twee variabelen zal nauwelijks veranderen.

7.2

Niet-lineaire en niet-parametrische regressie

In het lineaire regressiemodel is de voorwaardelijke verwachting van de afhankelijke variabele Y gegeven de verklarende variabele X = x, x 7→ E(Y | X = x), een lineaire functie van de parameter β. Deze regressiefunctie kan worden vervangen door een algemenere functie, E(Y | X = x) = f (x).

7: Opgaven

69

Bovenstaande vergelijking wordt ook wel de regressievergelijking genoemd. Als f = fθ bekend is op de parameter θ na en fθ een niet-lineaire functie is van θ, dan spreken we van niet-lineaire regressie. Net als bij lineaire regressie kan de parameter θ geschat worden met de kleinste kwadratenmethode. De KK-schatter voor θ minimaliseert het criterium θ 7→

n X i=1

2 Yi − fθ (xi ) .

Veelal is een numeriek algoritme noodzakelijk voor de bepaling van de kleinste kwadratenschatting. Wanneer de meetfouten in het model Y = fθ (x) + e normaal verdeeld zijn, is de KK-schatter voor θ tevens de maximum likelihood-schatter. Een voorbeeld van een geparametriseerd niet-lineair regressiemodel is de tijdscurve E(Y | X = x) = gθ (x) = θ0 + θ1 x + θ2 e−θ3 x . Bij waarnemingen y1 , . . ., yn op tijdstippen x1 , . . ., xn wordt de KK-schatter voor θ = (θ0 , θ1 , θ2 , θ3 ) gevonden door n X i=1

yi − θ0 − θ1 xi − θ2 e−θ3 xi

2

te minimaliseren naar θ.

Opgaven 1. Men wil onderzoeken of er verband is tussen (gemiddelde) eindexamencijfers en studieduur (in maanden). De volgende gegevens, betrekking hebbend op een steekproef van in 1965 in Nijmegen aangekomen wiskundestudenten, zijn beschikbaar. cijfer studieduur

8 82

7 80

7 66

7 77

7.5 79

7 75

6.5 58

9 46

7 58

9 56

8 70

7.5 55

(i) Bepaal schattingen voor de regressiecoëfficiënten in het lineaire regressiemodel m.b.v. de methode der kleinste kwadraten. (ii) Bereken de fractie verklaarde variantie. (iii) Maak een plaatje met daarin de waarnemingen en de aangepaste rechte lijn. 2. Beschouw het enkelvoudige lineaire regressiemodel, maar veronderstel dat a priori bekend is dat α = 0. (i) Bepaal de kleinste kwadratenschatter voor β. (ii) Ga na of deze schatter zuiver is. (iii) Bepaal de variantie van deze schatter. (iv) Veronderstel dat de waarden x1 , . . ., xn door de onderzoeker kunnen worden ingesteld binnen een gegeven interval [0, b]. Wat is dan een optimale instelling om de beste schatting van β te krijgen? (v) Zie je een reden om in de praktijk deze optimale instelling toch niet te kiezen? 3. Beschouw het enkelvoudige regressiemodel Yi = α + βxi + ei , waarin α en β onbekende parameters zijn, x1 , . . ., xn bekende constanten en e1 , . . ., en onderling onafhankelijk normaal verdeeld zijn met Eei = 0 en Ee2i = zi σ 2 voor bekende positieve getallen z1 , . . ., zn . Een dergelijk model met fouten van verschillende nauwkeurigheid noemt men heteroskedastisch. Bepaal de maximum likelihoodschatters van α en β.

70


4. Veronderstel dat we de valversnelling experimenteel willen bepalen met behulp van een slingerproef. Daartoe bevestigen we een massieve bol aan een dun koord. Vervolgens slingeren we de bol en meten we de slingertijd. De slingertijd is de tijdsduur tussen twee momenten waarop de bol zich op hetzelfde uiteinde bevindt. Als de uitwijking van de bol niet te groot is, is de slingertijd constant en dus onafhankelijk √ √ van de grootte van de uitwijking van de bol. Volgens de Wet van Newton geldt: T = 2π l/ g, waarbij T (in s) staat voor de slingertijd, g (in m/s2 ) voor de valversnelling en l (in m) voor de lengte van het koord. We gaan ervan uit dat in onze opzet aan deze wet is voldaan. We voeren de slingerproef n maal uit met gekozen koordlengten l1 , . . ., ln . Dit geeft de waarnemingen (l1 , T1 ), . . ., (ln , Tn ). Bij het uitvoeren van het experiment worden kleine meetfouten gemaakt, waardoor de metingen niet exact de Wet van Newton zullen volgen, maar slechts bij benadering. (i) Beschrijf een geschikt enkelvoudig lineair regressiemodel. (ii) Bepaal de kleinste kwadratenschatter voor de valversnelling g. √ (iii) Bepaal de variantie van de kleinste kwadratenschatter van de regressieparameter 2π/ g die voor de helling staat. (iv) Hoe zou je de koordlengtes in het experiment kiezen, zodat de variantie in het vorige onderdeel minimaal is. (v) Beredeneer (of bereken) hoe je de koordlengtes in het experiment zou kiezen, zodat de variantie van de kleinste kwadratenschatter van de valversnelling g minimaal is.

8 Appendix A: Elementen uit de Kansrekening

Deze appendix bevat een aantal onderdelen uit de kansrekening die van belang zijn bij het lezen van dit boek. Het doel is om deze stof kort weer te geven. Voor nadere toelichting, inclusief bewijzen van stellingen, voorbeelden en toepassingen, wordt verwezen naar tekstboeken over waarschijnlijkheidsrekening, zoals S. Ross, A First Course in Probability, Prentice-Hall en R. Meester, A Natural Introduction to Probability Theory, Birkh¨ auser.

8.1

Verdelingen

De basis van alle statistische procedures is een waarneming X waarbij onzekerheid, toeval of een andere vorm van willekeur een rol speelt. Net als in de kansrekening wordt in de statistiek de onzekerheid wiskundig vertaald door aan X een kansverdeling toe te kennen, Definitie 8.1. Een stochastische variabele (of stochastische grootheid, of stochast) is een waarneming onderhevig aan onzekerheid, beschreven door een kansverdeling. De verzameling van alle mogelijke uitkomsten van X wordt de uitkomstenruimte genoemd en wordt genoteerd als Ω. Een kansverdeling is een voorschrift dat aangeeft wat de kans is om de waarneming X in A te vinden voor (bijna) iedere deelverzameling A ⊆ Ω, P(X ∈ A). Kansverdelingen hebben drie eigenschappen (i) P(X ∈ Ω) = 1, (ii) 0 ≤ P(X ∈ A) ≤ 1 voor alle A ⊆ Ω, (iii) (σ-additiviteit) voor A1 , A2 , . . ., ⊆ Ω disjunct (Ai ∩ Aj = ∅ als i 6= j), geldt P(X ∈

∞ [

i=1

Ai ) =

∞ X i=1

P(X ∈ Ai ).

Uit deze drie eigenschappen kunnen alle andere algemeen geldende eigenschappen van een kansverdeling direct worden afgeleid. Een voorbeeld is de eigenschap P(X ∈ A) ≤ P(X ∈ B) voor A, B ⊆ Ω met A ⊆ B. Voor andere eigenschappen van kansverdelingen wordt verwezen naar de boeken over kansrekening, zoals eerder genoemd. In sommige gevallen zijn we niet ge¨ınteresseerd in de stochast X zelf, maar in een functie g van deze stochast, bijvoorbeeld g(X) = X 2 . In veel gevallen is de functie g gedefinieerd op de hele reële rechte en is g(X) daarom goed gedefinieerd voor alle reëelwaardige stochastische grootheden X † . Stel dat de kansverdeling van X bekend is, dan ligt de verdeling van Y = g(X) † De voorwaarde dat X met kans 1 ligt in het domein van g zorgt er in het algemene geval voor dat g(X) goed gedefinieerd is. Strikt genomen bestaat er een voorwaarde voor g (meetbaarheid), maar in dit boek komen niet-meetbare functies niet aan de orde.

72

8: Appendix A: Elementen uit de Kansrekening

vast door

P(Y ∈ A) = P(g(X) ∈ A) = P(X ∈ g −1 (A))

waarbij g −1 (A) het volledig origineel van A onder g genoemd wordt: g −1 (A) = {x: g(x) ∈ A}.

(De notatie g −1 zou kunnen suggereren dat de inverse van g vereist is voor de definitie van het volledig origineel. Dat is echter niet zo; in het geval dat g niet inverteerbaar is, is de verzameling in het rechterlid toch goed gedefinieerd.) 8.1.1

Discrete en continue verdelingen

Er bestaan twee basis soorten kansverdelingen: discrete verdelingen en continue verdelingen. Een discrete kansverdeling wordt gekenmerkt door een eindige of aftelbare verzameling van mogelijke uitkomsten van de stochastische variabele, terwijl een stochast met een continue verdeling uitkomsten kan hebben in een interval van de reële rechte. Met zowel iedere discrete als continue verdeling correspondeert een kansdichtheid (of dichtheid) en een verdelingsfunctie. In het geval van een discrete verdeling kent de kansdichtheid aan iedere mogelijke uitkomst een gewicht (kansmassa) toe, gelijk aan de kans op die uitkomst. Deze gewichten zijn nietnegatief en sommeren tot 1. De kans op een uitkomst in A, een deelverzameling van de uitkomstenruimte Ω, is gelijk aan X p(ω) P(X ∈ A) = ω∈A

waar p(ω) = P(X = ω). Voorbeelden van veelgebruikte discrete verdelingen zijn de Bernoulli-, de binomiale, de Poisson-, de geometrische, de hypergeometrische en de negatief binomiale verdeling. Als X continu verdeeld is over (een deel van) de reële rechte, hanteren we een kansdichtheidsfunctie f : R → R, die we ook kortweg kansdichtheid noemen. De sommatie bij discrete verdelingen wordt vervangen door integratie bij continue verdelingen. De kans op een uitkomst in A ⊆ R van de continu verdeelde stochastische variabele X met kansdichtheid f wordt gegeven door Z P(X ∈ A) = f (x) dx. A

Voorbeelden van bekende continue verdelingen zijn de uniforme, de normale, de exponentiële, de Cauchy-, de chikwadraat-, de t-, de Gamma- en de Bèta-verdeling.

8.1.2

Verdelingsfuncties

Kansdichtheden vormen een manier om een verdeling te specificeren. Een andere, equivalente manier om dit te doen, is door middel van een zogenaamde verdelingsfunctie. Definitie 8.2. Zij X een stochastische variabele, verdeeld volgens een bepaalde kansverdeling. De verdelingsfunctie F behorende bij die kansverdeling wordt gedefinieerd door F (x) = P(X ≤ x). De verdelingsfunctie is een monotoon stijgende functie, dat wil zeggen, als x ≤ y dan is F (x) ≤ F (y). De definitie van de verdelingsfunctie is in bovenstaande vorm geldig voor zowel discreet als continu verdeelde stochastische grootheden die reëelwaardig zijn. Voor een reëelwaardige, discreet verdeelde stochastische grootheid X kan de relatie tussen de kansdichtheid p en de verdelingsfunctie F als volgt worden uitgedrukt: X p(ω). F (x) = P(X ≤ x) = ω≤x

De verdelingsfunctie vertoont sprongen in alle punten die tot de mogelijke uitkomsten van X behoren. Tussen deze sprongen in is de verdelingsfunctie constant. De grootte van de sprong in

8.2: Verwachting en variantie

73

het punt ω is gelijk aan p(ω). Discrete verdelingen kunnen zodoende op twee manieren worden gespecificeerd: ofwel met de kansdichtheid p (de verdelingsfunctie F wordt gevonden door sommatie), ofwel met de verdelingsfunctie F (de kansdichtheid p volgt uit de spronggroottes). Voor een continu verdeelde stochastische grootheid X kan de relatie tussen de kansdichtheid f en de verdelingsfunctie F als volgt worden uitgedrukt: Z x F (x) = P(X ≤ x) = f (u) du. −∞

De verdelingsfunctie F kan daarom worden opgevat als de primitieve van de kansdichtheid f . Omgekeerd is f de afgeleide van F , f (x) = F 0 (x). Hieruit blijkt dat continue verdelingen eveneens kunnen worden vastgelegd op twee manieren: ofwel door de kansdichtheid f (de verdelingsfunctie F volgt uit integratie), ofwel door de verdelingsfunctie F (de kansdichtheid volgt uit differentiatie). Met de verdelingsfunctie is het eenvoudig om kansen voor intervallen van de vorm (c, d] uit te rekenen, P(c < X ≤ d) = P(X ≤ d) − P(X ≤ c) = F (d) − F (c).

Voor discrete verdelingen is het van belang of het interval open, gesloten of halfopen gekozen wordt. De kans P(c ≤ X ≤ d) is bijvoorbeeld groter dan P(c < X ≤ d) wanneer P(X = c) > 0, omdat P(c ≤ X ≤ d) = P(c < X ≤ d)+P(X = c). Aangezien voor continu verdeelde stochasten geldt dat P(X = c) = 0 voor alle c, speelt de keuze van open of gesloten intervallen daar geen rol.

8.2

Verwachting en variantie

De verwachting en variantie van een verdeling zijn eigenschappen die respectievelijk de locatie en de spreiding van de verdeling weergeven. De locatie is een punt waaromheen de verdeling zich centreert, terwijl de spreiding een maat is voor de breedte van de verdeling rondom zijn locatie. Er zijn meerdere eigenschappen die als locatie- of spreidingsbegrip kunnen dienen. Verwachting en variantie zijn voorbeelden die veel gebruikt worden. De verwachtingswaarde (of verwachting) E(X) van een stochastische grootheid X correspondeert met het begrip van gewogen gemiddelde. De wegingsfactoren zijn gebaseerd op de kansdichtheid. Wanneer er geen verwarring kan optreden schrijven we kortweg EX voor de verwachting. Voor een discreet verdeelde stochastische grootheid X met kansdichtheid p wordt EX gegeven door X EX = ω p(ω), ω∈Ω

Voor een continu verdeelde stochastische grootheid X met kansdichtheid f is de verwachtingswaarde gedefinieerd door Z ∞ EX = x f (x) dx. −∞

De verwachting van g(X), voor g een functie, is Z ∞ E(g(X)) = g(x) f (x) dx. −∞

De verwachtingswaarde van X of g(X) is niet altijd goed gedefinieerd. Het is mogelijk dat de integraal of som niet convergeert. De verwachting van een Cauchy-verdeelde stochastische grootheid bestaat bijvoorbeeld niet. In dit boek zullen we veronderstellen dat alle integralen die we gebruiken ook daadwerkelijk bestaan. De variantie is de verwachting van het kwadraat van de afstand van X tot zijn verwachtingswaarde, var(X) = E(X − EX)2 .

74


Het is eenvoudig na te gaan dat de variantie te schrijven is als var(X) = E(X 2 ) − (EX)2 . Deze schrijfwijze is in veel gevallen handig in de berekening van de variantie van een stochastische grootheid. De verwachting E(X 2 ) wordt gevonden uit E(g(X)) met g(X) = X 2 . De variantie is groot indien X met relatief grote kans op aanzienlijke afstand van EX aangetroffen wordt. Dit karakteriseert het spreidingsbegrip. De covariantie van twee stochastische grootheden X en Y is gelijk aan cov(X, Y ) = E (X − EX)(Y − EY ) = E(XY ) − EXEY.

Uit de definities van verwachtingswaarde en variantie kunnen de volgende rekenregels worden afgeleid E(a + bX) = a + b EX var(a + bX) = b2 var(X) E(X + Y ) = EX + EY var(X + Y ) = var(X) + var Y + 2 cov(X, Y ).

8.3

Standaard verdelingen

In deze paragraaf worden voorbeelden gegeven van discrete en continue verdelingen die vaak voorkomen. 8.3.1

Discrete verdelingen

Voorbeeld 8.3 (Bernoulli-verdeling). De stochastische grootheid X heeft de Bernoulliverdeling (of alternatieve verdeling) met parameter p ∈ [0, 1] als P(X = 0) = 1 − p

en

P(X = 1) = p.

Deze kansmassafunctie is ook te schrijven als P(X = x) = px (1 − p)1−x

x ∈ {0, 1}.

De verwachtingswaarde en variantie zijn in dat geval gelijk aan EX = p en Var(X) = p(1 − p). Als X1 , . . ., Xn onderling onafhankelijke Bernoulli-verdeelde stochasten zijn met parameter p, dan is X1 + . . . + Xn binomiaal verdeeld met parameters n en p. Voorbeeld 8.4 (Binomiale verdeling). De stochastische grootheid X heeft de binomiale verdeling met parameters n ∈ N en p ∈ [0, 1] als P(X = k) =

n k p (1 − p)n−k k

voor k ∈ {0, 1, . . ., n}. De verwachtingswaarde en variantie zijn in dat geval gelijk aan EX = np en Var(X) = np(1 − p). De binomiale verdeling met parameters n = 1 en p ∈ [0, 1] is gelijk aan de Bernoulli-verdeling met parameter p. Als X1 en X2 twee onafhankelijke binomiaal verdeelde stochasten zijn met respectievelijk parameters (n, p) en (m, p), dan is X1 + X2 weer binomiaal verdeeld, met parameters (n + m, p).

8.3: Standaard verdelingen

75

Voorbeeld 8.5 (Poisson-verdeling). De stochastische grootheid X heeft de Poisson-verdeling met parameter λ > 0 als λk e−λ P(X = k) = k! voor k ∈ {0, 1, . . .}. De verwachtingswaarde en variantie zijn in dat geval gelijk aan EX = λ en Var(X) = λ. Als X1 en X2 twee onafhankelijke Poisson-verdeelde stochasten zijn met respectievelijk parameters λ en µ dan is X1 + X2 weer Poisson-verdeeld, met parameter λ + µ.

Voorbeeld 8.6 (Geometrische verdeling). De stochastische grootheid X heeft de geometrische verdeling met parameter p ∈ (0, 1] als P(X = k) = p(1 − p)k−1 voor k ∈ {1, 2, . . .}. De verwachtingswaarde en variantie zijn in dat geval gelijk aan EX = 1/p en Var(X) = (1−p)/p2. Als X1 , . . ., Xr onderling onafhankelijke geometrisch verdeelde stochasten zijn met parameter p, dan is X1 + . . . + Xr negatief binomiaal verdeeld met parameters r en p.

8.3.2

Continue verdelingen

Voorbeeld 8.7 (Uniforme verdeling). De stochastische grootheid X heeft de (continue) uniforme verdeling op het interval [a, b] als de dichtheid van X gelijk is aan f (x) =

1 1[a,b] (x). b−a

De waarde van de indicatorfunctie 1[a,b] (x) = 1 als x ∈ [a, b] en 0 anders. De verwachtingswaarde en variantie worden in dat geval gegeven door EX = (a + b)/2 en Var(X) = (b − a) 2 /12. In het geval dat a = 0 en b = 1 is de dichtheid gelijk aan f (x) = 1[0,1] (x), de verwachting gelijk aan 1/2 en de variantie gelijk aan 1/12. Voorbeeld 8.8 (Normale verdeling). De stochastische grootheid X heeft de normale verdeling met parameters µ ∈ R en σ 2 > 0 als de dichtheid van X gelijk is aan f (x) = √

1 2πσ 2

2 1 (x−µ) σ2

e− 2

.

De verwachtingswaarde en variantie worden in dat geval gegeven door EX = µ en Var(X) = σ 2 . De standaard normale verdeling is de normale verdeling met parameters µ = 0 en σ 2 = 1. De dichtheid en de verdelingsfunctie van de standaard normale verdeling worden respectievelijk genoteerd als φ en Φ. Als X1 en X2 twee onafhankelijke normaal verdeelde stochasten zijn met respectievelijk parameters (µ, σ 2 ) en (ν, τ 2 ) dan is X1 + X2 weer normaal verdeeld, met parameters (µ + ν, σ 2 + τ 2 ). Voorbeeld 8.9 (Exponenti¨ ele verdeling). De stochastische grootheid X heeft de exponentiële verdeling met parameter λ > 0 als de dichtheid van X gelijk is aan f (x) = λe−λx ,

x ≥ 0.

De verwachtingswaarde en variantie worden in dat geval gegeven door EX = 1/λ en Var(X) = 1/λ2 . Als X1 , . . ., Xn onderling onafhankelijke exponentieel verdeelde stochasten zijn met parameter λ, dan is de som X1 + . . . + Xn Gamma-verdeeld met vormparameter n en inverse schaalparameter λ.

76


Voorbeeld 8.10 (Gamma-verdeling). De stochastische grootheid X heeft de Gammaverdeling met vormparameter α > 0 en inverse schaalparameter λ > 0 (of schaalparameter 1/λ) als de dichtheid van X gelijk is aan f (x) =

xα−1 λα e−λx , Γ(α)

x ≥ 0,

waar Γ de zogenaamde Gamma-functie is, Γ(α) =

Z

∞

xα−1 e−x dx.

0

Wanneer α ∈ N, dan geldt Γ(α) = (α − 1)!. De verwachtingswaarde en variantie van X worden gegeven door EX = α/λ en Var(X) = α/λ2 . De Gamma-verdeling met parameters α = 1 en λ > 0 is gelijk aan de exponentiële verdeling met parameter λ. Als X1 en X2 twee onafhankelijke Gamma-verdeelde stochasten zijn met respectievelijk parameters (α, λ) en (β, λ) dan is X 1 +X2 weer Gamma-verdeeld, met parameters α + β en λ. Voorbeeld 8.11 (Chikwadraat-verdeling). De stochastische grootheid de chikwadraatPnX heeft 2 verdeling met n vrijheidsgraden als X dezelfde verdeling heeft als Z voor Z1 , . . ., Zn i=1 i onderling onafhankelijke standaard normaal verdeelde stochasten. De verwachtingswaarde en variantie van X worden gegeven door EX = n en Var(X) = 2n. De chikwadraat-verdeling met n vrijheidsgraden wordt genoteerd als χ2n . Voorbeeld 8.12 (t-verdeling). De stochastische grootheid X bezit de t-verdeling (of Studentp verdeling) met n vrijheidsgraden, als X dezelfde verdeling heeft als Z/ Y /n waar Y en Z onafhankelijke stochastische grootheden zijn, Z een standaard normaal verdeling volgt en Y een χ2n -verdeelde stochast is. De t-verdeling met n vrijheidsgraden wordt genoteerd als t n . Voorbeeld 8.13 (F -verdeling). Een stochastische grootheid X bezit de F -verdeling met m en n vrijheidsgraden, als X dezelfde verdeling heeft als (U/m)/(V /n) waar U en V onafhankelijke chikwadraat-verdeelde stochastische grootheden zijn met respectievelijk m en n vrijheidsgraden. De F -verdeling met m en n vrijheidsgraden wordt genoteerd als Fm,n .

8.4

Multivariate en marginale verdelingen

In veel gevallen is men niet ge¨ınteresseerd in een enkele waarneming, maar wil men meerdere gemeten grootheden tegelijk beschouwen. In de kansrekening worden voor dergelijke situaties stochastische vectoren gebruikt. Het eenvoudigste geval is dat van twee stochasten X en Y die samengevoegd worden tot een vector (X, Y ). Noteer de uitkomstenruimtes van X en Y als Ω X en ΩY respectievelijk. De mogelijke uitkomsten van (X, Y ) zijn punten (x, y) ∈ Ω = ΩX × ΩY . Wanneer X en Y reëelwaardig zijn is de gezamenlijke uitkomstenruimte Ω gelijk aan (een deel van) het vlak, dat is Ω ⊆ R2 . De simultane verdeling van X en Y beschrijft de kansen van de vorm P (X, Y ) ∈ A , waarbij A een deelverzameling van Ω is. Ook bij stochastische vectoren is er onderscheid tussen discrete en continue verdelingen. Als de vector (X, Y ) een discrete verdeling heeft, dan ligt die verdeling vast door de simultane kansdichtheid p(ω1 , ω2 ), waar p(ω1 , ω2 ) = P (X, Y ) = (ω1 , ω2 ) , voor alle mogelijke uitkomsten (ω1 , ω2 ) ∈ Ω. In dat geval is de kans op een uitkomst binnen een deelverzameling A van Ω gelijk aan de som X P (X, Y ) ∈ A = p(ω1 , ω2 ). (ω1 ,ω2 )∈A

8.5: Onafhankelijkheid en conditionering

77

Wanneer de vector (X, Y ) continu verdeeld is, hanteren we een multivariate kansdichtheidsfunctie f : R2 → R, kortweg ook multivariate kansdichtheid genoemd. De kans op een uitkomst binnen een verzameling A is in dat geval gegeven door de integraal Z P (X, Y ) ∈ A = f (x, y) dx dy. A

Voorbeelden van multivariate dichtheden en berekeningen daarmee zijn te vinden in de tekstboeken, zoals eerder genoemd. Uit de multivariate verdeling van de stochastische vector (X, Y ) kunnen de marginale verdelingen van X en Y verkregen worden. Een marginale verdeling wordt vastgelegd door de bijbehorende marginale kansdichtheid. In het geval van discrete stochasten wordt de marginale kansdichtheid van X, pX , als volgt berekend uit de multivariate kansdichtheid, p: X pX (ω) = p(ω, ω2 ). ω2 ∈ΩY

Voor continue stochasten wordt de sommatie vervangen door integratie, en wordt de marginale kansdichtheidsfunctie fX van X berekend uit de multivariate kansdichtheid f , Z ∞ fX (x) = f (x, y) dy. −∞

Voor de marginale verdeling van Y gelden soortgelijke formules. Voor stochastische vectoren (X1 , . . ., Xn ) met n > 2 laat het bovenstaande zich gemakkelijk generaliseren. In dat geval wordt de marginale dichtheid van X1 bijvoorbeeld gevonden door de multivariate dichtheid te integreren over x2 , . . ., xn .

8.5

Onafhankelijkheid en conditionering

Onafhankelijkheid en conditionering van stochastische grootheden speelt binnen de statistiek een grote rol. Vanwege de analogie wordt een aantal definities en belangrijke stellingen omtrent onafhankelijke en conditionele eventualiteiten hieronder eveneens gegeven. Definitie 8.14. Twee stochastische grootheden X en Y heten onafhankelijk wanneer voor alle eventualiteiten A ⊆ ΩX en B ⊆ ΩY geldt P X ∈ A, Y ∈ B = P(X ∈ A) P(Y ∈ B).

De volgende stelling laat zien hoe onafhankelijkheid van twee stochasten in de simultane verdeling tot uiting komt. Voor het bewijs van deze stelling en van Stelling 8.17 wordt verwezen naar de tekstboeken. Stelling 8.15. Als de stochasten X en Y een discrete simultane verdeling hebben met kansdichtheid p, dan zijn X en Y onafhankelijk dan en slechts dan als p(ω1 , ω2 ) = pX (ω1 ) pY (ω2 ),

voor alle ω1 , ω2 .

Als (X, Y ) een continue simultane verdeling hebben met kansdichtheid f , dan zijn X en Y onafhankelijk dan en slechts dan als f (x, y) = fX (x) fY (y),

voor alle x, y.

Onafhankelijkheid van X en Y betekent dat informatie over de realisatie van Y geen invloed heeft op de verdeling van X en vice versa. Deze heuristische interpretatie kan worden onderbouwd door voorwaardelijke kansen te beschouwen.

78


Definitie 8.16. Voor de stochasten X en Y is de voorwaardelijke kans op X ∈ A gegeven Y ∈ B gedefinieerd door P X ∈ A, Y ∈ B P X∈A Y ∈B = , P Y ∈B voor A ⊆ ΩX , B ⊆ ΩY en P Y ∈ B > 0.

Stelling 8.17 (Regel vanSBayes). Veronderstel dat A1 , . . ., An een partitie is van Ω, ofwel n Ai ∩ Aj = ∅ voor i 6= j en i=1 Ai = Ω, en dat P(X ∈ Ai ) > 0 voor alle i. Dan geldt voor een willekeurige eventualiteit B met P(X ∈ B) > 0 X ∈ Ai P X ∈ Ai P X ∈ B . P X ∈ A i X ∈ B = Pn j=1 P X ∈ B X ∈ Aj P X ∈ Aj

Voor onafhankelijke stochastische grootheden X en Y vereenvoudigt de kans op X ∈ A gegeven Y ∈ B tot P X∈A Y ∈B =P X∈A ,

voor alle A ⊆ ΩX en B ⊆ ΩY , omdat de productvorm uit Definitie 8.14 in dat geval geldt. Deze berekening laat zien dat, als het gebruik van extra informatie over Y door conditionering de verdeling van X niet be¨ınvloedt, dan zijn X en Y onafhankelijk.

Definitie 8.18. Voor continu verdeelde stochasten X en Y is de voorwaardelijke dichtheid van X gegeven Y gelijk aan f (x, y) . fX|Y (x|y) = fY (y) Ook hier geldt dat als X en Y onafhankelijk zijn, de voorwaardelijke dichtheid van X gegeven Y vereenvoudigt: fX|Y (x|y) = fX (x) voor willekeurige y. De verwachtingswaarde en variantie van de som van twee onafhankelijke stochastische grootheden zijn gelijk aan E(X + Y ) = EX + EY var(X + Y ) = var(X) + var(Y ) omdat in dat geval Cov(X, Y ) = E(XY ) − EXEY = 0. Bovenstaande uitdrukkingen zijn eenvoudig uit te breiden naar sommen van n stochastische grootheden. Veronderstel dat X1 , . . ., Xn stochastische grootheden zijn met eindige verwachting µ, dan is n n X X EXi = nµ. Xi = E i=1

i=1

2

Als X1 , . . ., Xn eindige variantie σ hebben én onderling onafhankelijk zijn, dan geldt tevens dat n n X X var Xi = nσ 2 . var Xi = i=1

i=1

De verwachting en de variantie van het steekproefgemiddelde n

X= zijn dan gelijk aan

1X Xi n i=1

1 X E Xi = µ, n i=1 n

EX =

X σ 2 1 var Xi = . n2 n i=1 n

var X =

8.6: Limietstellingen en de normale benadering

8.6

79

Limietstellingen en de normale benadering

In het geval van onderling onafhankelijke en identiek continu verdeelde stochasten X 1 , X2 , . . ., Xn met marginale kansdichtheid fX , is de simultane kansdichtheid f gelijk aan het product van de marginale dichtheden n Y fX (xi ). f (x1 , . . ., xn ) = i=1

Deze simultane dichtheid komt veelvuldig voor in statistische vraagstukken. De volgende belangrijke stellingen uit de kansrekening zijn op rijen van onderling onafhankelijke en identiek verdeelde stochasten van toepassing. Aangezien de limietstellingen de limiet voor n → ∞ betreffen, wordt het steekproefgemiddelde X in deze paragraaf genoteerd met X n waarin de afhankelijkheid van n expliciet naar voren komt.

Stelling 8.19. (Zwakke wet van de grote aantallen) Stel dat X1 , X2 , . . . onderling onafhankelijk en identiek verdeeld zijn volgens een marginale verdeling met eindige verwachtingswaarde µ. Dan geldt voor iedere ε > 0 lim P |X n − µ| > ε = 0. n→∞

Stelling 8.20. (Sterke wet van de grote aantallen) Stel dat X1 , X2 , . . . onderling onafhankelijk en identiek verdeeld zijn volgens een marginale verdeling met eindige verwachtingswaarde µ. Dan geldt P lim X n = µ = 1. n→∞

Stelling 8.21. (Centrale Limietstelling) Stel dat X1 , X2 , . . . onderling onafhankelijk en identiek verdeeld zijn volgens een marginale verdeling met eindige verwachtingswaarde µ en eindige variantie σ 2 . Dan geldt √n(X − µ) n √ lim P ≤ z = Φ(z), n→∞ σ2 waar Φ de verdelingsfunctie van de standaard normale verdeling is.

De Centrale Limietstelling als zodanig kan niet in de praktijk worden toegepast omdat we nooit over oneindig veel data beschikken. Voor grote waarden van n is de kans in het linkerlid echter bij benadering gelijk aan de kans in het rechterlid. Hoe groot n precies moet zijn voor een redelijke benadering is onder andere afhankelijk van de scheefheid van de marginale verdeling. De stochastische grootheid in het linkerlid kan op een andere manier worden geschreven √ n(X n − µ) Xn − µ X n − EX n √ = p = p . 2 2 σ /n σ var X n

Met de Centrale Limietstelling volgt dus dat het gestandaardiseerde steekproefgemiddelde bij benadering de standaard normale verdeling volgt als het aantal waarnemingen groot is.

Voorbeeld 8.22 (Normale benadering van de binomiale verdeling). Veronderstel dat X1 , . . ., Xn een steekproef is uit de Bernoulli-verdeling met parameter p. De corresponderende verwachtingswaarde en variantie zijn beide eindig, respectievelijk p en p(1−p). Met de Centrale Limietstelling volgt voor grote waarden van n, dat X n − EX n Xn − p p =p p(1 − p)/n var X n

bij benadering standaard normaal verdeeld is. Het steekproefgemiddelde is dus bij benadering N (p, p(1 − p)/n) verdeeld. Hieruit kunnen we tevens een benadering voor de verdeling van Y = P n i=1 Xi , dat is de binomiale verdeling met parameters n en p, afleiden. Als X n bij benadering de N (p, p(1 − p)/n)-verdeling volgt, dan volgt Y = nX n bij benadering de N (np, np(1 − p))verdeling. Deze benadering is redelijk als n niet te klein is en p niet te dicht bij 0 of 1 in de

80


buurt ligt. Als vuistregel hanteert men de eis np(1 − p) ≥ 5. Aangezien de binomiale verdeling discreet is en de normale verdeling continu, past men in de regel een continu¨ıteitscorrectie toe om de benadering te verbeteren. De kansmassa op Y = i in de discrete verdeling wordt als het ware uitgesmeerd over het interval (i − 1/2, i + 1/2] in de continue verdeling, P(Y = i) = P(i − 1/2 < Y ≤ i + 1/2), voor alle i ∈ N. Dit levert P Y ≤ i = P Y ≤ i + 21 en P Y ≥ i = P Y > i − 21 .

De combinatie van de normale benadering en de continu¨ıteitscorrectie levert bijvoorbeeld voor n = 25 en p = 0.4 1.5 Y − 10 11.5 − 10 √ √ ≈ Φ √ = 0.730. P(Y ≤ 11) = P(Y ≤ 11.5) = P ≤ 6 6 6 De exacte kans is in dit geval gelijk√ aan 0.732. Ter vergelijking geven we ook de benaderde kans zonder continu¨ıteitscorrectie: Φ(1/ 6) = 0.658. De correctie zorgt hier dus duidelijk voor een verbetering van de benadering.

Voorbeeld 8.23 (Normale benadering van de Poisson-verdeling). Veronderstel dat X1 , . . ., Xn een steekproef is uit de Poisson-verdeling met parameter λ. De verwachtingswaarde en de variantie van deze verdeling zijn beide gelijk aan λ en dus eindig. De toepassing van de Centrale Limietstelling geeft nu dat X n − EX n Xn − λ p = p λ/n var X n

bij benadering de standaard normale verdeling volgt. Deze grootheid is ook te schrijven als Pn Xi − nλ i=1 √ . nλ Omdat de som van onafhankelijke Poisson-verdeeldePstochasten weer Poisson-verdeeld is, zie n Voorbeeld 8.5, is de stochastische grootheid Y = i=1 Xi Poisson-verdeeld met parameter µ = nλ. Er volgt dus dat Y −µ √ µ bij benadering de standaard normale verdeling volgt. Dit is equivalent met te zeggen dat de Poisson(µ)-verdeling te benaderen is met de N (µ, µ)-verdeling voor grote waarden van µ.

Opgaven 1. Bereken EX 2 als X Poisson-verdeeld is met parameter θ. 2. Veronderstel dat X en Y onderling onafhankelijk en exponentieel verdeeld zijn met verwachting 1. Bepaal de kansdichtheid en de verwachting van max(X, Y ). 3. Veronderstel dat X en Y onafhankelijk en N (0, 1) en N (1, 2)-verdeeld zijn. (i) Bepaal P(X + Y ≤ 2). (ii) Bepaal een getal ξ zodanig dat P(X + Y > ξ) = 0.95. 4. Veronderstel dat X1 , . . ., Xn onderling onafhankelijk en homogeen verdeeld zijn op het interval [0, 1]. Bepaal de verwachting en variantie van Y = max(X1 , . . ., Xn ). Hint: leid de dichtheid van Y af uit de verdelingsfunctie van Y , P(Y ≤ y), die gevonden kan worden met behulp van de verdelingsfuncties van X1 , . . ., Xn . 5. Veronderstel dat X1 , . . ., Xn onderling onafhankelijk en identiek verdeeld zijn met verwachting µ en variantie σ 2 . Bereken EX n , var X n , E(X n )2 en cov(Xi − X n , X n ).

6. Zij X binomiaal verdeeld met parameters 100 en 1/4. Bepaal een benadering voor P(X ≤ 30) met behulp van de Centrale Limietstelling.

7. Veronderstel dat X1 , . . ., X25 onderling onafhankelijk en Poisson(5) verdeeld zijn. Bepaal een benadering voor P(X n ≥ 4.5) met behulp van de Centrale Limietstelling.

ALGEMENE STATISTIEK. A.W. van der Vaart en anderen

Recommend Documents