Goodness-of-fit bij linksgecensureerde observaties. Sara Appeltants

Universiteit Antwerpen Faculteit Wetenschappen Departement Wiskunde-Informatica Academiejaar 2005-2006

Goodness-of-fit bij linksgecensureerde observaties

— Sara Appeltants —

Eindwerk ingediend met het oog op het behalen van de graad van Licentiaat in de Wiskunde Promotor: Prof. Dr. G. Molenberghs Copromotor: Prof. Dr. R. Braekers

Dankwoord Als eerste wil ik mijn promotor Geert Molenberghs en copromotor Roel Braekers bedanken omdat ze mij de kans gegeven hebben om deze thesis te maken. Roel Braekers wil ik nog extra bedanken voor de goede begeleiding, de snelle antwoorden op mijn mails met vragen en het vele geduld. Ook mijn ouders verdienen mijn dankbaarheid omdat ze altijd voor mij klaarstonden en me zijn blijven steunen, op alle mogelijke manieren. Natuurlijk wil ik ook mijn vrienden en medestudenten bedanken voor allerlei kleine en grote dingen. In het bijzonder wil ik Claudia en Carolien bedanken omdat ze steeds bereid waren om naar mij te luisteren en hun bemoedigende woorden, Mario voor het nalezen van mijn thesis, en ook Tamara voor de hulp met TeX. Tenslotte ben ik Peter erg dankbaar voor diverse hulp, van het uitleggen van foutmeldingen in mijn algoritmes tot het fungeren als publiek bij het oefenen van mijn presentatie. Maar vooral omdat hij steeds in mij is blijven geloven.

i

Inhoudsopgave Inleiding

1

1 Inleidende begrippen tot survivalanalyse 1.1 Survivalfunctie . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Dichtheid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Hazard rate functie . . . . . . . . . . . . . . . . . . . . . . . .

2 2 3 4

2 Gecensureerde gegevens 2.1 Rechtse censurering . . . . . . . . . . . . 2.1.1 Type I: vaste censureringstijd . . 2.1.2 Type I: random censureringstijd 2.1.3 Type II censurering . . . . . . . 2.2 Linkse censurering . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3 Goodness-of-fit tests 3.1 Likelihoodconstructie . . . . . . . . . . . . . . . . 3.2 Chi-kwadraat test . . . . . . . . . . . . . . . . . 3.2.1 Chi-kwadraat test voor niet gecensureerde 3.2.2 Chi-kwadraat test voor linksgecensureerde 3.3 QQ-plot . . . . . . . . . . . . . . . . . . . . . . . 3.4 EDF-statistieken . . . . . . . . . . . . . . . . . . 3.4.1 Empirische distributiefunctie . . . . . . . 3.4.2 Supremumstatistieken . . . . . . . . . . . 3.4.3 Kwadratische statistieken . . . . . . . . . 3.4.4 Berekenen van de statistieken: PIT . . . . 3.4.5 EDF-tests: case 0 . . . . . . . . . . . . . 3.4.6 EDF-tests met ongekende parameters . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . gegevens gegevens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .

8 8 8 10 13 15

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

17 18 19 19 23 27 31 31 32 33 34 37 45

Besluit

48

Bibliografie

50

A Enkele definities en eigenschappen

51

ii

B Algoritmes B.1 Kunstmatig censureren van simulaties . . . . . . . . . . . B.1.1 Linkse censurering . . . . . . . . . . . . . . . . . . B.1.2 Rechtse censurering . . . . . . . . . . . . . . . . . B.2 Construeren van voor simulaties . . . . . . . . . . . . . B.2.1 Linkse censurering . . . . . . . . . . . . . . . . . . B.2.2 Rechtse censurering . . . . . . . . . . . . . . . . . B.3 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3.1 likelihoodconstructie . . . . . . . . . . . . . . . . . B.3.2 Contourplot van de loglikelihoodfunctie . . . . . . B.4 Chi-kwadraat verdelingstest . . . . . . . . . . . . . . . . . B.4.1 Chi-kwadraattest voor ongecensureerde gegevens . B.4.2 Chi-kwadraattest voor linksgecensureerde gegevens B.5 QQ-plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.6 EDF-statistieken case 0 . . . . . . . . . . . . . . . . . . . B.6.1 Grafieken van de EDF en verdelingsfunctie . . . . B.6.2 EDF-testen voor ongecensensureerde gegevens . . . B.6.3 EDF-testen voor rechtsgecensensureerde gegevens . B.6.4 EDF-testen voor linksgecensensureerde gegevens . B.7 EDF-statistieken case 1 voor ongecensureerde gegevens . . B.8 EDF-statistieken case 2 voor ongecensureerde gegevens . . B.9 EDF-statistieken case 3 voor ongecensureerde gegevens . . B.10 Gebruikte data . . . . . . . . . . . . . . . . . . . . . . . .

iii

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

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

53 53 53 53 54 54 54 54 54 55 55 55 57 59 60 60 61 62 63 64 65 65 66

Inleiding De analyse van overlevingstijd, time to event data komt veelvuldig voor in de toegepaste statistiek. Voorbeelden hiervan kunnen we vinden in de medische wereld, biologie, psychologie, economie en levensduur van materialen, nl. tijd tot reageren op een geneesmiddel, tijd tot de dood van een laboratoriummuis, tijd tot het falen van een toestel e.a. Een probleem dat deze analyse nog bemoeilijkt, is het ontbreken van data, in het bijzonder censurering. Met gecensureerde data bedoelt men dat er naast de exacte ook gecensureerde observaties voorkomen. Men spreekt van een gecensureerde observatie als de lifetime van een bepaald individu of object (ook wel event time genoemd) niet weet wanneer dit exact plaatsheeft, maar enkel in welke tijdsperiode. Hierdoor zullen de klassieke statistische schatters en goodness-of-fit tests niet meer volstaan. Denk maar aan het steekproefgemiddelde als schatter voor het populatiegemiddelde, of de Chi-kwadraat verdelingstest. Na een beknopte inleiding over survivalanalyse zal er aan de hand van talrijke voorbeelden een overzicht gegeven worden van de verschillende soorten censurering. Vooral het verschil tussen rechtse en linkse censurering zal aan bod komen. Men spreekt van rechtsgecensureerde data wanneer we enkel weten dat dat de gebeurtenis heeft plaatsgehad voor een bepaalde tijd, bvb. wanneer een patiënt nog niet gereageerd heeft op het medicijn voor het einde van de studie. Bij linksgecensureerde data zal de gebeurtenis al plaatsgehad hebben voor een bepaalde tijd, bvb. wanneer de event time de tijd tot de eerste stapjes van een kind zijn, zijn de kinderen die al kunnen lopen linksgecensureerd. Ook het onderscheid tussen type I en type II bij rechtse censurering zal even aan bod komen. In het volgende deel zullen we zien hoe de likelihoodfunctie verandert wanneer we censurering in rekening brengen. Deze zullen we gebruiken in enkele van de latere verdelingstoetsen. We zullen de Chi-kwadraat test, de QQ-plot en verdelingstoetsen gebaseerd op EDF-statistieken uitvoerig bestuderen waarna we ze aanpassen voor linkse censurering met vaste censuringstijd, soms via een omweg langs de rechtste variant. Hiervoor heb ik telkens een algoritme geschreven in het programma R dat hier wordt uitgevoerd op enkele simulaties en voorbeelden. 1

Hoofdstuk 1

Inleidende begrippen tot survivalanalyse Zij X de tijd tot een bepaalde gebeurtenis. De stochastische veranderlijke X zal steeds positief zijn. Dit kan de tijd tot het falen van een elektrisch toestel zijn, tijd tot het hervallen van een patiënt met leukemie, sterftetijd van longkankerpatiënten, ... Event times worden meestal geassocieerd aan negatieve gebeurtenissen. In het Nederlands spreekt men dan soms ook van faaltijd. Maar er zijn ook voorbeelden zoals een onderzoek in een kinderdagverblijf waarbij men de tijd tot het leren van bepaalde vaardigheden, zoals stappen of spreken, waarneemt. We zullen de termen event time, lifetime en faaltijd door elkaar gebruiken. In dit hoofdstuk zullen we begrippen invoeren zoals dichtheid, survivalfunctie en hazard rate functie, die de event time volledig karakteriseren.

1.1

Survivalfunctie

De survivalfunctie geeft de kans weer dat de event time groter is dan een tijd x, dat het individu nog leeft na tijd x of de gebeurenis plaatsheeft na tijd x. We definiëren deze functie als: S(x) = P (X > x)

2

Figuur 1.1: Survival functie van een Weibull verdeling De survivalfunctie heeft nog volgende interessante eigenschappen. Eigenschap 1.1.1 Voor een continue stochastische veranderlijke X heeft de survivalfunctie volgende eigenschappen: 1. S(x) is rechtscontinu en niet-stijgend monotoon. 2. S(0) = 1 3. limx!1 S(x) = 0 4. S(x) = 1

P (X  x) = 1

F (x), met F (x) de verdelingsfunctie.

Figuur 1.1 geeft een voorbeeld van een survivalfunctie, nl. die van de Weibullverdeling met parameters = 1 en ↵ = 3. De helling van de survivalfunctie wordt bepaald door het risico dat het individu of object loopt: hoe meer risico, hoe sneller S(x) naar nul zal gaan.

1.2

Dichtheid

De dichtheid van een stochastische variabele definieert men als dS(x) dx

f (x) = waarbij

dS(x) dx

de afgeleide van S(x) naar x is.

Eigenschap 1.2.1 Voor een continue stochastische veranderlijke X heeft de dichtheidsfunctie volgende eigenschappen: 3

1. f (x) =

dF (x) dx

2. f (x) 0, voor alle x R1 3. 1 f (x)dx = 1 Rb 4. P (a  X  b) = a f (t)dt = F (b)

F (a) = S(a)

S(b)

Voor de dichtheid hebben we ook volgende verband met de survivalfunctie en verdelingsfunctie, dat we duidelijk kunnen zien op figuur 1.2. R1 • S(t0 ) = t0 f (t)dt Rt • F (t0 ) = 01 f (t)dt

Figuur 1.2: Dichtheidsfunctie van Weibulverdeling met

= 1 en ↵ = 3

We kunnen de dichtheid interpreteren als een ’benadering’ van de kans dat de gebeurtenis plaatsheeft op tijd x: f (x) x ⇡ P (x  X  x +

1.3

x) , voor kleine

x

Hazard rate functie

De hazard rate functie, ook wel faaltempo genoemd, is een populaire beschrijving van een lifetime X omdat ze het risico van het individu of object op tijdstip t voorstelt. We kunnen de definitie voor continue stochastische variabelen eenvoudig afleiden uit deze intu¨ıtieve voorstelling: 1. kans op gebeurtenis tussen tijd t en t + = P (t  T  t +

t) = S(t)

S(t +

wegens eigenschap 1.2.1 4

t)

t

2. kans op gebeurtenis tussen tijd t en t + op tijd t = P (t  T  t +

t) =

t|T

t, gegeven nog niets gebeurd

P (tT t+ t) P (T t)

=

S(t) S(t+ t) S(t)

wegens de definitie van voorwaardelijke kans en S survivalfuntie van een continue stochastische variabele 3. gemiddelde kans op gebeurtenis tussen tijd t en t + niets gebeurd op tijd t: =

P (tT t+ t|T t) t

=

S(t) S(t+ t) S(t)· t

4. het ogenblikkelijk gemiddelde (lim = lim

t!0

t, gegeven nog

P (tT t+ t|T t) t

= lim

t!0 ) t!0

:

S(t) S(t+ t) S(t)· t

=

dS(t)/dt S(t)

=

f (t) S(t)

wegens de definitie van afgeleide Men definieert hazard rate functie dan als P (x  X  x + x | X x!0 x

h(x) = lim

x)

Voor een continue stochastische variabele X geldt dan ook: h(x) = f (x)/S(x) =

d(ln[S(x)]) dx

Dankzij deze gelijkheid kunnen we kiezen of we de survivalfunctie of de hazard rate functie gebruiken om de event time van een individu of object te beschrijven, want beide zijn equivalent. Wanneer we H(x) definiëren als de cumulatieve hazard rate functie: Z x H(x) = h(u)du = ln [S(x)] 0

kunnen we voor continue survivalfuncties S(x) schrijven in functie van de (cumulatieve) hazard rate functie. S(x) = e

H(x)

=e

Rx 0

h(u)du

(1.1)

Hazard rate functies kunnen verscheidene vormen hebben. De meest eenvoudige hazard rate functie kunnen we vinden bij een exponentiële verdeling. Dit is de constante functie h(x) = waarbij > 0 de parameter van deze verdeling is. Hieruit volgt dat wanneer X een exponentiële verdeling volgt, het risico op ieder moment evenveel is (onafhankelijk van de tijd). Voor stijgende hazard rate functies geldt: hoe verder in de tijd, hoe meer risico. Zo’n model kan gebruikt worden bij het verouderen van personen of 5

slijtage van materialen. Voor dalende hazard rate functies geldt het omgekeerde, hoe verder in de tijd, hoe minder risico. Deze kunnen model staan voor de inlooptijd van een nieuwe werknemer of metalen die sterker worden door gebruik.

Figuur 1.3: Enkele voorbeelden van monotone hazard rate functies Maar er zijn ook niet-monotone hazard rate functies. Zoals de hazard rate functie van de log-logistische verdeling met parameters ↵ = 2 en = 0.135. Deze heeft eerst een stijgend verloop en daarna dalend. Dit kan een model zijn voor een risicovolle operatie waarbij in het begin kans is op complicaties en infectie, maar wanneer de patiënt herstellende is, zal het risico dalen. Tenslotte, de meest voorkomende hazard rate functie is badkuipvorm. Deze wordt gebruikt in experimenten waar men de populatie vanaf de geboorte volgt. Het hoge risico in het begin staat voor kindersterfte, tijdens de kinderjaren is het risico het laagst, maar naarmate men ouder wordt, zal het risico terug stijgen. Voor de badkuipcurve in figuur 1.4 gebruikten we de formule van Lee (zie A7), maar een badkuipcurve kan ook gemodelleerd worden als een stuksgewijze set van een dalende, een constante en een stijgende hazard rate functie [9].

6

Figuur 1.4: Enkele voorbeelden van niet-monotone hazard rate functies

7

Hoofdstuk 2

Gecensureerde gegevens In een steekproef van lifetimes kunnen naast exacte observaties ook een aantal gecensureerde observaties voorkomen. Een gecensureerde observatie bevat slechts gedeeltelijke informatie over de lifetime X. We weten enkel in welk interval deze zal voorkomen, niet de exacte waarde. Bijvoorbeeld: voor een longkankerpatiënt die sterft in een auto-ongeval kunnen we niet meten wanneer hij aan longkanker zal sterven, we weten enkel dat de exacte tijd hiervan later zou zijn dan het moment van het ongeval. We onderscheiden rechtse en linkse censurering, afhankelijk van of dit interval langs boven of langs onder begrensd is. Binnen deze categoriën kunnen we nog verder opdelen in Type I censurering vaste en random censureringstijd en bij rechtse censurering ook Type II.

2.1 2.1.1

Rechtse censurering Type I: vaste censureringstijd

Bij rechtse censurering wordt de gebeurtenis enkel exact geobserveerd als ze plaatsheeft voor een bepaalde tijd. Voorbeeld 2.1.1 Een onderzoek van 100 patiënten die een behandeling tegen leukemie ondergaan. Om tijd en kosten te besparen wordt de studie na 6 maanden stopgezet. De patiënten die nog niet op de behandeling gereageerd hebben, geven rechtsgecensureerde observaties. Voorbeeld 2.1.2 Een producent van gloeilampen stelt een onderzoek naar de levensduur van lampen zodat hij een garantie op de verpakking kan plaatsen. Hij eist dat de studie niet langer dan 3600 uur mag duren.

8

Figuur 2.1: Schematische voorstelling rechtse censurering: de observaties voor tc zijn exacte lifetimes Zij lifetimes X1 , X2 , .., Xn onafhankelijk en gelijk verdeeld met verdeling F en tc een (voorafbepaalde) tijd, die we vaste censureringstijd noemen. We zijn ge¨ınteresseerd in X1 , X2 , .., Xn , maar deze observeren we echter niet. In plaats daarvan observeren we T1 , T2 , .., Tn met ( Xi als Xi  tc Ti = . tc als Xi > tc We kunnen dit ook noteren met (Ti , i ) waarbij Ti = min(Xi , tc ) en i = 1(Ti Ci ) . Indicator i , die de informatie over de censurering bevat, kunnen we ook noteren als: 8 > 0 als gecensureerd > > > < als X > t , i c . i = >1 als niet gecensureerd > > > : als X  t i c

Figuur 2.2: Een geobserveerde survivalfunctie

9

Van deze geobserveerde verdeling kunnen we de survivalfunctie bepalen: STi (t) = = = =

P (Ti > t) P (min(Xi , tc ) > t) P (Xi > t en tc > t) P (Xi > t) · 1(t
definitie S definitie Ti

Een tekening hiervan zien we in figuur 2.1.1.

2.1.2

Type I: random censureringstijd

Tot nu toe hebben we enkel voorbeelden gezien met een vaste censureringstijd, maar dit kunnen we ook veralgemenen naar een random censureringstijd waarbij we het vast getal tc vervangen door een stochastische variabele Ci met een bepaalde verdeling G. Aan iedere lifetime Xi associëren we dan zo’n een stochastische variabele Ci (1  i  n). We kunnen voorbeeld 2.1.1 uitbreiden: Voorbeeld 2.1.3 Een onderzoek van 100 patiënten die een behandeling tegen leukemie ondergaan. Om tijd en kosten te besparen wordt de studie na 6 maanden stopgezet. Maar nu laten we ook andere factoren toe die censurering veroorzaken: het kwijtraken van gegevens, terugtrekkingen uit de studie (bvb. wegens teveel nevene↵ecten van de behandeling), ’loss to follow-up’ (bvb. door verhuis van de patiënt, gestorven door andere oorzaken dan diegene waarin we ge¨ınteresseerd zijn), ... We kunnen ook nog een extra veralgemening beschouwen die toelaat dat de individuen op verschillende tijden de studie binnenkomen en dus niet enkel op de starttijd van de studie t = 0. Dit komt vaak voor bij klinische studies waar men niet enkel de patiënten gebruikt die voor handen zijn bij het begin, maar men ook in de loop van de studie patiënten toevoegt. Bijvoorbeeld in een onderzoek met leukemiepatiënten neemt men patiënten op waarvan de diagnose pas tijdens de duur van de studie gesteld werd. Om deze observaties gelijk te behandelen zal men niet de kalendertijd gebruiken, maar herschalen naar tijd in de studie. Figuur 2.1.2 toont een triviaal voorbeeld van een mogelijke deelsteekproef van deze studie. In dit geval geeft patiënt 1 een ongecensureerde observatie want hij reageert voor het einde van de studie. Patiënt 2 komt wat later binnen in de studie en geeft een gecensureerde observatie, hij heeft nog niet gereageerd voor het einde van de studie maar verdween ook niet wegens andere redenen 10

uit de studie. Patiënt 3 daarentegen is gecensureerd, maar bvb. doordat hij verhuisd is of niet meer wil meedoen aan het onderzoek. Merk op dat gecensureerde gegevens genoteerd worden met een + .

Figuur 2.3: Schematische voorstelling van random censurering bij vb. 2.1.3 In het algemeen observeren we bij random censurering de paren (T1 , waarbij voor alle i in {1, .., n}: Ti = min(Xi , Ci )8 > 0 > > > < i = 1(Xi Ci ) = > 1 > > > :

als gecensureerd als Xi > Ci , als niet gecensureerd als Xi  Ci

1 ), .., (Tn , n )

.

We observeren dus Xi en Ci niet rechtstreeks, maar enkel Ti en de bijbehorende indicator voor censurering i . We kunnen opmerken dat vaste censurering een speciaal geval is van random censurering want we kunnen getal tc beschouwen als een ontaarde verdeling C in tc . Merk op dat we stochastische (of random) variabelen steeds aanduiden met hoofdletters en vaste getallen of bepaalde waarden van stochastische variabelen met kleine letters.

11

Figuur 2.4: De verdelingsfunctie van de ontaarde verdeling in tc Bij random censurering maken we wel de belangrijke onderstelling dat voor alle i de lifetimes Xi en censureringstijden Ci onafhankelijk zijn, omdat we anders weinig conclusies kunnen trekken aangezien de klassieke theoriën niet meer gelden. Dit lijkt een zwakke onderstelling wanneer de observaties random gecensureerd worden. Maar in werkelijkheid zal dikwijls een relatie bestaan tussen deze stochastische variabelen. In de context van voorbeeld 2.1.3 zien we dit duidelijk. Het vertrek van een patiënt uit het onderzoek kan te maken hebben met het feit dat de behandeling geen e↵ect heeft, of dat er teveel nevene↵ecten zijn, of wanneer de patiënt ontevreden is over het onderzoek. Er zal in werkelijkheid dus wel een afhankelijkheid kunnen optreden. Laten we nu nog enkele andere voorbeelden bekijken. Voorbeeld 2.1.4 Een onderzoek bestaande uit 26 psychiatrische patiënten die tussen 1935 en 1948 opgenomen werden in universitaire ziekenhuizen in Iowa. De data (zie tabel 2.1) bestaat voor iedere patiënt uit geslacht, leeftijd bij opname en het aantal jaren follow-up. Met het aantal jaren follow-up bedoelen we de tijd dat de patiënt in de studie is, de tijd tussen binnenkomen en dood of censurering. De event time is hier de tijd tot het overlijden van de patiënt. Er is hier geen algemene stop van de studie. Uit tabel 2.1 zien we dat de event time voor de eerste patiënt exact is en gelijk aan 1 jaar, dwz. 1 jaar na opneming sterft deze patiënt. Voor de laatste patiënt in deze tabel is er een rechtsgecensureerde event time 39+, zodat we enkel weten dat hij 39 jaar na opname nog in leven is. Een van de doelen van deze studie was het vergelijken van de sterftecoefficiënt van de psychiatrische patiënten met de standaard sterftecoefficiënt van inwoners van Iowa om te bepalen of de psychiatrische patiënten een duidelijk kortere levensduur hadden.

12

Sex F F F F M M F F F F F M M

Age 51 58 55 28 21 19 25 48 47 25 31 24 25

Time of follow-up 1 1 2 22 30+ 28 32 11 14 36+ 31+ 33+ 33+

Sex F F M M M F F F M M M F M

Age 30 33 36 30 41 43 45 35 29 35 32 36 32

Time of follow-up 37+ 35+ 25 31+ 22 26 24 35+ 34+ 30+ 35 40 39+

Tabel 2.1: survival data voor psychiatrische patienten Voorbeeld 2.1.5 Een studie over beenmergtransplantaties bij leukemie. De tijd tot afstoting wordt niet geobserveerd doordat de patiënt kort na de operatie sterft of een terugval heeft voordat de behandeling e↵ect kan hebben. Voorbeeld 2.1.6 Studie bij brandwondenpatiënten waarbij tijd tot het optreden van infectie de event time is en er censureringstijden zijn zoals overlijden of ontslag uit het ziekenhuis (en dus verdwijnen uit het onderzoek) voordat er infectie waargenomen wordt.

2.1.3

Type II censurering

Men kan inzien dat het niet altijd even makkelijk is om op voorhand het tijdstip te bepalen waarop de studie moet eindigen om de kost (in tijd en geld) van de studie niet te hoog te doen oplopen en toch voldoende ongecensureerde gegevens te verkrijgen (een goede waarde voor de vaste censureringstijd tc bepalen). Daarom beschouwen we ook een andere vorm van censurering: type II. De studie duurt tot het falen (of anders gezegd: plaatshebben van de gebeurtenis) van het re individu of object, met r een voorbepaald positief getal met r < n waarbij n de grootte van de steekproef is. Nu zullen alle objecten of individuen op hetzelfde moment de studie binnenkomen (op t = 0) en de studie stopt wanneer r ervan hun exacte event time bereikt hebben. Voorbeeld 2.1.7 Bij een studie over de levensduur van 300 gloeilampen kunnen we voor de twee types verschillende experimenten opzetten: • Type I: na 3600u de studie stopzetten. 13

• Type II: studie stopzetten na het falen van 200 lampen.

In studies naar de levensduur van materialen zal men meestal een type II studie verkiezen, omdat zonder voorkennis de keuze voor de vaste censureringstijd erg slecht kan zijn. In het type I-experiment uit voorgaand voorbeeld kan het gebeuren men de faaltijd van de lampen onderschat zodat er nog geen enkele gloeilamp heeft gefaald voor het einde van de studie. Dit betekent dat alle gegevens gecensureerd zullen zijn waardoor de klassieke statistische technieken niet bruikbaar kunnen zijn bij de analyse ervan. In plaats van de werkelijke event times X1 , X2 .., Xn observeren we T(1) = X(1) T(2) = X(2) .. . T(r) = X(r) T(r+1) = X(r) .. . T(n) = X(r) waarbij X(1) < X(2) < .. < X(n) de ordestatistieken van X1 , .., Xn zijn. Merk op dat het aantal falingen r (of ook het aantal event times dat we exact observeren) en het aantal gecensureerde observaties n r gekend zijn aan de start van de studie. De censureringstijd T(r) is stochastisch. In het geval van type I was het aantal gecensureerde observaties onbekend aan de start van de studie. Aangezien we enkel de r kleinste event times overhouden is de statistische analyse van Type II gecensureerde gegevens eenvoudiger. Want hierdoor kunnen we de theorie over ordestatistieken gebruiken om survivalfuncties, likelihood e.a. te bepalen. Bij wijze van voorbeeld bepalen we de verdelingsfunctie van de observeerde gegevens. FT(k) (x) = P (T(k)  x) = P (minstens k van de T(i) ’s  x) = P (#{Ti  x} minstens k)

definitie F definitie ordestatistiek (*)

Hierin herkennen we binomiaalverdeling B(n, FT (x)): n trekkingen met succeskans P (Ti  x) = FT (x) op succes. Uitdrukking (*) kunnen we dan zien als de kans op k successen uit n. Zodat wegens Appendix A5, waar we de uitdrukking voor de kans op j successen uit n in een binomiaal verdeeld experiment vinden: FT(k) (x) =

n ✓ ◆ X n j=k

j

(FT (x))j · (1 14

FT (x))n

j

2.2

Linkse censurering

Een event time is linksgecensureerd wanneer de gebeurtenis al heeft plaatsgehad vóór de censureringstijd. Een event time X is linksgecensureerd wanneer deze kleiner is dan de random censureringstijd C, dwz. dat de gebeurtenis al plaatsgehad heeft voor tijd C. We kennen dus enkel de exacte event time wanneer deze groter of gelijk aan C is. Ook hier noteren we de gegevens waarin we ge¨ınteresseerd zijn, maar die we niet waarnemen wegens censurering, met X1 , X2 , .., Xn . We stellen deze opnieuw onafhankelijk en gelijk verdeeld volgens een verdeling F . Hieraan associëren we random censureringstijden C1 , C2 , .., Cn , onafhankelijk en gelijk verdeeld volgens een verdeling G. We gaan er weer van uit dat voor alle i de stochastische variabelen Xi en Ci onafhankelijk zijn. De geobserveerde linksgecensureerde data noteren we met paren (Ti , i ) , 1  i  n, met: ( Xi als Xi Ci Ti = = max(Xi , Ci ) Ci als Xi < Ci 8 > 0 als gecensureerd > > > < als X < C , i i (2.1) i = 1(Xi Ci ) = > 1 als niet gecensureerd > > > : als X Ci i Voorbeeld 2.2.1 In kinderdagverblijven worden regelmatig testen uitgevoerd om te bepalen wanneer een kind een bepaalde vaardigheid aanleert. Dit zijn dan vaardigheden zoals: leren lopen, eerste woordjes, zich optrekken, een cirkel kunnen tekenen, ... De event time kunnen we bijvoorbeeld de leeftijd waarop het kind leert lopen, nemen. Er zullen kinderen zijn die al kunnen lopen voor aanvang van de studie, dit zijn dan de linksgecensureerde gegevens. Want van deze kinderen weten we enkel dat hun exacte event time kleiner is dan de huidige tijd. Een kind van twee jaar dat al kan lopen zal een linksgecensureerde observatie 24+ geven (uitgedrukt in maanden). Van een kind van 10 maanden dat we tijdens de studie voor het eerst zien lopen, kennen we de exacte event time, nl. 10. Voorbeeld 2.2.2 Een onderzoek over het marijuanagebruik van jongens uit middelbare scholen in Californië. Men stelde de vraag: ’Wanneer gebruikte je voor het eerst marijuana?’ De event time is dan de leeftijd waarop voor het eerst marijuana gebruikt werd. Als een jongen antwoordde ’al gebruikt, maar ik weet niet meer wanneer’, ligt zijn exacte event time voor zijn huidige leeftijd. Laat deze jongen bijvoorbeeld 15

13 jaar zijn, dan zal dit resulteren in een linksgecensureerde observatie 13+.

leeftijd 10 11 12 13 14 15 16 17 18 >18

# exacte obs 4 12 19 24 20 13 3 1 0 4

# nog nooit gebruikt 0 0 2 15 24 18 14 6 0 0

# al eerder gebruikt 0 0 0 1 2 3 2 3 1 0

Tabel 2.2: Marijuanagebruik bij schooljongens Wanneer men in een studie rechts- en linksgecensureerde gegevens toelaat, noemt men de bijbehorende event times ook wel dubbel gecensureerd. We kunnen de observaties weer voorstellen door (T, ), maar nu zal T = max[min(X, Crechts ), Clinks ] 1 en 8 > <1 als niet gecensureerd = 0 als rechtsgecensureerd > : 1 als linksgecensureerd Voorbeeld 2.2.3 Wanneer in vorig voorbeeld er ook een jongen antwoordde: ’nog nooit gebruikt’, zal de event time rechtsgecensureerd zijn. In tabel 2.2 zien we dat er 2 twaalfjarige jongens zijn die een rechtsgecensureerde observatie geven. In dit onderzoek komen zowel links- als rechtsgecensureerde observaties voor, zodat deze dubbel gecensureerd is.

1

merk op dat T = min[max(X, Clinks ), Crechts ] hetzelfde resultaat geeft

16

Hoofdstuk 3

Goodness-of-fit tests Wanneer we willen nagaan of een dataset uit een bepaalde verdeling komt, hebben we verschillende verdelingstoetsen en grafische methodes tot onze beschikking. Denk maar aan de Chi-kwadraat verdelingstoets, de QQ-plot, EDF-statistieken,... Maar in het geval van gecensureerde gegevens zullen de klassieke toetsen een verkeerd resultaat geven. Dit is omdat de geobserveerde gegevens niet aan de klassieke verdeling voldoen, enkel de onderliggende event times. In dit hoofdstuk zullen we een aantal methodes aanpassen voor gecensureerde gegevens met een vaste censureringstijd. We behandelen het geval van vaste censureringstijd omdat we hier het voordeel hebben dat we de onderliggende verdeling makkelijk kunnen herkennen in de geobserveerde censureerde verdeling. Dit kunnen we zien in figuur 3.1: de rechtse grafiek is de verdelingsfunctie van een normale verdeling met gemiddelde 2 en variantie 9, de linkse grafiek is de empirische verdelingsfunctie van een steekproef van grootte 100 uit deze verdeling, maar dan (kunstmatig) gecensureerd op vaste censureringstijd gelijk aan 1.

Figuur 3.1: De verdelingsfunctie linksgecensureerde en de verdelingsfunctie van dezelfde niet-gecensureerde gegvens 17

3.1

Likelihoodconstructie

Eerst besteden we aandacht aan de constructie van de likelihoodfunctie want deze speelt een belangrijke rol bij het schatten van parameters. Herinner dat we de event times en censuringstijden onafhankelijk veronderstellen. Als dit niet het geval is, gaat deze methode niet meer op. We gaan na welke informatie iedere observatie ons geeft. Een observatie die hoort bij een exacte event time geeft infomatie over de kans dat de gebeurtenis plaatsheeft op deze tijd, wat ongeveer overeenstemt met de dichtheid van X op deze tijd. Voor een rechtsgecensureerde observatie weten we enkel dat de event time groter is dan deze tijd, wat overeenkomt met S(Crechts ). Wanneer we denken in termen van een kankerpatiënt dan weten we dat hij nog in leven is na tijd Crechts (met Crechts de random censureringstijd geassocieerd aan X). Analoog weten we voor een linksgecensureerde observatie dat de gebeurtenis al heeft plaats gehad voor deze tijd, zodat de informatiebijdrage gelijk is aan 1 S(Clinks ), of ook F (Clinks ). Wanneer we al deze informatie samenvoegen, kunnen we de likelihoodfunctie construeren: L/

X i2E

f (xi )

X

S(Crechts )

i2R

X

(1

S(Clinks ))

i2L

waarbij xi een exacte event time is wanneer i 2 E, een rechtsgecensureerde observatie wanneer i 2 R, en linksgecensureerd wanneer i 2 L. Dit kunnen we verder specifiëren tot een experiment waarin geen rechtse censurering voorkomt (we kunnen dit ook voor geen linkse censurering, de werkwijze is analoog). De geobserveerde data worden dan voorgesteld door paren (T, ) met = 1 wanneer de exacte event time geobserveerd wordt, en = 0 wannneer ze linksgecensureerd is. We kunnen dan deze twee gevallen bekijken: Voor

= 0 (de linksgecensureerde gegevens): P (T, = 0) = P (T = Clinks | = 0) · P ( = 0) = 1 · P ( = 0) = P (X  Clinks ) = 1 S(Clinks ) = F (Clinks )

Voor

= 1 (de niet gecensureerde gegevens): P (T, = 1) = P (T = X | = 1) · P ( = 1) = P  (X = T | X > Clinks ) · P (X > Clinks ) f (t) = · [S(Clinks )] = f (t) S(Clinks )

Dit kunnen we samenvoegen tot: P (T, ) = f (t) F (t)1

18

.

Voor een steekproef van grootte n met (Ti , i ), i = 1, .., n wordt dan de likelihood n n Y Y L= P (ti , i ) = f (ti ) i F (ti )1 i (3.1) i=1

i=1

of ook, wegens f (t) = h(t)·S(t) en formule 1.1 die de survivalfunctie uitdrukt in functie van de cumulatieve hazard rate functie. L=

n Y

h(ti ) i e

H(ti )

i

(1

e

)

H(ti ) 1

i

(3.2)

i=1

Deze likelihood zullen we later gebruiken om de parameters van de onderliggende verdeling te schatten door deze te maximaliseren.

3.2 3.2.1

Chi-kwadraat test Chi-kwadraat test voor niet gecensureerde gegevens

Beschouw steekproef X1 , .., Xn uit een populatie X, deze heeft een onbekende verdelingsfunctie F . Door middel van verdelingstoetsen kunnen we testen of deze verdelingsfunctie gelijk is aan een bepaalde continue theoretische verdeling. Hier onderzoeken we steeds of de steekproef uit een normale verdeling komt. In dit deel bespreken we de Chi-kwadraat verdelingstest. Deze is gebaseerd op het verschil tussen de geobserveerde aantallen en de verwachte (theoretische) aantallen. We doen een eenzijdige test: wanneer dit verschil groot is zullen we de nulhypothese verwerpen. De nulhypothese H0 is dan: steekproef X1 , .., Xn komt uit een normale verdeling met parameters µ ˆ en ˆ , waarbij deze laatste schattingen zijn aan de hand van de data. Hiervoor worden het steekproefgemiddelde en steekproefvariantie het meest gebruikt. We toetsen de nulhypothese door de theoretische verdeling op te splitsen in k klassen. Voor iedere klasse bepalen we de geobserveerde aantallen Nj , berekenen we theoretische kans pj en verwachte aantallen n · pj . Hierbij zorgen we er wel voor dat de verwachte aantallen voor alle j in {1, .., k} groter of gelijk zijn aan 5. Dit is nodig opdat de teststatistiek de gewenste verdeling zou hebben. Voor meer details over de afleiding van deze test zie Braekers[1]. Onder H0 hebben we dan voor de teststatistiek dat deze Ci-kwadraat verdeeld is met k 1 v vrijheidsgraden, waarbij k het aantal klassen en v het aantal geschatte parameters van de onderliggende verdeling is: TS =

k X (Nj j=1

npj )2 ⇡ npj

19

2

(k

1

v)

(3.3)

Dit gebruiken we als basis voor de test. Wanneer het verschil niet significant is op niveau ↵, zal de teststatistiek kleiner zijn dan kritisch punt (of ook wel kwantiel genoemd) 2(k 1 v;1 ↵) , in het andere geval zullen we H0 verwerpen. Dit kunnen we ook in p-waarden uitdrukken. Noteer p = P (T S ts), de kans onder nulhypothese dat de toetsingsgrootheid een waarde aanneemt die even extreem of extremer is dan de geobserveerde waarde. Wanneer p  ↵ zullen we de nulhypothese verwerpen.

Figuur 3.2: Opdeling in klassen voor Chi-kwadraat verdelingstest Bij de constructie van het algoritme hebben we gekozen (zie B.4.1) om eerst het aantal klassen en de theoretische kansen vast te leggen, en daaruit de begin- en eindpunten van de klasse-intervallen te berekenen. Omdat we ook de voorkeur geven om met kleine datasets te werken, kozen we voor 10 klassen met elk een theoretische kans pj gelijk aan 0,1 (1  j  10). Met deze waarden zien we dat de verwachte aantallen groter zullen zijn dan 5 bij steekproeven vanaf grootte 50. Zij n de steekproefgrootte, dan: n · pj

5 , n · 0, 1

5 , n

50

We willen testen of de steekproef uit een normale verdeling komt dus schatten we eerst het gemiddelde en de variantie door het steekproefgemiddelde en -variantie. n 1X x ¯ = xi n i=1 n 1 X s2 = (xi x)2 n 1 i=1

20

Nu kunnen we aan de hand van de theoretische kansen de begin- en eindpunten van de klasse-intervallen a1 , .., a9 berekenen. Want we weten dat onder de nulhypothese: P (X  a1 ) = P (a1 < X  a2 ) = .. .

(y1 ) = 0, 1 (y2 ) (y1 ) = 0, 1

P (a7 < X  a8 ) = (y8 ) (y7 ) = 0, 1 P (a8 < X  a9 ) = (y9 ) (y8 ) = 0, 1 P (a9 < X) = 1 (y9 ) = 0, 1 met yj =

x ¯

aj s

,1j9

(3.4)

Uit deze 10 vergelijkingen kunnen we de 9 onbekende waarden yj bepalen, 1  j  9 (d.m.v. tabellen voor de standaardnormale verdeling), waarna we waarden voor begin- en eindpunten van de klasse-intervallen aj (1  j  9) hieruit kunnen berekenen door lineaire vergelijkingen 3.4 op te lossen. De geobseerveerde waarden nj (1  j  10) bepalen we eenvoudigweg door het aantal in iedere klasse te tellen. Dit zien we ook terug in algoritme B.4.1. Nu kunnen we de waarde voor de teststatistiek ts (formule 3.3) berekenen en vergelijken met kritisch punt 2(7;1 ↵) . In het algoritme wordt gebruik gemaakt van de p-waarde p = P (T S ts) waarbij TS Chi-kwadraat verdeeld is met 7 vrijheidsgraden. Als deze groter is dan ↵ verwerpen we de nulhypothese onder significantieniveau ↵, anders nemen we aan dat de steekproef normaal verdeeld is met gemiddelde x ¯ en 2 variantie s onder niveau ↵. Voorbeeld 3.2.1 We passen het algoritme toe op een steekproef die we getrokken hebben uit een normale verdeling met gemiddelde 2 en standaarddeviatie 3, grootte 100 (voor de precieze waarden zie B.10). Wanneer we het algoritme toepassen krijgen we als output: > chitest10k(data,0.05) [1] "teststatistiek is" [1] 5 [1] "cutoff" [1] 14.06714 [1] "pwaarde is" [1] 0.6599632 [1] "de data komt uit een normale verdeling met mu en sigma gelijk aan:" [1] 1.913235 [1] 3.214045

21

De p-waarde is kleiner dan het kritisch punt wat betekent dat we op het 5% significantieniveau niet kunnen verwerpen dat de dataset uit een normale verdeling komt met gemiddelde gelijk aan 1.913235 en standaarddeviatie gelijk aan 3.214045. Voorbeeld 3.2.2 Als we het algoritme toepassen op een dataset van grootte 100 die getrokken is uit een weibull verdeling met parameter ↵ = 0, 5 en  = 3 krijgen we volgende output: > chitest10k(dataweib,0.05) [1] "teststatistiek is" [1] 141.8 [1] "cutoff" [1] 14.06714 [1] "pwaarde is" [1] 0 [1] "We werwerpen de nulhypothese op significantieniveau alpha" [1] "de data komt niet uit een normale verdeling met mu en sigma gelijk aan:" [1] 3.248130 [1] 5.173186 [1] 0

Zoals we verwachtten, wordt de nulhypothese verworpen op het 5%-niveau. In wat volgt simuleren we het significantieniveau. Met het significantieniveau houden we de type 1-fout onder controle. Een type 1-fout maken betekent dat we H0 verwerpen wanneer H0 juist is. We spreken van een toets op significantieniveau ↵ wanneer de kans op deze fout kleiner of gelijk is aan ↵. We kijken nu in hoeveel percent van de gevallen een goede steekproef verworpen wordt. teller2<-vector(length=20) for(j in 1:20) { teller<-0 for(i in 1:500) { data<-rnorm(500,2,3) p<-chitest10k(data,0.05) if(p<=0.05){teller<-teller+1} } teller2[j]<-teller/500 } teller2 #20 sec

We hebben 20 maal 500 datasets uit een normale verdeling met gemiddelde 2 en standaardeviatie 3 genomen en hierop het algoritme toegepast. Iedere keer noteren we het percent van de gevallen die verworpen worden. We zien 22

dat deze waarden dicht tegen het significantieniveau 0,05 liggen en dikwijls lager zijn. Merk op dat we hiervoor ook omgekeerd te werk kunnen gaan door een een dataset de nemen niet uit een normale verdeling komt. >teller2 # 10 klassen [1] 0.056 0.056 0.054 0.052 0.042 0.048 0.066 0.042 0.096 0.058 0.046 0.056 [13] 0.046 0.048 0.068 0.048 0.052 0.050 0.042 0.034

Wanneer we de steekproeven kleiner nemen, zullen deze percentages een beetje hoger liggen omdat de schattingen voor de parameters minder goed zijn en we testen t.o.v. een verdeling met deze schattingen als parameters en niet de exacte. Ook het aantal klassen speelt een rol. Wanneer we het algoritme herschrijven naar 5 klassen, is dit wat te weinig voor een steekproef van grootte 500. We zien dat de percentages dan hoger zullen liggen dan bij de test met 10 klassen, zodat deze laatste een beter resultaat geeft. >teller2 # 5 klassen [1] 0.070 0.086 0.080 0.070 0.080 0.088 0.072 0.080 0.066 0.096 0.096 0.072 [13] 0.080 0.080 0.072 0.076 0.100 0.076 0.088 0.054

3.2.2

Chi-kwadraat test voor linksgecensureerde gegevens

We beschouwen een steekproef T1 , .., Tn van linksgecensureerde gegevens met een vaste censureringstijd c, deze is afkomstig van een onderliggende steekproef X1 , .., Xn met verdeling F zoals in sectie 2.2. We kunnen nu niet gewoon van de geobserveerde waarden testen of ze uit een normale verdeling komen. Alle lifetimes die we niet observeren omdat ze in werkelijkheid kleiner zijn dan de censureringstijd, zullen verzameld worden in een piek op deze tijd. Onder de nulhypothese zal het de onderliggende verdeling zijn die normaal verdeeld is, en niet de geobserveerde. De klassen zullen daarom niet meer gelijkmatig verdeeld worden zoals bij de Chi-kwadraat verdelingstest voor ongecensureerde gegevens. We groeperen de linksgecensureerde gegevens in één klasse zodat het eerste interval ( 1, 1] wordt, de andere verdelen we gelijkmatig. We moeten wel opletten dat de verwachte aantallen van iedere klasse niet kleiner worden dan 5. In dat geval kunnen we klassen samennemen zodat kans pj om in deze nieuwe klassen terecht te komen zal verhogen en daarmee ook de verwachte aantallen npj . Ook zijn er andere schattingen nodig voor het gemiddelde en de variantie dan bij de Chi-kwadraat verdelingstest voor ongecensureerde gegevens. Want we willen niet de parameters van de geobserveerde verdeling bepalen maar wel die van de onderliggende. Hiervoor gebruiken we de likelihoodfunctie uit sectie 3.1. In formule 3.1 vullen we dan de normale verdelingsfunctie en dichtheidsfunctie in.

23

n  Y

L =

i=1

p

1 e 2⇡

1 xi ( 2

i

µ 2 )

·



(

xi

µ

1

i

)

Om schattingen te vinden voor de parameters maximalizen we de loglikelihoodfunctie naar parameters µ en . n ✓ X

ln L =

i=1

=

n ✓ h X

=

i

i=1 n ✓ X i=1

i

ln(e



 1 p ln e i 2⇡

1 xi ( 2

1 xi ( 2

µ 2 )

+ (1

)

p

ln( 2⇡)

ln

)2

1 ln(2⇡) 2

ln

µ 2 )

µ

1 xi ( 2

i

i ) ln



(

xi

+ (1

i ) ln

(

+ (1

i ) ln (

µ

)

c

µ

c

µ

◆ ◆ )

◆ ) (3.5)

Deze functie maximaliseren naar µ en is analytisch niet mogelijk wegens de normale verdelingsfunctie in de laatste term. Nummeriek kan dit wel. Dit doen we door op de functie uit B.3.1 in het programma R een niet-lineaire minimalizatie toe te passen met als startwaarden het steekproefgemiddelde en steekproefvariantie. We kiezen deze startwaarden omdat deze in de buurt zullen liggen van de parameters, om zo het aantal iteraties in de minimalizatie te beperken. We testen dit eens uit op kunstmatig linksgecensureerde steekproef. Voorbeeld 3.2.3 Neem dezelfde ongecensureerde steekproef als in vb.3.2.1. Hieruit construeren we een linksgecensureerde dataset d.m.v. algoritme B.1.1. We minimaliseren ln L met: >nlm(loglikfct(cendata,delta(cendata,1),1),stwa,print.level=1) $estimate [1] 1.790869 3.359120 In figuur 3.3 wordt de contourplot van de loglikelihoodfunctie getoond. Het rode punt staat voor de parameterschatting, het snijpunt van de rechten voor de werkelijke waarde. Deze plot werd gemaakt met het algoritme uit B.3.2. We zien dat de schatting in de buurt van de werkelijke waarde ligt.

24

Figuur 3.3: Contourplot van de loglikelihood in vb.3.2.3 Bij de constructie van het algoritme kozen we voor een klasse met linksgecensureerde gegevens en zes andere klassen. Een schematische voorstelling hiervan zien we in figuur 3.4. Deze maakt gebruik van vaste censureringstijd gelijk aan 1.

Figuur 3.4: Opdeling in klassen voor Chi-kwadraat verdelingstest bij linksgecensureerde gegevens Aan de hand van de theoretische kansen kunnen we dan de begin- en eind25

punten van de klasse-intervallen berekenen. We weten: P (X  c) = P (c < X  a1 ) =

P (a1 P (a2 P (a3 P (a4

< X  a2 ) < X  a3 ) < X  a4 ) < X  a5 ) P (a5 < X)

= = = = =

(

µ ˆ

c

) ⌘ Pcens c µ ˆ ( ) = (1 Pcens )/6 ˆ (y2 ) (y1 ) = (1 Pcens )/6 (y3 ) (y2 ) = (1 Pcens )/6 (y4 ) (y3 ) = (1 Pcens )/6 (y5 ) (y4 ) = (1 Pcens )/6 1 (y5 ) = (1 Pcens )/6 ˆ (y1 )

met yj =

µ ˆ

aj ˆ

,1j5

(3.6)

Vaste censureringstijd c is gekend en schattingen µ ˆ en ˆ zijn eerder berekend. Hierdoor kunnen we Pcens berekenen. De overige klassen krijgen een theoretische kans die gelijk is aan de overblijvende kans gedeeld door het aantal overige klassen. Dan kunnen de andere punten berekend worden zoals eerder. Wat verder volgt in het algoritme is analoog aan het algoritme voor niet gecensureerde gegevens. Alleen nog even opletten bij het aantal vrijheidsgraden van de chi-kwadraatverdeling, dat wordt nu k v 1 = 7 1 2 = 4. Er werd ook een algoritme met 11 klassen geschreven, zodat ook voor grotere steekproeven een goed resultaat bereikt kan worden. Deze algoritmes vinden we in B.4.2. We kunnen weer het significantieniveau simuleren. We gebruiken hiervoor 20 keer 500 datasets van grootte 500 met een onderliggende normale verdeling met gemiddelde 2 en standaardeviatie 3 en vaste censureringstijd gelijk aan 1. We noteren het percentage steekproeven waarvan op niveau 0,05 de nulhypothese van normaliteit verworpen wordt. De resultaten zijn dan: # voor 11k n=500 10 min >teller2 [1] 0.032 0.036 0.044 0.044 0.048 0.048 0.048 0.048 0.050 0.050 [11] 0.050 0.052 0.052 0.054 0.054 0.058 0.058 0.060 0.062 0.068 # voor 7k n=500 6 min20 >teller2 [1] 0.038 0.040 0.042 0.042 0.046 0.048 0.052 0.056 0.058 0.058 [11] 0.060 0.062 0.062 0.062 0.062 0.062 0.066 0.068 0.070 0.086

Hieruit kunnen we aflezen dat het algoritme voor linksgecensureerde gegevens heel wat langer duurt dan hetgeen voor ongecensureerde. Dit is te

26

wijten aan de kostelijke minimalisatieprocedure om de schatters te bepalen in het geval van gecensureerde gegevens. We zien, net zoals bij ongecensureerde gegevens, dat voor grote steekproeven het algoritme met meer klassen kleinere percentages opleverd die dichter bij het gewenste significantieniveau 0,05 liggen.

3.3

QQ-plot

De QQ-plot is een verdelingstest waarbij men visueel kan nagaan of een steekproef X1 , .., Xn overeenstemt met een theoretische verdeling F . Deze is gebaseerd op het feit dat het menselijk oog het verschil kan nagaan tussen een rechte lijn en een curve. De QQ-plot is dan gedefinieerd als de grafiek van de verwachte waarden onder de nulhypothese (wanneer de steekproef uit verdeling F komt) t.o.v. de geobserveerde waarden van de ordestatistieken. ⇥ ⇤ (E X(i) , x(i) ) | 1  i  n . Wanneer de punten dicht bij de bissectrice (op 45 graden) liggen, zullen de verwachte waarden ruwweg overeenstemmen met de feitelijke waarnemingen. In het geval van linksgecensureerde gegevens met vaste censureringstijd willen we testen of de verdeling van de onderliggende steekproef X1 , .., Xn overeenkomt met een verdeling F (zie formule 2.1). We berekenen de verwachte waarden onder de nulhypothese door te vertrekken van deze verdeling en hieruit de theoretische verdeling van de exacte (geobserveerde) steekproef T1 , .., Tn te berekenen. Van deze berekenen we de ordestatistiek en zodoende de verwachtingswaarde. We specifiëren verdeling F hier weer door een normale verdeling. Als we een andere verdeling kiezen voor Xi verkrijgen we bij het opstellen van de QQ-plot op analoge manier andere verwachte waarden voor deze verdeling. Verder specifiëren we hier de censureringstijd als 1. Voor een andere waarde zijn de berekeningen analoog. We bevinden ons dus in de situatie dat onder de nulhypothese de onderliggende stochastische variabele X verdeeld is volgens N (µ, 2 ). De geobserveerde variabele is dan T = max(X, 1). Hiervan berekenen we de verdelingsfunctie: F T (y) = = = =

P (T  y) P (max(X, 1)  y) P (X  y en 1  y) P (X  y) · 1[1,1[ 27

( F X (y) als y = 0 als y < 1

1

met F X de normale verdelingsfunctie. We weten uit de theorie over ordestatistieken dat de verdeling van de k-e ordestatistiek gelijk zal zijn aan volgende uitdrukking T F(k) (y)

=

n ✓ ◆ X n j=k

j

F T (y)j (1

F T (y))n

j

.

(3.7)

Event times zijn steeds positieve stochastische variabelen. Wegens appendix B9 is voor een positieve stochastische variabele Y , met verdeling F en dichtheid f de verwachtingswaarde gedefinieerd als: Z 1 Z 1 E [Y ] = yf (y)dy = [1 F (y)] dy. 1

0

Zodat de verwachtingswaarde voor T(k) wordt: ⇥ ⇤ E T(k) =

Z

0

1h

i T F(k) (x) dx

1

(3.8)

We herschrijven het argument van door middel van het biPndezenintegraal n i n i nomium van Newton: (a + b) = i=1 i x y : 1 = (F T (x) + (1 F T (x)))n n ✓ ◆ X n = F T (x)i (1 F T (x))n i

i

i=1

We zien dat de uitdrukking uit 3.7 en vorige formule sterk op elkaar gelijken. Het argument wordt dan:

1

T F(k)

=

n ✓ ◆ X n i=1

=

i

k 1✓ ◆ X n j=1

j

F (x) (1 T

i

F (x)) T

n i

n ✓ ◆ X n j=k

F T (x)j (1

F T (x))n

j

F T (x)j (1

j

Hierin herkennen we een binomiaalverdeling. Stel Y binomiaalverdeeld met parameters n en p, dan geldt ✓ ◆ n P (Y = y) = py · (1 p)(n y) y 28

F T (x))n

j

zodat P (Y  y) = Als we p = F T (x) en y = k 1

y ✓ ◆ X n j=0

j

pj · (1

p)(n

j)

1 nemen dan is T F(k) = P (Y  k

1)

Zodat wanneer we dit invullen in 3.8 we een uitdrukking krijgen voor de gewenste verwachtingswaarde. Voor 1  k  n: ⇥ ⇤ E T(k) =

Z

0

1

P (Y  k

1)dx met Y ⇡ B(n, F T (x))

Nu kunnen we de QQ-plot opstellen. In het algoritme (zie B.5) maken we gebruik van voorgaande berekeningen om twee vectoren van lengte n op te stellen. Op de horizontale as zetten we de verwachte waarden, op de vertikale de observaties van groot naar klein. We gaan visueel na of de punten zich in de buurt van de bissectrice bevinden. Maar ’dichtbij’ is nogal subjectief, zodat we best nog een bijkomende andere verdelingstest uitvoeren bij twijfel.

Figuur 3.5: QQ-plot van steekproef uit normale verdeling N (2, 32 ) t.o.v. normale verdelingen Voorbeeld 3.3.1 We beschouwen weer dezelfde normaal verdeelde dataset als eerder (met parameters 2 en 3). We kunnen nu kiezen hoe we de parameters van de theoretische verdeling specifiëren. In de linkse grafiek van figuur 3.5 kozen we de exacte parameters, in de rechtse schatten we ze met de likelihoodmethode. 29

We zien erg weinig verschil tussen beide grafieken. Zo goed als alle punten liggen voor beide erg dicht bij de bissectrice, enkel de punten in de staart liggen wat verder. We moeten nu echter wel opletten tegen welke normaalverdeling we testen. In het geval van ongecensureerde gegevens was het niet nodig om µ en te kennen wegens: Zij Y1 , .., Yn een steekproef van grootte n. Als Yi ⇡ N (µ, 2 ) ) Zi = Yi µ ⇡ N (0, 1) ) QQ-plot met theoretische verdeling N (0, 1) van y(i) lineair ) QQ-plot met theoretische verdeling N (0, 1) van z(i) = y(i) + µ lineair.

Maar nu geldt dit niet meer aangezien het niet meer de geobserveerde stochastische variabele is die normaal verdeeld is, maar de onderliggende. Voorbeeld 3.3.2 We testen de normale dataset met parameters µ = 2 en = 3 tegen een standaardnormale verdeling. In figuur 3.6 zien we duidelijk dat de punten de bissectrice helemaal niet benaderen.

Figuur 3.6: QQ-plot van steekproef uit normale verdeling N (2, 32 ) t.o.v. standaardnormale verdeling Tenslotte bekijken we nog een gecensureerde steekproef uit een Weibullverdeling. Voorbeeld 3.3.3 We beschouwen dezelfde dataset als in vb. 3.2.2. Hieruit constueren we een linksgecensureerde dataset met vaste censureringstijd gelijk aan 1. De 30

parameters voor de theoretische normale verdeling schatten we met de likelihoodmethode. Figuur 3.7 geeft dan de QQ-plot, die sterk verschilt van de bissectrice.

Figuur 3.7: QQ-plot N(0.9969194,69.01996)

3.4

van

steekproef

uit

Weib(1/2,3)

tegen

EDF-statistieken

In dit deel gaan we wat dieper in op goodness-of-fit testen die gebaseerd zijn op de empirische verdelingsfunctie, ook wel empirische distributiefunctie (EDF) genoemd. Deze trapfunctie is een consistente schatter voor de populatieverdelingsfunctie. EDF statistieken zijn een maat voor het verschil tussen de EDF en een gegeven theoretische verdelingsfunctie. We maken een onderscheid tussen supremumstatistieken en kwadratische statistieken. Beide zijn gebaseerd op vertikale afstanden tussen de EDF Fn en de verdelingsfunctie F , maar verschillen in de gebruikte norm.

3.4.1

Empirische distributiefunctie

Onderstel dat er een steekproef van grootte n X1 , ..Xn gegeven is, dan is de EDF gedefinieerd als: aantal observaties  x voor n 8 > <0 als x < X(1) Fn (x) = i als X(i)  x  X(i+1) > : 1 als X(n)  x

Fn (x) =

31

1<x<1

Voor een x kleiner dan de kleinste observatie zal deze trapfunctie nul zijn, voor een x groter dan de grootste observatie één. Ertussen stijgt ze steeds met 1/n wanneer index i stijgt, met n de steekproefgrootte. Als interpretatie heeft Fn (x) de proportie van alle observaties kleiner dan x en F (x) de kans dat een observatie kleiner is dan x. Wanneer n ! 1 zal | Fn (x) F (x) | naar nul gaan met kans één, zodat Fn een consistente schatter is voor F . Voorbeeld 3.4.1 In figuur 3.8 zien we de EDF van de Leghorn Chicken dataset (zie B.10). Dit is een steekproef bestaande uit het gewicht (in gram) van 20 kippen van 21 dagen oud.

Figuur 3.8: EDF van Leghorn Chicken dataset

3.4.2

Supremumstatistieken

Hieronder verstaan we de D+ , D en de Kolmogorov-statistiek D. Statistiek D+ (resp. D ) is de grootste vertikale afstand tussen de EDF en de verdelingsfunctie wanneer Fn (x) groter is dan F (resp. kleiner). Formeel definiëren we deze als D+

=

sup(Fn (x)

F (x))

D

=

sup(F (x)

Fn (x))

x x

32

Maar aangezien deze soms niet zo’n goed resultaat resultaat geven, beschouwen we ook nog de meest bekende supremumstatistiek D, de Kolmogorov statistiek. D = sup | Fn (x)

F (x) |= max(D+ , D ) (Kolmogorov)

x

Deze statistiek is gebaseerd op de L1 norm. In figuur 3.9 zien we de EDF van de Leghorn Chicken dataset Fn en de normale verdelingsfunctie met gemiddelde 200 en standaarddeviatie 35. Statistieken D+ en D zijn aangeduid. We zien dat de Kolmogorov statistiek D gelijk zal zijn aan D .

Figuur 3.9: Supremumstatistieken van Leghorn Chicken dataset

3.4.3

Kwadratische statistieken

Een tweede grote klasse van statistieken is die van de Cramér-von Mises familie. Deze is gebaseerd op de L2 -norm en is gedefinieerd als: Z 1 Q = n (Fn (x) F (x))2 · '(x) dF (x). 1

Functie '(x) geeft gewichten aan de kwadratische verschillen. We bekijken twee statistieken van deze vorm: • de Cramér-von Mises statistiek W 2 met '(x) = 1 • de Anderson-Darling statistiek A2 met '(x) = [F (x) · (1

33

F (x))]

1

3.4.4

Berekenen van de statistieken: PIT

In dit deel construeren we formules voor de statistieken uit sectie 3.4.3 die makkelijker te implementeren zijn, en dus meer bruikbaar voor het schrijven van algoritmes. Dit doet men door gebruik te maken van de Probability Integral Transformation (PIT): Z ⌘ F (X). In figuur 3.10 zien we hiervan een schematische voorstelling.

Figuur 3.10: Schematische voorstelling PIT Zij F de echte (continue) verdeling van X dan zal de stochastische variabele Z uniform verdeeld zijn tussen nul en één. Dus Z heeft als verdelingsfunctie F Z (z) = z , 0  z  1. Dit kunnen we inzien met volgende afleiding: F Z (z) = = = =

P (Z  z) P (F (X)  z) P (X  F 1 (z)) F (F 1 (z)) = z

als z 2 [0, 1] als z 2 [0, 1]

Beschouw steekproef X1 , .., Xn . Deze geeft waarden Zi = F (Xi ) voor 1  i  n. We noteren de empirische verdelingsfunctie van deze laatste steekproef met FnZ (z). Dan gelden volgende gelijkheden: Fn (x)

F (x) = Fnz (z)

F z (z) = Fnz (z)

z

want n

Fn (x)

F (x) =

1X 1(Xi n i=1

 x)

F (x)

34

def initie EDF X

(3.9)

n

= = =

1X 1(Xi  F 1 (z)) z n i=1 n 1X 1(F (Xi )  z) z n i=1 n 1X 1(Zi  z) z n

P IT en z = F F niet

1

(x)

dalend monotoon

P IT

i=1

= FnZ (z) = FnZ (z)

z F Z (z)

def initie EDF Z Z unif orm verdeeld.

Dankzij formule 3.9 geven EDF-statistieken voor steekproef Z1 , .., Zn hetzelfde resultaat als voor X1 , .., Xn . Dit leidt tot volgende belangrijke formules. D+ = =

D D

W

2

=

n X i=1

A2 =

n

max (

1in

i n

max (Zi

Zi ) 1

i

n = max(D+ , D ) 1in

)

(2i 1) 2 1 ) + 2n 12n n 1 X⇥ (2i 1) ln Z(i) + (2n + 1 n

(3.10)

(Z(i)

i=1

2i) ln(1

⇤ Z(i) ) (3.11)

De uitdrukkingen voor D+ en D volgen via de definitie van EDF (dit is makkelijk te zien op figuur 3.9) en het feit dat het maximum steeds bestaat omdat de steekproef eindig is uit formule 3.9. De formule voor W 2 kunnen we als volgt afleiden:

W

2

= n

Z Z

+1

[Fn (x)

F (x)]2 dF (x)

1 1⇥

⇤2 FnZ (z) F Z (z) dz 0 #2 Z 1" X n 1 = n ( 1(Z(i)z ) ) z dz n 0 i=1 3 # 2 Z 1" X n n X 1 1 = n ( 1(Z(i)z ) ) z · 4( 1(Z(j)z ) ) z 5 dz n n 0 j=1 2 i=1 3 Z 1 n X n n X X 1 2 4 = n 1(Z(i)z ) 1(Z(j)z ) 1(Z(i)z ) · z + z 2 5 dz 2 n n 0 = n

i=1 j=1

i=1

35

n

=

n

1 XX n i=1 j=1

= = =

Z

0

1h

1 n

i=1 j=1

z |1max{Z

n n 1 XX⇥ 1 n

2

1(max{Z(i) ,Z(j) }z) dz

n n Z 1 XX 1 dz n max{Z(i) ,Z(j) } i=1 j=1 n X n X

i

(i) ,Z(j) }

2

⇤

i=1 j=1

0

i=1

n Z X

1

n 3 n 1 Xh

2·

1h

2

i

1(Z(i)z ) · z dz + n

z dz + n ·

j=1 Z(i) n X z2 2 |1 2 Z(i) i=1

max Z(i) , Z(j)

n Z X

z3 1 | 3 0

+

i n 2 Z(i) + 3

1

i=1

Omdat Z(i) , Z(j) ordestatistieken zijn, geldt max Z(i) , Z(j) = Z(j) als i  j en Z(i) anders. We splitsen de som over j op in een som van 1 tot i en een som van i + 1 tot n zodat in de eerste i steeds groter is dan j, en in de tweede j groter dan i. 2 3 n i n X 1 X 4X Z(i) + Z(j) 5 n

= n

i=1

j=1

n+

j=i+1

n X

2 Z(i) +

i=1

n 3

Als we nu de eerste som over j uitwerken, en bij de tweede de sommen omwisselen krijgen we: n j 1

n

1X i · Z(i) n

=

i=1

1 Z n (1)

= =

j=2 i=1

n

1X iZ(i) n i=2

n 1X (2i n

n

X 1 XX n 2 Z(j) + Z(i) + n 3

1)Z(i) +

i=1

i=1

n

1X (j n

1)Z(j) +

j=2

n X

2 Z(i) +

i=1

n X

2 Z(i) +

i=1

n 3

n 3

En dit is gelijk aan: =

n X i=1

n

2 Z(i)

1X (2i n i=1

n X 2i 1 2 1 1)Z(i) + ( ) + 2n 12n

(3.12)

i

Want na het uitwerken van het kwadraat hebben we met volgende gelijkheden: n X

i2 =

i=1 n X i=1

i =

n(n + 1)(2n + 1) 6 n(n + 1) 2 36

Z

0

1

z 2 dz

dat n n X X 2i 1 2 1 1 ( ) + = 2 (4 i2 + n 2n 12n 4n i=1

i=1

4

n X i=1

i) +

1 n = . 12n 3

Als we in 3.12 de som naar voor halen, kunnen we deze uitdrukking via het merkwaardig product (a + b)2 = a2 + 2ab + b2 verder uitwerken tot 3.11: n  X 1 2i 1 2 1 2 2 W = Z(i) · Z(i) · (2i 1) + ( ) + n 2n 12n i=1 n  X 2i 1 2i 1 2 1 2 = Z(i) 2 · Z(i) · +( ) + 2n 2n 12n i=1  n 2 X 2i 1 1 = Z(i) + (3.13) 2n 12n i=1

Dit is de formule die we verder zullen gebruiken voor de Cramér-von Mises statistiek. Ook voor de Anderson-Darling statistiek kan men zo’n afleiding vinden, maar die wordt aan de ge¨ınteresseerde lezer overgelaten.

We zullen formules 3.10 en 3.11 in het volgende deel gebruiken om goodnessof-fit tests gebaseerd op EDF-statistieken te construeren. We testen weer de nulhypothese dat X1 , .., Xn een steekproef is uit een bepaalde theoretische verdeling. De algoritmes zijn gescheven voor een normale verdeling. Voor andere verdelingen kan men een analoge redenering volgen. We maken het onderscheid tussen een volledig gespecifieerde verdeling, wat we case 0 noemen, en de situatie dat één of meerdere parameters onbekend zijn.

3.4.5

EDF-tests: case 0

Voor niet gecensureerde gegevens We willen de nulhypothese testen dat steekproef X1 , .., Xn uit een welgespecifieerde normale verdeling N (µ, 2 ) komt. Hiervoor hebben we volgende werkwijze: Werkwijze case 0 1. Sorteer Xi : X(1) < .. < X(n) . 2. Bereken Z(i) = F (X(i) ; ✓), 1  i  n, met F de verdelingsfunctie van de theoretische normale verdeling. 3. Bereken de gewenste statistiek met voorgaande formules 3.10, 3.11. 37

4. Vergelijk de aangepaste statistiek (deze kunnen we berekenen met formule 3.14) met de juiste lijn van kritische punten in tabel 3.1. Als de statistiek de waarde van het kritisch punt voor een gegeven ↵ overschrijdt, dan verwerpen we H0 op significantieniveau ↵. In de eerste twee stappen berekenen we de PIT zodat we in de derde stap de gewenste statistiek kunnen berekenen. De kritische punten zijn de kwantielen G 1 (1 ↵) waarbij G de verdeling is van de toetsingsgrootheid onder de nulhypothese en ↵ het gewenste signigicantieniveau. Opdat de verdeling van de toetsingsgrootheid beter zou aansluiten bij de asymptotische verdeling werken we met aangepaste statistieken: p p (D+ )⇤ = D+ (pn + 0.12 + 0.11/pn) (D )⇤ = D p ( n + 0.12 + 0.11/ p n) ⇤ (D) = D( n + 0.12 + 0.11/ n) (W 2 )⇤ = (W 2 0.4/n + 0.6/n2 )(1.0 + 1.0/n) A⇤ = A voor alle n 5 (3.14) De kritische punten voor eindige steekproefgrootte werden gevonden door Monte Carlo studies en de aangepaste statistieken door onderzoek hoe deze punten naar het asymptotisch punt convergeerden (zie Stephens[3]). Zoals bij de Chi-kwadraat verdelingstest is de EDF-test een eenzijdige verdelingstest en verwerpen we de nulhypothese op niveau ↵ als de waarde van de teststatistiek groter is dan cut-o↵ waarde voor die ↵. T* (D± )⇤

(D)⇤ (W 2 )⇤ A⇤

signf. 0.25 0.828 1.019 0.209 1.248

lvl 0.15 0.973 1.138 0.284 1.610

↵ 0.10 1.073 1.224 0.347 1.933

0.05 1.224 1.358 0.461 2.492

0.025 1.358 1.480 0.581 3.070

0.01 1.518 1.628 0.743 3.857

0.005 1.628 1.731 0.869 4.500

0.001 1.859 1.950 1.167 6.00

Tabel 3.1: Aangepaste statistieken en kritische punten voor EDFstatistieken: ongecensureerd, case 0 Deze werkwijze hebben we in twee algoritmes gegoten, die we vinden in B.6.2. Het eerste algoritme berekent de gewone statistieken, het tweede de aangepaste. De output kunnen we dan vergelijken met de kritische punten uit de tabel. Dit passen we toe op enkele voorbeelden. Voorbeeld 3.4.2 Beschouw de Leghorn Chicken dataset (B.10). Hierop passen we achtereenvolgens de algoritmes voor het berekenen van de EDF-statistieken voor case 0 en voor de aangepaste EDF-statistieken (zie B.6.2) toe. Wanneer we µ = 200 en = 35 kiezen krijgen volgende output: 38

> ongecensEDF(X,200,35) Dplus Dmin D Wkwadr Akwadr [1,] 0.04437504 0.1712159 0.1712159 0.1874563 1.016849 > ongecensEDFster(X,200,35) Dplusster Dminster Dster Wkwadrster Akwadrster [1,] 0.2048677 0.7904581 0.7904581 0.1774041 1.016849

Als we deze waarden vergelijken met de kritische punten zien we dat, voor alle significantieniveaus, alle aangepaste statistieken kleiner zijn dan deze cuto↵-waarden. Op elk niveau aanvaarden we de nulhypothese voor iedere statistiek. We kunnnen niet verwerpen dat de steekproef normaal verdeeld is met parameters 200 en 35. Voorbeeld 3.4.3 We nemen weer de data uit N (2, 32 ) en berekenen de aangepaste EDFstatistieken. > ongecensEDFster(data,2,3) Wkwadrster Akwadrster [1,] 0.7623949 0.2606439 0.7623949 0.07946807

0.4813691

Voor niveau 0,25 liggen de waarden van alle statistieken goed onder de kritische punten. Enkel (D+ )⇤ is wat groter, maar nog niet zo groot dat we de nulhypothese verwerpen. Op figuur 3.11 zien we goed dat de EDF vooral positieve afwijkingen vertoond t.o.v. de verdelingsfunctie en bijna geen negatieve, wat de grote waarde voor (D+ )⇤ en de kleine voor (D )⇤ verklaart. Figuren van deze vorm werden getekend met de algoritmes uit B.6.1.

Figuur 3.11: EDF van de Xi en Zi uit vb. 3.4.3 t.o.v. hun verdelingen EDF-statistieken zijn meestal sterker dan de Chi-kwadraat statistiek. De reden hiervoor is dat we informatie verliezen bij het groeperen in klassen bij de Chi-kwadraat verdelingstest. Ook zijn (volgens Stephens[3]) de kwadra39

tische statistieken W 2 en A2 dikwijls sterker dan de Kolmogorov statistiek D. De Anderson-Darling statistiek A2 en de Cramér-von Mises statistiek W 2 gedragen zich dikwijls gelijkaardig, maar wanneer het verschil tussen de EDF en de theoretische verdeling zich manifesteert in de staarten (i.h.b. in het geval van uitschieters) zal A2 een beter resultaat geven. Supremumstatistieken D+ en D worden bijna nooit gebruikt wanneer men gewoon wil testen of de steekproef uit een bepaalde verdeling komt, omdat deze snel een verkeerd resultaat kunnen geven, zoals we in volgend voorbeeld kunnen zien. Voorbeeld 3.4.4 We hernemen vb. 3.4.3, maar nu testen tegen de normale verdeling met parameters µ = 6 en = 5: > ongecensEDFster(data2,6,5) Dplusster Dminster Dster Wkwadrster Akwadrster [1,] 6.37446 0.06812284 6.37446 16.84782 79.97709

We beschouwen significantieniveau 0,05. Hierop wordt de nulhypothese voor alle statistieken sterk verworpen, behalve voor (D )⇤ . We zien op figuur 3.12 de EDF bijna overal ver boven de normale verdelingsfunctie (resp. de uniforme) ligt zodat statistiek D enkel dicht bij 0 verschillend is van nul. En voor kleine waarden van de statistiek wordt de nulhypothese niet verworpen. Statistiek D+ geeft in dit geval wel het juiste resultaat.

Figuur 3.12: EDF van de Xi en Zi uit vb. 3.4.4 t.o.v. hun verdelingen

Voor rechtsgecensureerde gegevens We beschouwen een rechtsgecensureerde steekproef T1 , .., Tn , met vaste censureringstijd c. We noemen r het aantal ongecensureerde observaties. Dan 40

zijn er r observaties exact, wat betekent dat ze kleiner zijn dan c. De andere zijn rechtsgecensureerd en dus gelijk aan c We willen de nulhypothese testen dat de onderliggende verdeling van steekproef uit een welgespecifieerde normale verdeling N (µ, 2 ) komt. De theoretische verdeling F is volledig gespecifieerd zodat we de Probability Integral Transformation kunnen gebruiken om steekproef Zi , .., Zn te construeren. Deze steekproef zal ook gecensureerd zijn want Z(1)  Z(2)  ..  Z(r) < t met t = F (c) We passen eerst de Kolmogorov statistiek D aan voor rechtse censurering. We nemen het supremum nu niet meer over z 2 [0, 1] maar beperken tot [0, t]. In vb. 3.4.5 kunnen we intu¨ıtief zien waarom de statistiek gelijk is aan volgende uitdrukking. Dt = =

sup | FnZ (z)

0zt ⇢

max

1ir

i n

z|

Z(i) , Z(i)

1

i n

, t

r n

(3.15)

Voorbeeld 3.4.5 We kunnen de Leghorn Chicken dataset kunstmatig rechts censureren door alle observaties groter dan 235 gram gelijk te stellen aan de vaste censureringstijd 235. We bekomen dan volgende steekproef: [1] 156 162 168 182 186 190 190 196 202 210 [11] 214 220 226 230 230 235 235 235 235 235

deze zetten we met de PIT om in de Z steekproef. Waarde t wordt dan 0.8413447 en waarden Z(1) , .., Z(r) : [1] 0.1043510 0.1388027 0.1802834 0.3035261 0.3445783 0.3875485 0.3875485 [8] 0.4545057 0.5227843 0.6124515 0.6554217 0.7161454 0.7712159 0.8043170 [15] 0.8043170

Deze gebruiken we om de EDF voor de Z-steekproef op te stellen. Voor de ongecensureerde observaties (1  i  r) is er geen verschil zodat we de eerste twee argumenten kunnen overnemen van het ongecensureerde geval. Voor de gecensureerde observaties zal de uniforme verdeling steeds hoger liggen dan de EDF zodat we voor observaties met i > r argument t nr toevoegen. Merk op dat dit overeenkomt met de steekproef beperken tot enkel de ongecensureerde observaties, en daarna een observatie Z(r+1) ⌘ t toevoegen.

41

Figuur 3.13: De EDF van de Zi t.o.v. de uniforme verdeling voor vb. 3.4.5 Ook voor de Cramér-von Mises statistiek W 2 en de Anderson-Darling statistiek A2 bestaan er aanpassingen voor rechtse censurering (zie Stephens[3]). Hiervoor beschouwen we Z(1) , .., Z(r) en voegen we het punt Z(r+1) ⌘ t er aan toe. Vertrekkende van de definitie van deze kwadratische statistieken (zoals in 3.4.3), kan men op een analoge manier als de afleiding van formule 3.13 voor W 2 volgende uitdrukkingen bekomen. Wt2 =

r+1 X

(Z(i)

i=1

A2t

=

r+1

1X (2i n i=1

2i 1 2 r+1 n ) + + (Z(r+1) 2 2n 12n 3 ⇥ 1) ln Z(i)

1⇥ (r + 1 n

ln(1

n)2 ln(1

t)

⇤ Z(i) )

2

r+1 X

r+1 3 ) n ln(1

Z(i) )

i=1

(r + 1)2 ln t + n2 t

⇤

(3.16)

Net als bij ongecensureerde gegevens bestaan er voor de asymptotische verdeling van deze teststatistieken ook tabellen met kritische punten en een aanpassing van de statistieken voor eindige steekproeven. De werkwijze voor de EDF-test voor rechtsgecensureerde gegevens is analoog aan die voor ongecensureerde. Werkwijze case 0: rechtsgecensureerde gegevens 1. Bereken Z(i) = F (X(i) ; ✓), i = 1, .., r, met F de theoretische normale verdeling waartegen we willen testen. En stel Z(r+1) gelijk aan vaste censureringstijd t van de Z steekproef. 42

2. Bereken de gewenste statistiek met formule 3.15 voor Dt of 3.16 voor Wt2 of A2t . 3. Bereken in het geval van statistiek Dt de aangepaste statistiek Dt⇤ = p p nDt + 0.19/ n voor n 25 en t 0, 25. Voor Wt2 of A2t is er geen aanpassing nodig 4. Vergelijk de aangepaste statistiek met de juiste tabel van kritische punten. Als de statistiek de waarde van het kritisch punt voor een gegeven ↵ overschrijdt, dan verwerpen we H0 op significantieniveau ↵ De voorwaarde t 0, 25 betekent dat de steekproef minstens 25 percent ongecensureerde observaties moet bevatten opdat deze methode nog bruikbaar is. Voor extremere waarden en steekproeven kleiner dan 25 gegevens bestaan er ook oplossingen. Hiervoor en voor de tabellen voor Wt2 en A2t verwijzen we naar de litereatuur: Stephens[3]. D⇤ r/n .2 .3 .4 .5 .6 .7 .8 .9 1.0

signf. 0.50 .4923 .5889 .6627 .7204 .7649 .7975 .8183 .8270 .8273

lvl 0.25 .6465 .7663 .8544 .9196 .9666 .9976 1.0142 1.0190 1.0192

↵ 0.15 .7443 .8784 .9746 1.0438 1.0914 1.1208 1.1348 1.1379 1.1379

0.10 .8155 .9597 1.0616 1.1334 1.1813 1.2094 1.2216 1.2238 1.2238

0.05 .9268 1.0868 1.1975 1.2731 1.3211 1.3471 1.3568 1.3581 1.3581

Tabel 3.2: Kritische punten voor

p

0.025 1.0282 1.2024 1.3209 1.3997 1.4476 1.4717 1.4794 1.4802 1.4802

0.01 1.1505 1.3419 1.4696 1.5520 1.5996 1.6214 1.6272 1.6276 1.6276

0.005 1.2361 1.4394 1.5735 1.6583 1.7056 1.7258 1.7306 1.7308 1.7308

p nD + 0.19/ n

De kritische punten verschillen voor verschillende percentages van ongecensureerde observaties nr . In de laatste rij van tabel 3.2 herkennen we de kritische punten voor ongecensureerde steekproeven. Voor een percentage r n dat niet voorkomt in de tabel kunnen we lineaire interpolatie gebruiken. Het algoritme om deze statistieken te berekenen vinden we in B.6.3. Dit gebruiken we in volgend voorbeeld. Voorbeeld 3.4.6 We beschouwen de normaal verdeelde steekproef uit B.10. Hiervan maken we een rechtsgecensureerde dataset met vaste censureringstijd gelijk aan 4. Hierop passen we het algoritme toe. 43

> rechtsgecensEDF(cenrechts(data,4),deltarechts(data,4),4,2,3) D Dster Wkwadr Akwadr pct ongecensureerd [1,] 0.07525366 0.7715366 0.07892602 0.8077139 0.76

Wanneer we a.d.h.v. statistiek (Dt )⇤ op significantieniveau 0,05 willen testen of de onderliggende verdeling normaal verdeeld is, moeten we in tabel 3.2 in de bijbehorende kolom naar de kritische waarden kijken. We benaderen het kritisch punt voor het gewenste percentage van ongecensureerde observaties d.m.v. lineaire interpolatie tussen 0,7 en 0,8: (y

y1 ) · (x2

x1 ) = (x

x1 ) · (y2

) (y 1, 3471) · (0.8 ) y = 1.35292

y1 )

0.7) = (0, 76

0.7) · (1.3568

1.3471)

wat groter is dan de waarde van Dt⇤ . Voor deze statistiek aanvaarden we de nulhypothese op het 0,05-niveau. Ook voor de andere statistieken ligt de waarde onder het kritisch punt zodat we ook voor deze de nulhypothese niet kunnen verwerpen op het 0,05-niveau (voor tabellen zie Stephens[3]).

Voor linksgecensureerde gegevens Bij het construeren van EDF-tests voor linksgecensureerde steekproeven herbruiken we de methode voor rechtsgecensureerde steekproeven door de de linksgecensureerde steekproef te transformeren naar een rechtsgecensureerde. Zij T1 , .., Tn een linksgecensureerde steekproef en Z1 , .., Zn de bijbehorende steekproef na Probability Integral Transformation, met vaste censureringstijd t. We veronderstellen dat de r grootste observaties ongecensureerd zijn. Z˜(i)

= 1 Z(n+1 t˜ = 1 t

i)

, i = 1, .., r

De waarden Z˜(i) zullen dan rechtsgecensureerd zijn want ze zullen allen kleiner zijn dan 1 t, zoals we zien op figuur 3.14.

44

Figuur 3.14: De transformatie van links- naar rechtsgecensureerde gegevens Om het algoritme voor het berekenen van EDF-statistieken voor linksgecensureerde steekproeven (zie B.6.4) te bekomen, moeten we aan het algoritme in B.6.3 enkel deze transformatie toevoegen. Voorbeeld 3.4.7 We beschouwen de normaal verdeelde steekproef uit B.10. Hiervan maken we een rechtsgecensureerde dataset met vaste censureringstijd gelijk aan 4. Hierop passen we het algoritme toe. > linksgecensEDF(cen(data,1),delta(data,1),1,2,3) D Dster Wkwadr Akwadr pct ongecensureerd [1,] 0.07525366 0.7715366 0.06821092 0.9543142 0.61

We testen d.m.v statistiek (Dt˜)⇤ op signigicantieniveau 0,05 of de onderliggende verdeling normaal verdeeld is. Wegens de transformatie naar rechtsgecensureerde gegevens mogen we dezelfde tabellen van kritische punten gebruiken. Op een analoge manier als in vb. 3.4.6 (d.m.v. lineaire interpolatie) vinden we dat de waarde van de teststatistiek kleiner is dan het bijbehorend kritisch punt 1,3237. We aanvaarden de nulhypothese dus op significantieniveau 0,05.

3.4.6

EDF-tests met ongekende parameters

Voor ongecensureerde gegevens In het vorige deel specifieerden we steeds parameters µ en ↵. Maar wanneer we geen voorkennis hebben over de steekproef is het moeilijk om hiervoor een goede keuze te maken. Dan gebruiken we schatters voor de parameters. 45

De nulhypothese H0 is dan: steekproef X1 , .., Xn komt uit verdeling N (µ, met één of meer parameters onbekend.

2)

We onderscheiden dan drie gevallen: • Case 1 : µ onbekend, • Case 2 : µ bekend,

2

2

¯ bekend ) µ ˆ=X

onbekend ) ˆ2 =

1 n

P

i (Xi

¯ en ˆ2 = • Case 3 : allebei onbekend: ) µ ˆ=X

1 n 1

µ)2

P

i (Xi

µ ˆ)2

¯ bedoelen we het steekproefgemiddelde. Merk ook het verschil op Met X tussen de schatters voor 2 in case 2 en 3: in case 2 maken we gebruik van de gekende waarde voor µ. In de praktijk is case 3 de belangrijkste situatie. Om Z(i) = F (X(i) ; ✓) te berekenen gebruiken we voor alle gevallen volgende substitutie. Berekening Z(i) met substitutie: 1. Bereken wi , i = 1..n zodat: • wi = (X(i) • wi = (X(i) • wi = (X(i)

¯ X)/ (case 1) µ)/s1 (case 2) ¯ X)/s (case 3)

2. Bereken Z(i) = (wi ), 1  i  n, waarbij de standaardnormale verdeling is.

de verdelingsfunctie van

3. Bereken de gewenste statistiek met formules 3.10 of 3.11. 4. Vergelijk de (aangepaste) statistiek met de kritische punten in de bijbehorende tabel. Wanneer de waarde van de (aangepaste) statistiek groter is dan het kritisch punt op niveau ↵ verwerpen we H0 op significantieniveau ↵. Voor cases 1 en 2 zullen er geen aangepaste statistieken nodig zijn omdat de kritische punten voor eindige steekproefgroottes zeer snel naar de asymptotische punten convergeren. Voor case 3 is dit niet het geval. Voor meer details en de bijbehorende tabellen en aangepaste statistieken: zie Stephens[3]. De algoritmes voor deze situaties vinden we in B.7,B.8 en B.9. Voor linksgecensureerde gegevens In het geval van linksgecensureerde gegevens kunnen we µ en schatten met de loglikelihoodmethode. En zetten we de linksgecensureerde steekproef om naar een rechtsgecensureerde zodat we volgende werkwijze kunnen toepassen.

46

Beschouw de geobserveerde steekproef T1 , .., Tn met vaste censureringstijd gelijk aan c. Laat r het aantal ongecensureerde observaties zijn, dit zijn dan de r grootste: T(n+1 i) met 1  i  r. Werkwijze linksgecensureerde gegevens case 3: 1. Bereken schattingen µ ˆ en ˆ voor de parameters met de loglikelihoodmethode. 2. Bereken voor de ongecensureerde observaties wi zodat wi =

T(i) µ ˆ , i=n+1 ˆ

r, .., n.

3. Bereken Z(i) = (wi ) voor i = n + 1 r, .., n en t = ((c µ ˆ)/ˆ ), waarbij de verdelingsfunctie van de standaardnormale verdeling is. 4. Deze zetten we om naar een rechtsgecensureerde steekproef: Z˜(i) = 1 Z(n+1 i) , i = 1, .., r en t˜ = 1 t 5. Bereken de gewenste statistiek voor deze rechtsgecensureerde gegevens met formule 3.15 voor Dt˜ of 3.16 voor Wt˜2 of A2t˜ . 6. Vergelijk de aangepaste statistiek met de juiste tabel van kritische punten. Deze vinden we in Stephens[3] en zijn van ongeveer dezelfde vorm als tabel 3.2. De steekproefgrootte wordt er nog extra in rekening gebracht, aangezien de nauwkeurigheid van de parameterschatting afhangt van de steekproefgrootte. Als de statistiek de waarde van het kritisch punt voor een gegeven ↵ overschrijdt, dan verwerpen we H0 op significantieniveau ↵ Het opstellen van het algoritme wordt aan de ge¨ınteresseerde lezer overgelaten.

47

Besluit In het eerste hoofdstuk gaven we een korte inleiding over survivalanalyse. We lieten zien welke functies een event time karakteriseerden. In hoofstuk twee kwamen de verschillende soorten van censurering aan bod. Aan de hand van voorbeelden voerden we begrippen als rechtse censurering type I en type II, linkse censurering en vaste en random censureringstijd in. In het laatste hoofdstuk bespraken we uitvoerig enkele goodness-of-fit tests voor het testen t.o.v. een normale verdeling. Aangezien we in het geval van gecensureerde steekproeven niet meer willen testen t.o.v. de geobserveerde verdeling pasten we hier enkel klassieke methoden aan aan censurering met vaste censureringstijd. We startten met de constructie van de likelihoodfunctie en hoe deze wijzigd wanneer er gecensureerde observaties voorkomen. Deze gebruikten we in het verdere verloop van het hoofdstuk voor het schatten van parameters bij een gecensureerde steekproef. In de volgde sectie werden de Chi-kwadraattest besproken en aangepast aan linksgecensureerde gegevens, ook werd er een simulatie uitgevoerd voor het significantieniveau van deze testen. In de derde sectie van dit hoofdstuk kwam een visuele en daardoor soms nogal objectieve goodness-of-fit test aan bod: de QQ-plot. We stelden deze onmiddellijk op voor het linksgecensureerde geval. We bespraken ook waarom we bij een QQ-plot voor een linksgecensureerde steekproef de parameters van de normale verdeling wel moeten schatten, in tegenstelling tot een ongecensureerde situatie waar we enkel moeten testen t.o.v. een standaardnormale verdeling. In de laatste sectie bespraken we testen gebaseerd op verschillen met de empirische verdelingsfunctie, ook wel EDF-testen genoemd. Deze zijn dikwijls sterker dan de Chi-kwadraat verdelingstest omdat daarin informatieverlies is door het groeperen in klassen. We voerden twee soorten EDFstatistieken in: supremumstatistieken en kwadratische statistieken waarvan de Kolmogorov statistiek de meest bekende was maar dikwijls veel zwakker dan de kwadratische statistieken A2 en W 2 . We bespraken Probability Integral Transformation of hoe we een steekproef 48

kunnen omzetten naar een Z-steekproef zodat alle statistieken die maat zijn voor het verschil tussen de EDF en de verdelingsfunctie hetzelfde resultaat resultaat geven als de statistieken voor deze Z-steekproef t.o.v. de uniforme verdeling op [0, 1]. Uitgaande van de definities berekenden we ook formules voor de statistieken die eenvoudiger te implementeren zijn zodat we daarna enkele testen konden construeren. We bespraken voor het volledig gespecifieerde geval de EDF-test voor ongecensureerde en rechtsgecensureerde gegevens uitvoerig en leidden hier de test voor het linksgecensureerde geval uit af. Tenslotte bespraken we de werkwijze voor de situatie met ongekende parameters voor een ongecensureerde en een linksgecensureerde steekproef.

49

Bibliografie [1]

Prof. dr. Roel Braekers, (2005), Kanstheorie en statistiek 1, cursus 2e bach wiskunde, 162-175

[2]

Prof. dr. Roel Braekers, (2003-2004) Modelling ethanol-induced sleeping time in inbred strains of mice: A repeated measures design with left-censored data, eindwerk Master of Sciences in Biostatistics

[3]

Ralph B. DAgostino, Michael A. Stephens, (1986), Goodness-of-fit techniques, Marcel Dekker, New York

[4]

John P. Klein, Melvin L. Moeschberger, (2003) Survival Analysis: Techniques for Censored and Truncated Data, 2nd edition, New York Springer, 1-138

[5]

Rupert G. Miller Jr., (1981) Survival Analysis, John Wiley & Sons, Inc.,1-80,164-175

[6]

Prof. dr. N. Veraverbeke, (2003-2004) Stochastische processen cursus 2e kan wiskunde,81-98

[7]

http://www.sys-ev.com/reliability01.htm : voor de bathtubcurve

[8]

L. Lee, (1980) , Testing adequacy of the Weibull and log linear rate models for a Poisson process, Technometrics, 22(2), 195-199 (bathtubcurve)

[9]

http://en.wikipedia.org/wiki/Bathtub curve

50

Bijlage A

Enkele definities en eigenschappen 1. De exponentiële verdeling :

> 0 en x

• survivalfunctie S(x) = e

0

x

• dichtheidsfunctie f (x) = e

x

• hazard rate functie h(x) = 2. De Weibull verdeling : ↵,

> 0 en x

0

x↵

• survivalfunctie S(x) = e

• dichtheidsfunctie f (x) = ↵ x↵

• hazard rate functie h(x) = ↵ x↵ 3. De log-logistische verdeling : ↵, • survivalfunctie S(x) =

x↵

1e 1

> 0 en x

0

1 1+ x↵

• dichtheidsfunctie f (x) =

↵x↵ 1 (1+ x↵ )2

• hazard rate functie h(x) = 4. De Normale verdeling : µ,

↵x↵ 1 1+ x↵

> 0 en x

• survivalfunctie S(x) = 1

• dichtheidsfunctie f (x) =

(x p1 2⇡

e

µ

0 )

1 x µ 2 ( ) 2

• hazard rate functie h(x) = µ 5. De Binomiale verdeling : p • pj =

n j

pj (1

p)n

j

(j = 0, .., n)

6. Definitie voorwaardelijke kans: P (A | B) = 51

P (A\B) P (B)

7. De formule voor de badkuipfunctie voorgesteld door Lee[7] hbt (x) = k x

1 exp(↵x)

met ↵,k

8. De continue uniforme verdeling : • verdelingsfunctie F (x) = • dichtheidsfunctie f (x) =

0,0

<1

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

1

+ 1[b,1) (x)

9. Voor een continue niet-negatieve R 1stochastische variable X met eindig eerste moment geldt: E [X] = 0 [1 F (x)] dx

52

Bijlage B

Algoritmes B.1 B.1.1

Kunstmatig censureren van simulaties Linkse censurering

Input: Een ongecensureerde dataset en een vaste censureringstijd c. Werkwijze: Alle observaties kleiner dan c worden vervangen door c. Output: De linksgecensureerde dataset met vaste censureringstijd c. cen <-function(reeks,c) { x <-vector(length=length(reeks)) for (i in 1:length(reeks)) { if(reeks[i]<=c) {x[i] <-c} else{x[i] <-reeks[i]} } x }

B.1.2

Rechtse censurering

Input: Een ongecensureerde dataset en een vaste censureringstijd c. Werkwijze: Alle observaties groter dan c worden vervangen door c. Output: De rechtsgecensureerde dataset met vaste censureringstijd c. cenrechts <-function(reeks,c) { x <-vector(length=length(reeks)) for (i in 1:length(reeks)) { if(reeks[i]<=c) {x[i] <-reeks[i]} else{x[i] <-c} } x }

53

B.2

Construeren van

B.2.1

Linkse censurering

voor simulaties

Input: Een dataset en een vaste censureringstijd c. Werkwijze: Voor alle linksgecensureerde observaties (kleiner of gelijk aan c) wordt i nul, anders één. Output: Indicatorvector . delta <-function(reeks,c) { d <-vector(length=length(reeks)) for (i in 1:length(reeks)) { if(reeks[i]<=c) {d[i] <-0} else{d[i] <-1} } d }

B.2.2

Rechtse censurering

Input: Een dataset en een vaste censureringstijd c. Werkwijze: Voor alle rechtsgecensureerde observaties (groter of gelijk aan c) wordt i nul, anders één. Output: Indicatorvector . deltarechts<-function(reeks, c){1-delta(reeks,c)}

B.3 B.3.1

Likelihood likelihoodconstructie

Input: Een linksgecensureerde steekproef, bijbehorende en vaste censureringstijd c. Werkwijze: Met formule 3.5. Hier berekenen we echter ln L ipv. ln L Output: ln L waarbij ln L de log-likelihoodfunctie van deze steekproef voor een normale verdeling is. loglikfct<-function(cendata,di,c) { x<-cen(cendata,c) loglik<-function(p){ mu<-p[1] sigma<-p[2] -sum(di*((-1/2)*((x-mu)/sigma)^2-(1/2)*log(2*pi)-log(sigma)) +(1-di)*log(pnorm((c-mu)/sigma)))} loglik }

54

B.3.2

Contourplot van de loglikelihoodfunctie

Input: Een linksgecensureerde steekproef, de bijbehorende , de vaste censureringstijd c en het bereik van de assen. Werkwijze: Functie contourlikl berekent een matrix met waarden van de loglikelihood voor verschillende parameters. Deze worden bepaald door alle verschillende combinaties van elementen van vectoren x en y. Hier werd er voor de keuze van x en y ervan uitgegaan dat de schatting zich rond (2, 3) zou bevinden. Output: Een contourplot van de log-likelihoodfunctie van deze steekproef voor een normale verdeling en een rood punt waar deze haar maximum bereikt. x <- seq(1, 3, by=0.2) y <- seq(2, 4, by=0.2) contourlikl<-function(x,y,cendata,di,c) { z <-matrix(0,nrow=length(x),ncol=length(y)) for(i in 1:length(x)) { for(j in 1:length(y)) { z[i,j]<-sum(di*((-1/2)*((cendata-x[i])/y[j])^2-(1/2)*log(2*pi)-log(y[j]))+ +(1-di)*log( pnorm((c-x[i])/y[j]) )) } } z } fxy <-contourlikl(x,y,cendata2a,delta(cendata2a,1),1) contour(x,y,fxy,xlab=’mu’,ylab=’sigma’) points(mu,sigma,pch=par(’19’),col=’red’)

B.4 B.4.1

Chi-kwadraat verdelingstest Chi-kwadraattest voor ongecensureerde gegevens

Input: Een ongecensureerde dataset en een significantieniveau. Werkwijze: Er wordt getest of deze dataset normaal verdeeld is met als geschatte parameters het steekproefgemiddelde en de steekproefvariantie, a.d.h.v. een Chi-kwadraattest met 10 klassen. Output: De geobserveerde waarden, teststatistiek, cut-o↵ waarde, p-waarde, de geschatte parameters en of de nulhypothese aanvaard wordt. chitest10k<-function(data,alpha)

55

{ mu <-mean(data) sigma <-sqrt(var(data)) a1 <-qnorm(0.1)*sigma+mu geobs1 <-sum((data<=a1)*1) a2 <-qnorm(0.2)*sigma+mu geobs2 <-sum((data<=a2)*1)-sum((data<=a1)*1) a3 <-qnorm(0.3)*sigma+mu geobs3 <-sum((data<=a3)*1)-sum((data<=a2)*1) a4 <-qnorm(0.4)*sigma+mu geobs4 <-sum((data<=a4)*1)-sum((data<=a3)*1) a5 <-qnorm(0.5)*sigma+mu geobs5 <-sum((data<=a5)*1)-sum((data<=a4)*1) a6 <-qnorm(0.6)*sigma+mu geobs6 <-sum((data<=a6)*1)-sum((data<=a5)*1) a7 <-qnorm(0.7)*sigma+mu geobs7 <-sum((data<=a7)*1)-sum((data<=a6)*1) a8 <-qnorm(0.8)*sigma+mu geobs8 <-sum((data<=a8)*1)-sum((data<=a7)*1) a9 <-qnorm(0.9)*sigma+mu geobs9 <-sum((data<=a9)*1)-sum((data<=a8)*1) geobs10 <-sum((data > a9)*1) geobs<-c(geobs1,geobs2,geobs3,geobs4,geobs5,geobs6,geobs7,geobs8,geobs9,geobs10) print(geobs) p <-0.1 verwa <-length(data)*p t <-sum((geobs-verwa)^2/verwa) print("teststatistiek is") print(t) print("cutoff") print(qchisq(1-alpha,7)) print("pwaarde is") p<-1-pchisq(t,7) print(p) if(p<=alpha) {print("We werwerpen de nulhypothese op significantieniveau alpha") print("Data komt niet uit normale verdeling met mu en sigma gelijk aan:")} else {print("Data komt uit normale verdeling met mu en sigma gelijk aan:")} print(mu) print(sigma) p }

56

B.4.2

Chi-kwadraattest voor linksgecensureerde gegevens

7 klassen Input: Een linksgecensureerde dataset, bijbehorende deltavector, significantieniveau en de vaste censureringstijd. Werkwijze: Er wordt getest of de onderliggende verdeling (van exacte event times) van deze dataset normaal is met als geschatte parameters het steekproefgemiddelde en de steekproefvariantie, a.d.h.v. een Chi-kwadraattest met 7 klassen, waarvan 1 met enkel gecensureerde observaties. Output: De p-waarde, de geschatte parameters en of de nulhypothese aanvaard wordt. chilinks<-function(data,di,alpha,c) { mu <-mean(data) sigma <-sqrt(var(data)) stwa<-c(mu,sigma) stwa mu<-nlm(loglikfct(data,di,c),stwa,hessian=TRUE,print.level=0)$estimate[1] sigma<-nlm(loglikfct(data,di,c),stwa,hessian=TRUE,print.level=0)$estimate[2] #voor c Pcens<-pnorm((c-mu)/sigma) gecens<-sum((data<=c)*1) #tss c en a1 a1 <-qnorm(Pcens+(1-Pcens)/6)*sigma+mu geobs1 <-sum((data<=a1)*1)-sum((data<=c)*1) #tss a1 en a2 a2 <-qnorm(Pcens+2*(1-Pcens)/6)*sigma+mu geobs2 <-sum((data<=a2)*1)-sum((data<=a1)*1) #tss a2 en a3 a3 <-qnorm(Pcens+3*(1-Pcens)/6)*sigma+mu geobs3 <-sum((data<=a3)*1)-sum((data<=a2)*1) #tss a3 en a4 a4 <-qnorm(Pcens+4*(1-Pcens)/6)*sigma+mu geobs4 <-sum((data<=a4)*1)-sum((data<=a3)*1) #tss a4 en a5 a5 <-qnorm(Pcens+5*(1-Pcens)/6)*sigma+mu geobs5 <-sum((data<=a5)*1)-sum((data<=a4)*1) #groter dan a5 geobs6 <-sum((data>=a5)*1) pvect<-c(Pcens,(1-Pcens)/6,(1-Pcens)/6,(1-Pcens)/6,(1-Pcens)/6, +(1-Pcens)/6,(1-Pcens)/6) n <-length(data) verwa<-n*pvect print("het verwachte aantal gecensureerde onder Ho") print(verwa[1]) print("andere verwachte aantallen")

57

print(verwa[2:7]) #teststat t = sum_(1^k)(geobs-verwa)^2/verwa

geobs= aantal in intv

geobs<-c(gecens,geobs1,geobs2,geobs3,geobs4,geobs5,geobs6) t <-sum((geobs-verwa)^2/verwa) t <= qchisq(1-alpha,4) print("pwaarde is") p<-1-pchisq(t,4) print(p) #met pwaarde if(p<=alpha) {print("We werwerpen de nulhypothese op significantieniveau alpha") print("de data komt niet uit normale verdeling met mu en sigma gelijk aan:")} else {print("niet significant (niv. alpha) verschillend van normale met mu en sigma")} print(mu) print(sigma) p }

11 klassen chilinks11<-function(data,di,alpha,c) { mu <-mean(data) sigma <-sqrt(var(data)) stwa<-c(mu,sigma) stwa mu<-nlm(loglikfct(data,di,c),stwa,hessian=TRUE,print.level=0)$estimate[1] sigma<-nlm(loglikfct(data,di,c),stwa,hessian=TRUE,print.level=0)$estimate[2] #voor c Pcens<-pnorm((c-mu)/sigma) gecens<-sum((data<=c)*1) #tss c en a1 a1 <-qnorm(Pcens+(1-Pcens)/10)*sigma+mu geobs1 <-sum((data<=a1)*1)-sum((data<=c)*1) #tss a1 en a2 a2 <-qnorm(Pcens+2*(1-Pcens)/10)*sigma+mu geobs2 <-sum((data<=a2)*1)-sum((data<=a1)*1) #tss a2 en a3 a3 <-qnorm(Pcens+3*(1-Pcens)/10)*sigma+mu geobs3 <-sum((data<=a3)*1)-sum((data<=a2)*1) #tss a3 en a4 a4 <-qnorm(Pcens+4*(1-Pcens)/10)*sigma+mu geobs4 <-sum((data<=a4)*1)-sum((data<=a3)*1) #tss a4 en a5 a5 <-qnorm(Pcens+5*(1-Pcens)/10)*sigma+mu geobs5 <-sum((data<=a5)*1)-sum((data<=a4)*1)

58

#tss a5 en a6 a6 <-qnorm(Pcens+6*(1-Pcens)/10)*sigma+mu geobs6 <-sum((data<=a6)*1)-sum((data<=a5)*1) #tss a6 en a7 a7 <-qnorm(Pcens+7*(1-Pcens)/10)*sigma+mu geobs7 <-sum((data<=a7)*1)-sum((data<=a6)*1) #tss a7 en a8 a8 <-qnorm(Pcens+8*(1-Pcens)/10)*sigma+mu geobs8 <-sum((data<=a8)*1)-sum((data<=a7)*1) #tss a8 en a9 a9 <-qnorm(Pcens+9*(1-Pcens)/10)*sigma+mu geobs9 <-sum((data<=a9)*1)-sum((data<=a8)*1) #groter dan a9 geobs10 <-sum((data>=a9)*1) pvect<-c(Pcens,(1-Pcens)/10,(1-Pcens)/10,(1-Pcens)/10,(1-Pcens)/10, +(1-Pcens)/10,(1-Pcens)/10,(1-Pcens)/10,(1-Pcens)/10, +(1-Pcens)/10,(1-Pcens)/10) n <-length(data) verwa<-n*pvect print("het verwachte aantal gecensureerde onder Ho") print(verwa[1]) print("het verwachte aantal andere") print(verwa[2:11]) #teststat t = sum_(1^k)(geobs-verwa)^2/verwa geobs<-c(gecens,geobs1,geobs2,geobs3,geobs4,geobs5,geobs6,geobs7, +,geobs8,geobs9,geobs10) t <-sum((geobs-verwa)^2/verwa) t <= qchisq(1-alpha,11-2-1) print("pwaarde is") p<-1-pchisq(t,11-2-1) print(p) #met pwaarde if(p<=0.05) {print("We werwerpen de nulhypothese op significantieniveau alpha") print("de data komt niet uit een normale verdeling met mu en sigma gelijk aan:")} else {print("niet significant (niv alpha) verschillend van normale met mu en sigma:")} print(mu) print(sigma) p }

B.5

QQ-plot

Input: Een linksgecensureerde dataset met vaste censureringstijd gelijk aan 1. Op voorhand moeten ook nog de parameters van de theoretische verdeling gespecifieerd worden als µ en . 59

Werkwijze: De geobserveerde waarden van de ordestatistieken worden uitgezet tegen de verwachte waarden. Deze laatste worden berekend zoals in sectie 3.3. Output: De QQ-plot waar getest wordt tegen een normale verdeling met parameters µ en . qqplot<-function(data) #om te testen tegen normaalverdeling { Fnorm<-function(x) pnorm(x,mu,sigma) Fc<-function(x) { Fc<-vector(length=length(x)) for(i in 1:length(x)) { if(x[i]<1)Fc[i]<-0 else Fc[i]<-Fnorm(x[i]) } Fc } Evect<-function(n,theoverdel) { E<-vector(length=n) for(k in 1:n) { arg<-function(x) { pbinom(k-1,n,theoverdel(x)) } E[k]<-integrate(arg,0,Inf)$value } E } n<-length(data) Evect(n,Fc) plot(Evect(n,Fc),sort(data),xlab="E[Z_(k)]",ylab="X_(i)",main="qq-plot") abline(0,1) }

B.6 B.6.1

EDF-statistieken case 0 Grafieken van de EDF en verdelingsfunctie

Input: Een ongecensureerde dataset en parameters µ en van de theoretische normale verdeling. Output: Een grafiek met de verdelingsfunctie van N (µ, 2 ) en de EDF van de steekproef. tekennorm<-function(data,mu,sigma)

60

{ datas<-sort(data) G<-function(x)pnorm(x,mu,sigma) plot(G,min(data),max(data),main="EDF en DF N(mu,sigma^2)") #verdelingsfct N(mu,sigma^2) lines(datas,(1:length(data))/length(data), type=’s’) #empirische verdelingsfct abline(1,0) }

Input: Een ongecensureerde dataset en parameters µ en van de theoretische normale verdeling. Output: Een grafiek met de uniforme verdelingsfunctie (op [0, 1]) en de EDF steekproef die we bekomen na de PIT. tekenunif<-function(data,mu,sigma) { datas<-sort(data) G<-function(x)pnorm(x,mu,sigma) Z<-G(datas) #PIT plot(Z,(1:length(data))/length(data),type=’s’, +main="EDF van Z_i tov uniforme verdeling",xlab="Z_(i)",ylab="i/n") abline(0,1) #empirische zjes tov uniforme (rechte) abline(1,0) }

B.6.2

EDF-testen voor ongecensensureerde gegevens

Input: Een ongecensureerde dataset en parameters µ en van de theoretische normale verdeling, zodat deze volledig gespecifieerd is Werkwijze: Het algoritme berekent de formules 3.10 en 3.11. Output: De EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . ongecensEDF<-function(data,mu,sigma) { n<-length(data) G<-function(x)pnorm(x,mu,sigma) Z<-G(sort(data)) Dplus<-max((1:n)/n-Z) Dmin<-max(Z-((1:n)-1)/n) D<-max(Dplus,Dmin) Wkwadr<-sum((Z-(2*(1:n)-1)/(2*n))^2)+1/(12*n) Akwadr<- -n-(1/n)*sum((2*(1:n)-1)*log(Z)+(2*n+1-2*(1:n))*log(1-Z)) EDF<-matrix(c(Dplus,Dmin,D,Wkwadr,Akwadr),byrow=T,ncol=5) EDFnames<-list(c(),c("Dplus","Dmin","D","Wkwadr","Akwadr")) dimnames(EDF)<-EDFnames

61

EDF #om EDFstatn tov normale verdeling met mu en sigma te berekenen }

Input: Een ongecensureerde dataset en parameters µ en van de theoretische normale verdeling, zodat deze volledig gespecifieerd is Werkwijze: Na het berekenen van de EDF-statistieken wordt formule 3.14 toegepast. Output: De aangepaste EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . ongecensEDFster<-function(data,mu,sigma) { n<-length(data) Dplusster<-ongecensEDF(data,mu,sigma)[1]*(sqrt(n)+0.12+0.11/sqrt(n)) Dminster<-ongecensEDF(data,mu,sigma)[2]*(sqrt(n)+0.12+0.11/sqrt(n)) Dster<-ongecensEDF(data,mu,sigma)[3]*(sqrt(n)+0.12+0.11/sqrt(n)) Wkwadrster<-(ongecensEDF(data,mu,sigma)[4]-0.4/n+0.6/(n*n))*(1.0+1.0/n) Akwadrster<-ongecensEDF(data,mu,sigma)[5] #want in orde noor n>=5 EDFster<-matrix(c(Dplusster,Dminster,Dster,Wkwadrster,Akwadrster),byrow=T,ncol=5) EDFsternames<-list(c(),c("Dplusster","Dminster","Dster","Wkwadrster","Akwadrster")) dimnames(EDFster)<-EDFsternames EDFster #om EDFsterstatn tov normale verdeling met mu en sigma te berekenen }

B.6.3

EDF-testen voor rechtsgecensensureerde gegevens

Input: Een rechtsgecensureerde dataset met vaste censureringstijd c en parameters µ en van de theoretische normale verdeling, zodat deze volledig gespecifieerd is Werkwijze: Berekening van de aangepaste statistiek van Dt en statistieken Wt2 en A2t d.m.v. formule 3.16. Output: De EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . rechtsgecensEDF<-function(cendata,deltadata,c,mu,sigma) { n<-length(cendata) G<-function(x)pnorm(x,mu,sigma) Z<-G(sort(cendata)) r<-sum(deltadata) t<-G(c) Znieuw<-vector(length=r+1) Znieuw<-c(Z[1:r],t) Zshort<-Z[1:r]

62

D<-max((1:r)/n-Zshort,Zshort-((1:r)-1)/n,t-r/n) Dster<-sqrt(n)*D+0.19/sqrt(n) Wkwadr<-sum((Znieuw-((2*(1:(r+1))-1)/(2*n)))^2)+(r+1)/(12*n*n)+(n/3)*(t-(r+1)/n)^3 Akwadr<-(-1/n)*sum((2*(1:(r+1))-1)*(log(Znieuw)-log(1-Znieuw))) -2*sum(log(1-Znieuw))-(1/n)*((r+1-n)*(r+1-n)*log(1-Znieuw[r+1]) -r*r*log(Znieuw[r+1])+n*n*Znieuw[r+1]) perc<-r/n uitk<-matrix(c(D,Dster,Wkwadr,Akwadr,perc),byrow=T,ncol=5) uitknames<-list(c(),c("D","Dster","Wkwadr","Akwadr","pct ongecensureerd")) dimnames(uitk)<-uitknames uitk #om EDFstatn tov verdeling met mu en sigma te berekenen voor rechtsgecens gegns }

B.6.4

EDF-testen voor linksgecensensureerde gegevens

Input: Een linksgecensureerde dataset met vaste censureringstijd c en parameters µ en van de theoretische normale verdeling, zodat deze volledig gespecifieerd is. Werkwijze: Berekening van de aangepaste statistiek van Dt˜ en statistieken Wt˜2 en A2t˜ d.m.v. een transformatie naar een rechtsgecensureerde steekproef. Output: De EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . linksgecensEDF<-function(cendata,deltadata,c,mu,sigma) { n<-length(cendata) G<-function(x)pnorm(x,mu,sigma) Z<-G(sort(cendata)) r<-sum(deltadata) t<-1-G(c) Zster<-vector(length=r) Zster<-Z[(n-r+1):n] Zster<-1-sort(Zster,decreasing=T) Znieuw<-c(Zster,t) D<-max((1:r)/n-Zster,Zster-((1:r)-1)/n,t-r/n) Dster<-sqrt(n)*D+0.19/sqrt(n) Wkwadr<-sum((Znieuw-((2*(1:(r+1))-1)/(2*n)))^2)+(r+1)/(12*n*n)+(n/3)*(t-(r+1)/n)^3 Akwadr<-(-1/n)*sum((2*(1:(r+1))-1)*(log(Znieuw)-log(1-Znieuw))) -2*sum(log(1-Znieuw))-(1/n)*((r+1-n)*(r+1-n)*log(1-Znieuw[r+1]) -r*r*log(Znieuw[r+1])+n*n*Znieuw[r+1]) perc<-r/n uitk<-matrix(c(D,Dster,Wkwadr,Akwadr,perc),byrow=T,ncol=5) uitknames<-list(c(),c("D","Dster","Wkwadr","Akwadr","pct ongecensureerd"))

63

dimnames(uitk)<-uitknames uitk #om EDFstatn tov normale met mu en sigma te berekn voor linksgecens gegns }

B.7

EDF-statistieken case 1 voor ongecensureerde gegevens

Input: Een ongecensureerde dataset en parameter van de theoretische normale verdeling, zodat deze volledig gespecifieerd is Werkwijze: Na het berekenen van de EDF-statistieken wordt formule 3.14 toegepast. Output: De EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . ongecensEDFcase1<-function(data,sigma) { #sigma gekend, mu ongekend n<-length(data) mu<-(1/n)*sum(data) data<-sort(data) w<-(data-mu)/sigma Z<-pnorm(w) #pnorm geeft verdelingsfct #Dplus<-max((1:n)/n-Z) #Dmin<-max(Z-((1:n)-1)/n) #D<-max(Dplus,Dmin) Wkwadr<-sum((Z-(2*(1:n)-1)/(2*n))^2)+1/(12*n) Akwadr<- -n-(1/n)*sum((2*(1:n)-1)*log(Z)+(2*n+1-2*(1:n))*log(1-Z)) EDF<-matrix(c(Wkwadr,Akwadr),byrow=T,ncol=3) EDFnames<-list(c(),c("Wkwadr","Akwadr")) dimnames(EDF)<-EDFnames print("geschatte mu") print(mu) EDF #om EDFstatn tov normale }

met onbek mu en bekende sigma te berekenen

64

B.8


Input: Een ongecensureerde dataset en parameter µ van de theoretische normale verdeling, zodat deze volledig gespecifieerd is Werkwijze: Na het berekenen van de EDF-statistieken wordt formule 3.14 toegepast. Output: De aangepaste EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . ongecensEDFcase2<-function(data,mu) { #mu gekend, sigma ongekend n<-length(data) sigmakwadr<-(1/n)*sum((data-mu)^2) #want de var in R geeft 1/(n-1), want schat mu nog data<-sort(data) w<-(data-mu)/sqrt(sigmakwadr) Z<-pnorm(w) #pnorm geeft verdelingsfct Wkwadr<-sum((Z-(2*(1:n)-1)/(2*n))^2)+1/(12*n) Akwadr<- -n-(1/n)*sum((2*(1:n)-1)*log(Z)+(2*n+1-2*(1:n))*log(1-Z)) EDF<-matrix(c(Wkwadr,Akwadr),byrow=T,ncol=3) EDFnames<-list(c(),c("Wkwadr","Akwadr")) dimnames(EDF)<-EDFnames print("geschatte sigma") print(sqrt(sigmakwadr)) EDF #om EDFstatn tov normale met onbek mu en bekende sigma te berekenen }

B.9


Input: Een ongecensureerde dataset. Werkwijze: Na het berekenen van de EDF-statistieken wordt formule 3.14 toegepast. Output: De aangepaste EDF-statistieken voor deze dataset voor een verdelingstest met theoretische normale verdeling met parameters µ en . ongecensEDFcase3<-function(data) { #alle param ongekend n<-length(data) mu<-(1/n)*sum(data) sigmakwadr<-(1/(n-1))*sum((data-mu)^2)

65

data<-sort(data) w<-(data-mu)/sqrt(sigmakwadr) Z<-pnorm(w) #pnorm geeft verdelingsfct Dplus<-max((1:n)/n-Z) Dmin<-max(Z-((1:n)-1)/n) D<-max(Dplus,Dmin) Wkwadr<-sum((Z-(2*(1:n)-1)/(2*n))^2)+1/(12*n) Akwadr<- -n-(1/n)*sum((2*(1:n)-1)*log(Z)+(2*n+1-2*(1:n))*log(1-Z)) Dster<-D*(sqrt(n)-0.01+0.85/sqrt(n)) Wkwadrster<-Wkwadr*(1+0.5/n) Akwadrster<-Akwadr*(1+(0.75/n)+(2.25/(n*n))) EDF<-matrix(c(Dplus,Dmin,D,Wkwadr,Akwadr,Dster,Wkwadrster,Akwadrster), byrow=T,ncol=10) EDFnames<-list(c(), c("Dplus","Dmin","D","Wkwadr","Akwadr","Dster","Wkwadrster","Akwadrster")) dimnames(EDF)<-EDFnames print("geschatte mu") print(mu) print("geschatte sigma") print(sqrt(sigmakwadr)) EDF #om EDFstatn tov normale verdeling met onbek mu en }

B.10

sigma te berekenen

Gebruikte data

In dit deel bevinden zich de datasets die in de voorbeelden gebruikt werden. Deze data komt uit een normale verdeling met µ = 2 en

= 3:

data<-c(1.84513612,5.44572642,-0.87601091,-3.42768277,4.25569255,0.15233678, 10.76629967, 1.69936886, 1.54534191, 2.58541346, 3.69206743, 4.08872675, -0.31922706, 2.26565922, 1.41843088, 0.20615212, 3.42359721, 2.06647916, 2.11091601, 2.74631771, 3.28728209, 0.39965569, 4.61853127, 1.15733454, -0.06511262, 0.74334873, 1.79541485, 1.88721725, 7.12278313, -2.98314629, -0.20041048, 1.74175380, 11.12450951, -1.25010944, 3.04450157, 4.65842532, -0.70476539, 2.54178948, 0.65089573, -6.03752387, 1.07387780, 5.75667301, 2.20562184, -1.73800510, 1.98689564, -1.34444796, 6.02982133, 0.35507592, -2.21412885, 1.85941637, -3.40645253, -0.13219305, 1.04489677, 0.12643212, -1.76360092, 0.87470514, 1.74870460, -0.04490868, 7.05302555, 1.78749684, 3.04646997, -0.54195618, 1.90914327, 6.75994902, 3.65808570, 2.43288601, 0.44157121, -0.49440859, -4.74123958, 5.00689397, 5.61972199, -0.15563538, 1.92434639, 6.51210247, 4.24233620, -1.10119084, 3.77817458, -4.41894989, -0.19503947, 0.81052497, -1.88998982, 3.62415648, -1.16723139, 6.89276665, 4.19060323, 1.09393475, 3.39943636, 3.23702779, 1.97867974, -3.76432297, 5.39141932, 0.19819723, 5.01232856, 8.56912984, -2.33134499, 4.71189708, 9.05188265, 5.96865280, -0.49158922, 2.67207838)

66

Deze data komt uit een weibull verdeling met parameters ↵ = 1/2 en  = 3 dataweib<-c(1.842478e-02, 1.455713e+00, 7.048937e-01, 3.499931e+00, 3.256789e-03,1.103787e+00, 3.450415e+00, 1.130465e-01, 6.118869e-02, 3.780170e-01, 5.588736e-01, 7.117366e+00, 9.424159e-01, 5.461072e-01, 2.720056e+01, 1.841125e+00, 8.894903e+00, 4.769068e-02, 2.191561e+00, 5.096049e-01, 1.851164e+01, 2.069112e-02, 3.163609e-01, 5.402820e-02, 9.201491e+00, 3.615361e-01, 7.933895e-01, 2.602043e+00, 1.281132e-03, 5.284573e-01, 5.491591e+00, 1.351198e+01, 1.047990e+01, 3.507491e+00, +5.823014e+00, 1.246389e+00, 1.677134e-04, 5.896115e-02, 4.873803e+00, 1.948938e-01, 5.777831e-01, 6.462059e-01, 5.992197e-01, 1.417918e+00, 1.325229e-02, 4.030842e+00, 3.317959e+00, 1.883690e+00, 3.625706e-02, 1.324459e-01, 3.979752e+00, 6.525171e+00, 1.042883e+01, 7.064847e+00, 2.058648e+00, 6.476457e-01, 9.938603e-02, 7.344764e-01, 2.197217e+00, 1.842943e-03, 8.071963e-02, 2.008953e-02, 2.596304e-01, 6.160634e-01, 2.749551e+00, 2.127994e+01, 4.162498e+00, 1.466304e+01, 1.126960e+01, +9.154395e-02, 1.130415e+01, 2.096326e-01, 2.733877e+00, 2.179884e+00, 4.316987e-02, 1.635736e+01, +3.393410e-02, 1.084282e+01, 5.198009e-01, 1.851786e-01, 4.645803e-01, 2.032836e-01, 1.885048e-02, 1.209044e+00, 1.283897e+01, 1.947998e+00, 1.311148e+01, 1.552580e+00, 2.634471e-03, 1.780850e+00, 5.549896e-01, 4.769368e-01, 4.650112e-03, 3.093434e-01, 6.337369e-01, 3.166860e-02, 1.603347e+00, 5.765788e-03, 2.851725e+00, 9.966618e-01)

Dit is de Leghorn Chicken dataset die gebruikt wordt in sectie 3.4. Deze data vonden we in Stephens[3] pagina 98. X<-c(156,162,168,182,186,190,190,196,202,210,214,220,226,230,230,236,236,242,246,270)

67

Goodness-of-fit bij linksgecensureerde observaties. Sara Appeltants

Recommend Documents