Bacheloreindproject Statistische Analyse Overkomstpercentage Poststukken

Bacheloreindproject Statistische Analyse Overkomstpercentage Poststukken Fieke Taal 1236466 31 maart 2011

Inhoudsopgave 1 Inleiding

4

2 Achtergrond informatie 2.1 TNT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 5 7

3 Introductie gegeneraliseerde lineaire modellen 3.1 Het lineaire model . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Het gegeneraliseerde lineaire model . . . . . . . . . . . . . . . 3.3 Een speciaal geval: de Bernoulli-verdeling . . . . . . . . . . .

9 9 10 12

4 Analyse van gegeneraliseerde lineaire modellen met behulp van p-waarden 4.1 Methoden voor modelselectie binnen gegeneraliseerde lineaire modellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Fitten van een logistisch regressie model in R . . . . . . . . . 4.2.1 Kleine dataset . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Freedman Experiment . . . . . . . . . . . . . . . . . . 4.2.3 Nabootsing dataset TNT . . . . . . . . . . . . . . . . 4.3 Conclusie vooronderzoek . . . . . . . . . . . . . . . . . . . . . 5 Bayesiaanse theorie 5.1 Termen en Bayesiaanse notatie 5.2 Bayes Theorema . . . . . . . . 5.3 Bayesian Information Criterion 5.4 Voorbeeld . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

6 Onderzoek dataset TNT 6.1 Analyse in Bayesian Model Averaging 6.2 Opzet onderzoek . . . . . . . . . . . . 6.3 Resultaat . . . . . . . . . . . . . . . . 6.4 Testen data 2008 . . . . . . . . . . . . 6.5 Conclusie onderzoek TNT . . . . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

14 15 15 17 19 22

. . . .

23 23 23 24 26

. . . . .

34 34 34 35 42 42

7 Conclusie A Script vooronderzoek, datasets A.1 Kleine dataset . . . . . . . . . A.1.1 Test 1 . . . . . . . . . A.1.2 Test 2 . . . . . . . . . A.1.3 Test 3 . . . . . . . . . A.2 Freedman Experiment . . . . A.3 Nabootsing TNT dataset . .

14

44 en . . . . . . . . . . . .

resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

45 45 46 47 49 52 52

A.4 Voorbeeld Bayes . . . . . . . . . A.4.1 Berekening MLE . . . . . A.4.2 Script in R alleen B21 . . A.4.3 Script in R, B21 met BIC B Script onderzoek

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

53 53 54 58 63

3

1

Inleiding

Bij TNT is men ge¨ınteresseerd in het overkomstpercentage van poststukken. Dit percentage is gedefinieerd als het aantal poststukken dat binnen een dag op het correcte afleveradres wordt bezorgd ten opzichte van het totaal aantal verzonden poststukken. Dit wordt gemeten door een steekproef, waarbij ongeveer 200 poststukken worden gevolgd. Een belangrijke reden om dit overkomstpercentage te meten schuilt in het feit dat de OPTA (een wettelijk toezichtsorgaan) TNT verplicht om een percentage van tenminste 95% te halen. Ook door de toenemende concurrentie is TNT gebaat bij een zo goed mogelijke presentatie naar de consument toe. Van verschillende factoren wordt vermoed dat zij invloed hebben op het overkomstpercentage. Hierbij moet gedacht worden aan bijvoorbeeld het gewicht of formaat van de envelop, is het poststuk handmatig of machinaal verwerkt, maakt de dag van de week uit, etc. Het doel van dit onderzoek is om met behulp van historische data de belangrijkste factoren te achterhalen die het overkomstpercentage be¨ınvloeden. Hierbij zal gebruik gemaakt worden van logistische regressie, waarbij naar verschillende methoden voor modelselectie wordt gekeken. In deze scriptie zal eerst in hoofdstuk 2 een korte beschrijving gegeven worden van hoe de structuur binnen TNT eruit ziet en hoe het sorteerproces in zijn werk gaat. Ook wordt in dit hoofdstuk de gegeven data beschreven. In hoofdstuk 3 zal de theorie achter lineaire modellen en gegeneraliseerde lineaire modellen beschreven worden. Als voorbeeld voor deze theorie zal de verwachting en variantie van de Bernoulli-verdeling uitgewerkt worden. Vervolgens wordt in hoofdstuk 4 een aantal methoden van modelselectie beschreven. In dit hoofdstuk wordt ook een drietal experimenten gedaan om kennis te maken met het statistisch software pakket R om te achterhalen welke problemen zich voordoen bij de verschillende methoden van analyse van (gegeneraliseerde) lineaire modellen. In hoofdstuk 5 wordt een deel van de Bayesiaanse theorie uitgelegd en aan de hand van een voorbeeld wordt gekeken naar de resultaten van deze Baysian Model Averaging (BMA) methode en het Bayesian Information Criterion (BIC) in R. Hierna zal in hoofdstuk 6 de theorie toegepast worden op de data van TNT en zullen de resultaten weergegeven worden. Er wordt afgesloten met een conclusie.

4

2

Achtergrond informatie

In dit hoofdstuk zal enige achtergrond informatie gegeven worden over TNT en het sorteerproces van de post. Ook wordt de historische data die gebruikt is uitgelegd en worden enkele variabelen verder toegelicht.

2.1

TNT

TNT Post is ontstaan op 15 januari 1799, wanneer de posterijen die er tot dusver waren worden omgevormd tot een nationale onderneming. Het heette toen nog PTT Post, Staatsbedrijf der Posterijen, Telegrafie en Telefonie. Later werd dit Koninklijke PTT Nederland (KPN) en in 1996 nam de KPN het Australische bedrijf ’Thomas Nationwide Transport (TNT)’ over en werd het TNT Post Groep oftewel TPG Post. Pas in 2006 wordt de naam wegens concurrentie in binnen en buitenland veranderd in TNT Post. De letters hebben geen betekenis meer, maar TNT is nog wel wereldwijd een merknaam. De diensten die TNT biedt kunnen opgedeeld worden in twee delen: Post en Express. Het Post gedeelte kan ook weer opgesplitst worden in: Mail Nederland, de verspreidingsorganisatie van de post in Nederland, European Mail Networks, het Europese netwerk van TNT, Cendris, het dienstenaanbod op het gebied van marketing en data en documentmanagement en Spring, het internationale samenwerkingsverband. Onder Express vallen de wereldwijde koeriersdiensten. In dit Bachelorproject is alleen het Post gedeelte van belang en dan in het bijzonder de verspreiding van de post. Nederland is verdeeld in zes regio’s, elk met een eigen sorteercentrum. Deze staan in Den Haag, Rotterdam, Amsterdam, Zwolle, Nieuwegein en Den Bosch. Zodra een poststuk in de postbus wordt gegooid, begint het collectie-sorteer-bezorgproces. Alle poststukken worden eerst verzameld en naar een van de 87 Businessbalies gebracht. Onder poststukken wordt alle post verstaan die door de brievenbus past, dus ook zakelijke post. In de Businessbalies worden de poststukken, pakketten of aangetekende post klaargemaakt om naar de sorteercentra gebracht te worden. Poststukken, pakketten en aangetekende post gaan naar aparte centra. In het sorteercentrum worden de poststukken door verschillende machines gesorteerd op grootte, worden de postzegels gestempeld, adressen worden gelezen en de post wordt op postcode gesorteerd. Dit gebeurt in twee sorteringen. Tijdens de eerste sortering wordt er per sorteercentra op 180 richtingen, dus gemeenteniveau, gesorteerd. Aangezien er zes regio’s zijn, betekent dit 30 richtingen per regio. Na de eerste sortering worden de verschillende richtingen per vrachtwagen naar de andere sorteercentra gebracht en een aantal richtingen blijft in het sorteercentrum. Vervolgens wordt er tijdens de tweede sortering weer op 180 richtingen, dus wijkniveau, gesorteerd. Hierna worden de gesorteerde 5

poststukken gebracht naar een van de 180 voorbereidingsgebieden per regio. Hier worden ze verder gesorteerd op loop-en huisnummervolgorde, waarna ze door de postbode rondgebracht worden. Dagelijks gaat het om 16 miljoen poststukken, waarvan 92% zakelijk en dus slechts 8% van de consument, voor 7,6 miljoen huishoudens en bedrijven. Dit zijn dus 644 poststukken per huishouden per jaar!

Figuur 1: Envelop met Pric en SIX

Tijdens het sorteerproces wordt er door de machine een zogenaamde Pric-code en SIX-code op het poststuk gespoten. Beide codes staan rechtsonder op het poststuk aan de voorkant, zie figuur (1). Hierbij is de Pric-code de zwarte reeks letters en nummers die maximaal uit 16 karakters bestaat, in dit geval dus R13VV 3062CZ300. Deze wordt tijdens de eerste sortering op het poststuk aangebracht. De eerste letter is het sorteercentrum waar het poststuk is geweest, dit kan zijn Amsterdam (A), ’s Hertogenbosch (H), ’s Gravenhage (G), Nieuwegein (N), Rotterdam (R) en Zwolle (Z). In dit geval is dit dus Rotterdam. De volgende 2 cijfers geven aan door welke Sorteermachine het poststuk is gegaan, in dit geval 13. De 2 a 3 letters die volgen staan voor de wijze waarop het poststuk herkend is, dit is bijvoorbeeld door de Klant Index code (KIX, de zwarte code op de adressticker) of door Video Coding Recognizer. In dit geval is het adres door de KIX herkend. Het laatste deel van de code bestaat uit de vier postcodecijfers en twee letters en het huisnummer van de geadresseerde. In dit geval is de postcode dus 3062CZ met huisnummer 300. De letter die voor de Pric-code staat heeft te maken met welk proces het poststuk tijdens de tweede sortering heeft doorlopen, namelijk gesorteerd met behulp van de Sorteer Index (I), Automatisch codeermiddel (A) of vi6

deocodering (V). Als het niet mogelijk is om het poststuk te sorteren omdat de codeerinformatie ontbreekt, wordt een ’N’ geprint. De fluorescerende code onder de Pric-code wordt de Sorteer Index code (SIX) genoemd en wordt ook tijdens de eerste sortering op het poststuk aangebracht. Deze code is voor elk huishouden uniek. Dankzij de SIX-code kan de sorteerinformatie in de tweede sortering en daarna bij de huisnummersortering foutloos worden teruggelezen en gesorteerd worden. Elk jaar worden er door verschillende onderzoeksbureaus ongeveer 60000 proefbrieven van verschillend formaat en gewicht rondgestuurd. Bij ongeveer 30000 van deze brieven zit ook een transponder, een soort chip, zodat gekeken kan worden hoe laat en waar de brief is in een sorteercentrum en of de brief op tijd aangekomen is of vertraging heeft gehad. Hieruit blijkt dat 96,6% van de post binnen een dag bezorgd wordt, waarbij de 95% norm die gesteld is door het wettelijke orgaan OPTA ’ruim’ gehaald wordt. Uiteraard is er nog ruimte voor verbetering om ook de laatste 3,4% op tijd te laten aankomen.

2.2

Data

De data die gebruikt is bij dit onderzoek is verstrekt door TNT en bevat informatie over de proefbrieven met transponder die in 2007 zijn rondgestuurd. Deze data bestaat uit 66 variabelen, maar deze zijn niet allemaal even relevant voor dit onderzoek. Onder deze variabelen moet gedacht worden aan het formaat van de brief, gewicht, dag van verzenden en ontvangen, verzonden en ontvangen postcode, welk sorteercentrum, is het adres hand of machinaal geschreven, staat er een stempel op, is het poststuk op tijd of vertraagd, etc. Definieer  0 als het i-de poststuk vertraagd is Yi = 1 als het i-de poststuk op tijd is dus heeft een Bernoulli-verdeling, Yi ∼ Ber(p) . In overleg met TNT is een aantal variabelen gekozen, waarvan men denkt dat ze invloed hebben op het overkomstpercentage van de poststukken. Deze zijn hieronder verder uitgelegd. Gewicht: Het gewicht van de brief in gram. Dit kan zijn 0-20, 20-50, 50100 en 100-2000. Aangezien in R deze waarden problemen geven bij het analyseren van het model, zijn de gemiddelde waarden genomen. Dus 10, 35, 75 en 1000. Formaat: Het formaat van de brief. Dit kan de waarde A4, C5, C6 en EA5 zijn. 7

Schrift: Het schrift van de brief. Dit is met de machine of met de hand geschreven. Verzenden DvdW: De dag van de week dat het poststuk is verzonden. We onderscheiden hier drie soorten: vrijdag verstuurd, zaterdag/zondag verstuurd en de rest van de week verstuurd. Deze verdeling heeft te maken met het feit dat er op zaterdag andere postbezorgers (vaak studenten) rondlopen dan de rest van de week. Uiteraard is er ook verschil in de hoeveelheid post die tijdens het weekend wordt gesorteerd en dus maandag door de brievenbus valt en die van de rest van de week. InterofIntra: Inter betekent dat de poststukken tussen de 1e en 2e sortering van sorteercentrum zijn veranderd. Intra betekent dat het poststuk tijdens de twee sorteringen in hetzelfde sorteercentrum is geweest. Poststroom: Deze bestaat uit drie soorten, namelijk groot, klein en buspakjes. Onder groot wordt A4 formaat verstaan en klein is C5, C6 en EA5. Buspakjes zijn alle pakjes die wel door de brievenbus passen, maar die te dik of te zwaar zijn voor de andere twee stomen. Deze drie stromen gaan door verschillende machines. Streepjescode: Hiermee wordt de roze SIX-code bedoelt die rechtsonderaan op het poststuk staat. Alleen poststukken die tijdens de 1e sortering door de machine zijn gegaan, hebben een SIX-code. InfoPRIC: Deze variabele geeft informatie over de Pric-code, die besproken is in paragraaf 2.2. De waarden die de variabele kan aannemen zijn: ’geen’, ’nietgemeten’, ’aanwjuist’ en ’aanwonjuist’. Deze waarden staan respectievelijk voor: geen Pric-code aanwezig, Pric-code is niet gemeten (NA in de dataset), Pric-code is aanwezig en juist, Pric-code is aanwezig, maar onjuist of onleesbaar. Om te bepalen welke factoren invloed hebben op het sorteerproces, wordt er gezocht naar een verband tussen de invloedsfactoren en P (Y = 1), de kans dat een poststuk op tijd op bestemming is. Dit kan gedaan worden met logistische regressie wat een bijzonder geval is van het gegeneraliseerde lineaire model.

8

3

Introductie gegeneraliseerde lineaire modellen

Bij dit onderzoek wordt er gebruik gemaakt van logistische regressie. In dit hoofdstuk zal de theorie achter de modellen die gebruikt zijn worden uitgelegd en er zal blijken dat er ook tegen enkele problemen wordt aangelopen.

3.1

Het lineaire model

Stel dat er data is (x1 , Y1 )...(xn , Yn ) met {xi }i vast en {Yi }i random variabelen. Hiervoor wordt een model gezocht dat de data het beste fit. Veronderstel het lineaire model Yi = β0 + β1 xi + i . Hierbij is • Yi de respons variabele of ook wel afhankelijke variabele. • xi de instelvariabele of ook wel onafhankelijke variabele. • β0 en β1 zijn onbekend, waarbij β0 de intercept is. • i de ruis die bij het model hoort met verdeling N (0, σ 2 ) en σ 2 onbekend en 1 ...n onafhankelijk. Volgens dit model geldt E(Yi ) = µi = β0 + β1 xi . Het model heet lineair omdat het lineair is in de parameter β. Als er meerdere variabelen zijn (bijvoorbeeld xi1 en xi2 ) is het overzichtelijker om het model in matrixvectorvorm te zetten: Y = Xβ +  Hierbij is X de modelmatrix en Y,    Y1 1 x11 x12  Y2   1 x21 x22     ..  =  .. .. ..  .   . . . Yn 1 xn1 xn2

β en  vectoren, uitgeschreven     . . . x1m β0 1     . . . x2m    β1   2  ..   ..  +  ..  .  .   .  . . . xnm

βn

n

Er wordt nu een β gezocht die het model zo goed mogelijk fit met de data die gegeven zijn. Er wordt hier gebruikt gemaakt van het kwadraat van de som van de residuen, S=

n X

(Yi − (Xβ)i )2 = kY − Xβk2

i=1

Het verschil tussen de werkelijke waarde en de verwachte waarde moet zo ˆ de klein mogelijk zijn om een goed fittende β te vinden. Aangezien X β 9

orthogonale projectie is van Y op de lineaire ruimte R met R = {Xβ : β ∈ Rp }, geldt dat Y − Xβ orthogonaal is op iedere vector in R en dus ˆ =0 X t (Y − X β) Omschrijven geeft ˆ = X tY X tX β Als voor de modelmatrix X geldt dat deze volle kolomrang heeft en dus is X t X inverteerbaar, dan ˆ = (X t X)−1 X t Y β

3.2

Het gegeneraliseerde lineaire model

Bij lineaire modellen wordt de respons gemodelleerd met een normale verdeling. Gegeneraliseerde Lineaire Modellen (GLM) zijn een uitbreiding op lineaire modellen waarbij de responsvariabele Yi uit de exponenti¨ele familie van verdelingen komt. Hiertoe behoren bijvoorbeeld de Poisson-verdeling, Gamma-verdeling en Bernoulli-verdeling. In deze paragraaf wordt uitgelegd hoe gegeneraliseerde lineaire modellen zijn opgebouwd. In het statistische pakket R zijn methoden ge¨ımplementeerd waardoor het eenvoudig wordt om met deze GLM’s te werken. Het GLM werd in 1972 door John Nelder en Robert Wedderburn geformuleerd. Dit model heeft een basisstructuur van g(µi ) = Xti β met • µi = E(Yi ) • g is een gladde monotone ’link-functie’ • X is de modelmatrix, gedefinieerd in paragraaf 3.1 • β is een vector met onbekende parameters Een verdeling behoort tot de exponenti¨ele familie van verdelingen als zijn kansdichtsheidfunctie te schrijven is als   yθ − b(θ) fθ (y) = exp + c(y, φ) (1) a(φ) In vergelijking (1) zijn b(θ), a(φ) en c(y, φ) functies. φ is een schalingsparameter en zal vaak afhangen van de σ. θ hangt af van de modelparameters β en wordt de canonieke parameter genoemd.

10

Het is mogelijk om algemene uitdrukkingen te vinden voor de verwachting en variantie voor de exponenti¨ele familie van verdelingen. Stel Y is een random variabele met een dichtheid uit de exponenti¨ele familie zoals in vergelijking (1), dan wordt de loglikelihood gegeven door l(θ|y) = log(fθ (y)) dan is yθ − b(θ) l(θ|y) = + c(y, φ) a(φ) en y − b0 (θ) ∂ l(θ|y) = (2) ∂θ a(φ) In het algemeen geldt dat Z 0 Z Z fθ (y) ∂ ∂ 0 0 E(l (θ|Y )) = fθ (y)dy = fθ (y)dy = fθ (y)dy = 1=0 fθ (y) ∂θ ∂θ Hierbij is aangenomen dat differentiatie en integratie verwisseld mag worden f 0 (y) ∂ log(fθ (y)) = fθθ (y) . Dus en is l0 (θ|y) = ∂θ E(l0 (θ|Y )) =

E(Y ) − b0 (θ) =0 a(φ)

Waaruit blijkt dat E(Y ) = b0 (θ)

(3)

Deze vergelijking is de sleutel tot het verband tussen de modelparameters β en de canonieke parameter θ. Als µi = E(Yi ) dan geldt µi = g −1 (Xti β) maar ook µi = b0 (θ) en dus zijn de canonieke parameter θ en de parameter β gerelateerd volgens b0 (θ) = g −1 (Xti β) Uiteraard kan ook een algemene formule voor de variantie afgeleid worden door de tweede afgeleide van de loglikelihood te nemen ∂2 b00 (θ) l(θ|y) = − ∂θ2 a(φ) 2

∂l(θ|Y ) 2 ∂ Aangenomen dat E( ∂θ 2 l(θ|Y )) = −E[( ∂θ ) ] geldt met behulp van vergelijking (2) en (3)

b00 (θ) E[(Y − b0 (θ))2 ] = − a(φ) a(φ)2 b00 (θ)a(φ) = E[(Y − b0 (θ))2 ] = E[(Y − E(Y ))2 ] −

= Var(Y )

(4) 11

3.3

Een speciaal geval: de Bernoulli-verdeling

Aangezien er in dit onderzoek wordt gekeken naar poststukken die of op tijd of vertraagd zijn, dus Yi ∼ Ber(pi ) waarbij pi = E(Yi ) en Var(Y ) = pi (1 − pi ) zal voor deze verdeling worden aangetoond dat deze behoort tot de exponenti¨ele familie van verdelingen. Ook zal bewezen worden dat de algemene formules voor de verwachting en variantie kloppen. De Bernoulliverdeling heeft kansdichtheidsfunctie fp (y) = py (1 − p)1−y ,

y ∈ {0, 1},

p ∈ [0, 1]

Deze kan omgeschreven worden naar de vorm van vergelijking (1) gegeven in paragraaf 3.2 fp (y) = exp[log(py (1 − p)1−y )] = exp[log(py ) + log((1 − p)1−y )] = exp[y log(p) − y log(1 − p) + log(1 − p)]   p = exp[y log + log(1 − p)] 1−p   p Aangezien log 1−p = θ kan herschreven worden als p=

eθ eθ + 1

en

(1 − p) =

1 eθ + 1





(5)

(6)

geldt voor log(1 − p) als functie van θ 1 log(1 − p) = log θ e +1 = −b(θ)

De Bernoulli-verdeling behoort dus tot de exponenti¨ele familie van verdelingen waarbij voor vergelijking (1) geldt   p θ = log en b(θ) = log(1 + eθ ) 1−p en a(φ) = 1

en c(y, φ) = 0

Voor de verwachting E(Y ) = b0 (θ) geldt met behulp van vergelijking (6) ∂ log(1 + eθ ) ∂θ eθ = 1 + eθ = p

b0 (θ) =

12

(7)

De verwachting van de Bernoulli-verdeling is dus p zoals eerder al gegeven en de linkfunctie die bij deze verdeling hoort is   p g(p) = log 1−p Dit GLM is ook ge¨ımplementeerd in R en wordt gebruikt bij dit onderzoek. Voor de variantie wordt eerst b00 (θ) berekend aan de hand van de eerder eθ berekende b0 (θ) = 1+e θ en hierbij wordt gebruik gemaakt van de vergelijkingen in (6) ∂ eθ ∂θ 1 + eθ eθ eθ 1 = = · θ θ 2 θ (e + 1) 1+e e +1 = p · (1 − p)

b00 (θ) =

Invullen van de eerder gegeven a(φ) = 1 en de gevonden b00 (θ) = p(1 − p) in vergelijking (4) uit paragraaf 3.2 wordt de variantie gegeven door Var(Y ) = b00 (θ)a(φ) = p(1 − p) Dus ook de variantie van de Bernoulli-verdeling kan de herleid worden uit de algemene formule voor de variantie van de exponenti¨ele familie van verdelingen.

13

4

Analyse van gegeneraliseerde lineaire modellen met behulp van p-waarden

Om uit te vinden welke factoren relevant zijn in het sorteerproces van poststukken, is analyse nodig op de modellen zelf, maar ook tussen de modellen. Er zijn vele methoden ge¨ımplementeerd in R voor het analyseren van modellen. Een aantal van deze methoden wordt in dit hoofdstuk uitgelegd en aan de hand van een aantal testen met verschillende kleine datasets wordt gekeken naar de werking van de commando’s in R.

4.1

Methoden voor modelselectie binnen gegeneraliseerde lineaire modellen

AIC: Het Akaike’s Information Criterion (AIC) en ook het Bayesian Information Criterion (BIC) zijn lid van een algemene familie van de zogeheten ”penalized”model fittende statistieken. Het AIC is een waarde van het model die afhangt van de likelihood en het aantal modelparameters (β’s). Er geldt: ˆ + p] AIC = 2[−l(β) ˆ de likelihood van het model en p het aantal modelparaHierbij is l(β) meters dat de ”penalty”vormt op het aantal parameters. Het model met de laagste AIC-waarde is het beste. Deviantie: De deviantie van een model wordt gegeven door: ˆ ˆmax ) − l(β)]φ D = 2[l(β ˆmax ) de maximale likelihood van het volledige model, Hierbij is l(β ˆ is de likelioftewel het model waarin alle variabelen voorkomen. l(β) hood van het model waarna gekeken wordt en φ is de schalingsparameter die gegeven is in paragraaf 3.2 in vergelijking (1). In R wordt met ANalysis-Of-VAriance of analysis-of-deviance (ANOVA) modellen vergeleken aan de hand van de variantie of deviantie. p-waarden: In een model hebben de modelparameters allemaal een pwaarde. Deze ligt tussen 0 en 1 en is de kans dat in de verdeling gegeven door de nulhypothese de waarde van de toetsingsgrootheid wordt behaald of overschreden (links, rechts of tweezijdig). De pwaarde wordt door middel van een statistische toets bepaald. Bij een p-waarde van 1 kunnen we aannemen dat het gevonden resultaat op toeval berust. Een lage p-waarde betekent dat de uitkomst significant is. Bij deze p-waarden kunnen verschillende significantieniveau’s gebruikt worden, zoals 0.05 of 0.1.

14

Bij modelselectie kunnen met p-waarden, AIC en ANOVA verschillende modellen vergeleken worden. De meest gangbare manieren zijn: • Begin met een model met alle variabelen erin, bereken de p-waarden van de modelparameters en verwijder de parameters met een waarde hoger dan een bepaald significantieniveau. Bereken vervolgens van het nieuwe nieuwe model opnieuw de waarden van de parameters en verwijder parameters met hoge p-waarden. • Kies het model waarvoor de AIC minimaal is. Uiteraard zijn er nog een heleboel andere manieren voor modelselectie, maar deze komen bij dit vooronderzoek verder niet aan de orde. Na het vooronderzoek is gebleken dat bovenstaande methoden voor dit onderzoek niet de voorkeur hebben en zal gebruik worden gemaakt van BIC. Dit zal in hoofdstuk 5 uitgelegd worden.

4.2

Fitten van een logistisch regressie model in R

Voordat er in R met de gegeven dataset van TNT gemodelleerd wordt, is er eerst een klein vooronderzoek gedaan om de werking van de commando’s te begrijpen. Er zijn drie experimenten gedaan, waarbij gekeken is naar de werking van het commando glm met verschillende kleine datasets, een dataset van 100 rijen en 50 kolommen en een dataset met 1000 rijen die de originele dataset van TNT nabootst. 4.2.1

Kleine dataset

Bij dit onderzoek zijn drie modellen gebruikt die op verschillende datasets zijn toegepast. De data bestaat uit een responsvariabele y die 0 of 1 is en 2 catagorische variabelen f1 en f2. De drie modellen die gebruikt zijn, zijn Model 1: y ∼ β0 + β1 f 1 + β2 f 2 Model 2: y ∼ β0 + β1 f 1 + β2 f 2 + β3 (f 1 : f 2) Model 3: y ∼ β1 f 1 + β2 f 2 + β3 (f 1 : f 2) Hierbij is • y de responsvariabele die gelijk is aan 0 of 1 • f1 een variabele die gelijk is aan ’G’, ’B’ of ’R’ • f2 een variabele die gelijk is aan ’Hoog’ of ’Laag’ • f1:f2 de interactie tussen de variabele f1 en f2

15

In de resultaten komen een paar termen voor die hieronder verder zijn uitgelegd. Estimate: De geschatte waarde van de parameters β die horen bij de intercept of variabelen f1, f2 en f1:f2 Std Error: Standaard fout van de parameters β die horen bij de intercept of variabelen f1, f2 en f1:f2 z value: De toetsingsgrootheid z. Pr(< |z|): p-waarde van de parameter. Deze is uitgelegd in paragraaf 4.1 AIC: De AIC waarde van het model. Deze is uitgelegd in paragraaf 4.1 Met verschillende datasets zijn de modellen getoetst. Hieronder is te zien hoe de resultaten eruit zien. De waarden zijn berekend voor de eerste dataset voor model 1. Er blijkt uit de waarden dat er geen combinatie van factoren eruit spingt als meest voorkomende combinatie. In dit geval is dit ook logisch aangezien de dataset in tabel 13 in Appendix A.1.1 zo geconstrueerd is dat alle mogelijke combinaties even vaak voorkomen. Het script dat gebruikt is in R, de datasets die zijn gebruikt en de resultaten staan in Appendix A.1. Waarden: zie tabel 1 AIC: 24.636 Model 1 Intercept f1G f1R f2Laag

Estimate 0 0 0 0

Std Error 1.155 1.414 1.414 1.155

z value 0 0 0 0

Tabel 1: Resulaten test 1 van model 1

16

Pr(> |z|) 1 1 1 1

4.2.2

Freedman Experiment

Een mogelijk nadeel in modelselectie met p-waarden kan de onjuiste significantie van variabelen zijn. In deze paragraaf wordt het Freedman experiment nagebootst. Het experiment bestaat uit een dataset met 100 rijen (datapunten) en 51 kolommen (variabelen). De eerste 50 rijen zijn de onafhankelijke variabelen X1 , ..., Xn en de 51e rij is de afhankelijke responsvariable Y . Alle waarden zijn uit een Normaal-verdeling getrokken en zijn dus onafhankelijke observaties. Er zou moeten blijken dat er geen verband is tussen Y en de X-en. Ook zou uit de standaard t-test moeten blijken dat de regressie co¨effici¨enten niet significant zijn. Het script dat gebruikt is in R om dit experiment na te bootsen, staat in Appendix A.2. Eerst is de dataset gegenereert, waarna er gekeken is welke van de 50 variabelen significant zijn op twee niveaus, namelijk 0.25 en 0.05. Met significant wordt bedoeld: een p-waarde kleiner dan 0.25 dan wel 0.05. R2 is het kwadraat van het correlatie co¨effici¨ent tussen het resultaat en de voorspelde waarden van de steekproef en P is de p-waarde van het volledige model. De volgende resultaten zijn te zien: • R2 = 0.4362 en P = 0.8337 • 9 van de 50 variabelen zijn significant op 0.25 niveau • 3 van de 50 variabelen zijn significant op 0.05 niveau De 41 variabelen die een p-waarde hoger dan 0.25 hadden zijn uit het model gehaald en met het nieuwe model met 9 variabelen is weer hetzelfde experiment gedaan. De resultaten zijn nu: • R2 = 0.1976 en P = 0.0148 • 8 van de 9 variabelen zijn significant op 0.25 niveau • 3 van de 9 variabelen zijn significant op 0.05 niveau Dit experiment is totaal tien keer uitgevoerd, waarbij telkens verschillende datasets zijn gegenereerd en gebruikt. De resultaten van alle tien de experimenten worden weergegeven in tabel 2. Wat blijkt uit dit experiment, vooral uit het tweede gedeelte, is dat er wel degelijk variabelen significant zijn, terwijl in feite dit niet te verwachten is. Of in de woorden van Freedman(1983): The results from the second pass are misleading indeed, for they appear to demonstrate a definite relationship between Y and the X’s, that is, between noise and noise. Het uitgevoerde experiment illustreert dus een nadeel van modelselectie met behulp van p-waarden. 17

Nr 1 2 3 4 5 6 7 8 9 10

1e run met 50 variabelen R2 P ≤ 0.25 0.4362 0.8337 9 0.4938 0.5632 12 0.6366 0.03020 19 0.4582 0.7447 7 0.6248 0.04431 19 0.5 0.5288 14 0.4506 0.7779 9 0.4445 0.8026 12 0.4667 0.7047 20 0.5123 0.4599 14

≤ 0.05 3 0 10 0 6 5 0 1 2 5

2e run met significante variabelen R2 P ≤ 0.25 ≤ 0.05 0.1976 0.01477 8 3 0.285 0.002115 11 4 0.3313 0.01238 12 5 0.0.227 0.001002 7 3 0.5028 0.000002 17 11 0.344 0.0004788 12 10 0.1987 0.01724 6 2 0.2707 0.003971 10 4 0.3416 0.01318 18 12 0.3147 0.001872 12 7

Tabel 2: Resultaten experiment Freedman

18

4.2.3

Nabootsing dataset TNT

Eenzelfde experiment als het Freedman experiment is gedaan voor de dataset van TNT. In dit experiment wordt de dataset van TNT nagebootst met 1000 rijen (datapunten) en 14 kolommen (variabelen). De eerste kolom heeft waarden van 1 en 0 en stelt de y-variabele voor uit de TNT-dataset. De andere dertien kolommen bestaan uit catagorische variabelen, waarbij het aantal variabelen per kolom verschilt van twee tot tien. Aangezien de dataset uit 1000 variabelen bestaat en het onoverzichtelijk wordt om alles weer te geven, zijn alleen de eerste tien rijen weergegeven in tabel 3. Rij 1 2 3 4 5 6 7 8 9 10

y 0 1 1 0 0 0 0 1 0 0

x1 A B I E F F F J F C

x2 A I G F J A G J I E

x3 C C B C B B C B A A

x4 F A F E E B F A F D

x5 E E A D B F A F E A

x6 B B A D D C A A C A

x7 D A C A D A B D B D

x8 B A B B A A B A A B

x9 B A A A B B B B A B

x10 A C A B B B C C B C

x11 A C B A A B B B C B

x12 B A A B B C A A C C

x13 C A A C C B B A B B

Tabel 3: Dataset nabootsing van de TNT dataset Het model waarmee het experiment begint bevat alle variabelen en ziet er dus als volgt uit: y ∼ β0 + β1 x1 + β2 x2 + β3 x3 + β4 x4 + β5 x5 + β6 x6 + β7 x7 +β8 x8 + β9 x9 + β10 x10 + β11 x11 + β12 x12 + β13 x13

(8)

Waarbij • y de responsvariabele is die gelijk is aan 0 of 1 • x1 , x2 , ... , x13 verklarende variabelen zijn die een waarde kunnen hebben van ’A’, ’B’, ’C’, ’D’, ’E’, ’F’, ’G’, ’H’, ’I’ of ’J’. Dit hangt af van de hoeveelheid catagorische variabelen binnen x1 , x2 , ... , x13 . Van dit model (8) worden de p-waarden van de verklarende variabelen berekend. Deze zijn te zien in tabel 4. In deze tabel zijn 6 kolommen te zien, waarbij alleen gekeken wordt naar de laatste kolom, de kolom met de p-waarden. Alle verklarende variabelen met p-waarden die groter zijn dan 0.25 worden uit model (8) gehaald en het nieuwe model wordt hierdoor nu y ∼ β0 + β3 x3 + β8 x8 + β12 x12 19

(9)

NULL x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13

Df

Deviance

9 9 2 5 5 3 3 1 1 2 2 2 2

7.9482 10.9476 2.8299 3.8096 1.1716 1.1223 1.3323 2.9000 0.1285 2.6367 2.6918 4.4992 1.5601

Resid. Df 999 990 981 979 974 969 966 963 962 961 959 957 955 953

Resid. Dev 1385.0 1377.0 1366.1 1363.3 1359.5 1358.3 1357.2 1355.8 1352.9 1352.8 1350.2 1347.5 1343.0 1341.4

P(> |Chi|) 0.53938 0.27932 0.24294 0.57714 0.94757 0.77170 0.72148 0.08858 0.71998 0.26758 0.26031 0.10544 0.45839

Tabel 4: Resultaten p-waarden van de verklarende variabelen van het volledige model

x3 x8 x12

Pr(> |Chi|) 0.23194 0.08729 0.13056

Tabel 5: Resultaten p-waarden nieuwe model (9) Opnieuw worden de p-waarden van de verklarende variabelen berekend. De waarden staan in tabel 5. Er zijn dus 3 variabelen significant op 0.25 niveau, alleen x8 is significant op 0.1 niveau en geen variabelen zijn significant op 0,05 niveau. Dit experiment is net zoals het Freedman experiment een aantal keer herhaald om een beter beeld te krijgen van deze methode van modelselectie. Hierbij is telkens een nieuwe dataset gegenereerd. De verschillende datasets die zijn gebruikt worden niet meer weergegeven, maar lijken op de dataset gegeven in tabel 3. Net zoals eerder, wordt er begonnen met model (8) en vervolgens alle p-waarden 0.25 verwijderd uit het model. De resultaten staan in tabel 6. Het script dat gebruikt is in R bij dit TNT nabootsingsexperiment is te vinden in Appendix A.3.

20

Nr 1 2 3 4 5 6 7 8 9 10

1e run met 13 variabelen ≤ 0.25 ≤ 0.05 Min. p-waarde 3 0 x8 : 0.08858 5 1 x1 : 0.01638 3 1 x6 : 0.01699 3 0 x10 : 0.06246 2 1 x10 : 0.01963 2 1 x11 : 0.04952 7 1 x9 : 0.03928 6 0 x3 : 0.1035 3 1 x2 : 0.03008 4 2 x1 : 0.005572

2e run met significante variabelen ≤ 0.25 ≤ 0.05 Min. p-waarde 3 0 x8 : 0.08729 5 1 x1 : 0.01638 2 1 x6 : 0.01989 2 0 x10 : 0.09188 2 1 x10 : 0.03176 2 0 x11 : 0.06026 6 1 x9 : 0.03527 6 0 x3 : 0.1035 3 1 x2 : 0.03397 4 2 x1 : 0.005572

Tabel 6: Resultaten TNT nabootsingsexperiment

21

4.3

Conclusie vooronderzoek

In dit hoofdstuk zijn er verschillende manieren van modelselectie gebruikt om te kijken hoe de methoden werken in het programma R en welke methoden geschikt zijn voor de dataset van TNT. In het eerste experiment is gekeken naar de AIC-waarden, deviantie en p-waarden voor drie kleine datasets. Aan de hand van de AIC-waarde is het beste model gekozen. Aangezien de dataset telkens weinig data bevat, is het makkelijk de verkregen resultaten te verklaren, maar zeggen deze resultaten niet zoveel als bij grote datasets. Voor kleine datasets maakt het dus niet uit welke methode van modelselectie wordt gebruikt. In het tweede experiment is een dataset gegenereerd met vijftig onafhankelijke variabelen getrokken uit een standaard Normale verdeling. Met deze dataset is door middel van p-waarden gekeken welke significant zijn op 0.25 en 0.05 niveau. Uit dit onderzoek bleek dat er wel degelijk een verband was tussen de variabele Y en de X-en, terwijl de variabelen allemaal random gegenereerd zijn. Dus voor grotere datasets is het verkrijgen van het beste model door middel van het vergelijken van p-waarden niet de geschiktste methode. In het derde experiment is de TNT dataset nagebootst en is op dezelfde manier als het Freedman experiment modelselectie gedaan Van het volledige model worden telkens de p-waarden berekend en vervolgens worden de verklarende variabelen met p-waarden hoger dan 0.25 uit het model verwijderd. Van het nieuwe model worden weer de p-waarden berekend. Zoals te zien is in de resultaten, is er verschil tussen de aantallen significante variabelen per experiment. Aangezien er geen rekening wordt gehouden met de verschillende levels van een verklarende variabele, heeft ook deze methode niet de voorkeur. Daarom is er voor gekozen om een andere methode te bekijken, namelijk de Bayesiaanse aanpak. Deze wordt in het volgende hoofdstuk verder uitgelegd en er wordt daar ook een voorbeeld gegeven hoe deze theorie werkt en hoe met behulp van BIC modelselectie kan worden gedaan.

22

5

Bayesiaanse theorie

In de Bayesiaanse aanpak worden de onbekende parameters gezien als random variabelen. Voordat de theorie van Bayes verder wordt uitgelegd, zullen eerst een aantal veel gebruikte termen worden toegelicht.

5.1

Termen en Bayesiaanse notatie

In dit hoofdstuk is steeds X de waarneming binnen het model en θ de parametervector van het model. Prior: De prior of a priori kansverdeling is de kansverdeling van een onbekende parameter die de kennis of geloof van deze onbekende a priori, dus voordat er data is gegeven, weergeeft. In wiskundige notatie p(θ). Likelihood: De likelihood van het model, die eerder ook al bekeken is, wordt gegeven als l(θ|X) = p(X|θ) Posterior: De posterior of posterior kansverdeling is de conditionele verdeling van de onbekende parameter gegeven de data, in wiskundige notatie p(θ|X). Normaal is de notatie voor de prior, likelihood en posterior, respectievelijk PΘ (θ), PΘ|X (X|θ) en PX|Θ (θ|X). In de Bayesiaanse notatie is dit p(θ), p(X|θ) en p(θ|X).

5.2

Bayes Theorema

Het Theorema van Bayes wordt gegeven door p(θ|X) =

p(X|θ)p(θ) p(X)

ook wel bekend als p(θ|X) ∝ p(X|θ)p(θ)

(10)

waarbij ∝ het proportie teken is. In woorden: De posterior verdeling is in grootte evenredig aan de likelihood maal de prior. Hierbij is de evenredigheidsconstante 1 1 =R p(X) p(X|θ)p(θ)dθ als X continu verdeeld is en 1 1 =P p(X) p(X|θ)p(θ) als X discreet verdeeld is.

23

Stel dat er twee modellen zijn, M1 en M2 met respectievelijk parametervectors θ1 en θ2 . De vraag is nu: Gegeven de data X, welk model heeft de grootste kans de data te fitten. Om hypotheses te testen wordt gebruik gemaakt van de posterior kansen. Volgens het Theorema van Bayes is de posterior kans dat M1 het juiste model is, gelijk aan: p(M1 |X) =

p(X|M1 )p(M1 ) p(X|M1 )p(M1 ) + p(X|M2 )p(M2 )

(11)

Hierbij is Z p(X|θ1 , M1 )p(θ1 |M1 )dθ

p(X|M1 ) =

(12)

hierna ook wel de ge¨ıntegreerde likelihood voor model M1 genoemd. Eenzelfde vergelijking geldt voor M2 : p(M2 |X) =

p(X|M2 )p(M2 ) p(X|M1 )p(M1 ) + p(X|M2 )p(M2 )

met

(13)

Z p(X|M2 ) =

p(X|θ2 , M2 )p(θ2 |M2 )dθ

(14)

Het is eenvoudig te zien dat p(M1 |X) + p(M2 |X) = 1. Om nu te bepalen welk model de data X het beste fit, wordt de posterior odds van M2 tegen M1 berekend, ook wel de ratio van de posterior kansen genoemd. Deze vergelijking is met behulp van vergelijkingen (11) en (13) bepaald: p(M2 |X) p(X|M2 ) p(M2 ) = · p(M1 |X) p(X|M1 ) p(M1 )

(15)

De eerste term aan de rechterkant van de vergelijking (15) is dus de ratio tussen de ge¨ıntegreerde likelihoods van de twee modellen, de Bayes factor B21 voor M2 tegen M1 genoemd. De tweede term aan de rechterkant van de vergelijking is de ratio van priors op de modellen en deze wordt vaak 1 genomen. Dit correspondeert met p(M1 ) = p(M2 ) = 21 . In woorden kan deze vergelijking geschreven worden als Posterior odds = Bayes factor · Prior odds Als B21 > 1 dan is M2 beter, voor B21 < 1 wordt de data beter gefit door M1 .

5.3

Bayesian Information Criterion

Aangezien de ge¨ıntegreerde likelihood van een model vaak voor grote en moeilijke berekeningen zorgt, is er voor het berekenen van de Bayes factor B21 de Bayesian Information Criterion (BIC) benadering gebruikt. Deze zal

24

hieronder worden afgeleid. Voor het gemak laten we M1 weg in vergelijking (12) en schrijven we θ in paats van θ1 zodat Z p(X) = p(X|θ)p(θ)dθ Z = eg(θ) dθ met g(θ) = loge p(X|θ)p(θ) = loge p(X|θ) + loge p(θ)

(16)

De setting is als volgt: X = (X1 , ..., Xn ) met Xi onafhankelijk en identiek verdeeld. De verdeling van Xi hangt af van een parameter θ ∈ Rd . Nu wordt ˜ waarbij θ ˜ = argmax g(θ) de g(θ) in een Taylor-reeks ontwikkeld rond θ, θ 0 posterior mode is. Hierbij is g de vector met parti¨ele afgeleiden en g 00 de Hessiaan. ˜ + (θ − θ) ˜ T g 0 (θ) ˜ + 1 (θ − θ) ˜ T g 00 (θ)(θ ˜ ˜ + o(kθ − θk ˜ 2) g(θ) = g(θ) − θ) 2 ˜ + 1 (θ − θ) ˜ T g 00 (θ)(θ ˜ ˜ ≈ g(θ) − θ) 2 ˜ het maximum is van g(θ) en de afgeleide hiervan in θ ˜ dus aangezien θ ˜ = 0. Ook wordt de fout o(kθ − θk ˜ 2 ) hierna gelijk is aan nul wordt g 0 (θ) verwaarloosd. Dus Z

eg(θ) dθ Z 1 ˜ ˜ T 00 ˜ ˜ g(θ) ≈ e e 2 (θ−θ) g (θ)(θ−θ) dθ

p(X) =

˜

d

1

= eg(θ) (2π) 2 |A|− 2 ˜ De laatwaarbij d het aantal parameters is in het model en A = −g 00 (θ). ste vergelijking wordt verkregen doordat de integraal de vorm heeft van de dichtheid van een multivariabele normaal verdeling. De bovenstaande manier van benaderen heet Laplace benadering. Nu geldt 1 d loge 2π − loge |A| 2 2 ˜ + log p(θ) ˜ + d log 2π − 1 log |A| = loge p(X|θ) e e e 2 2 ˜ ≈ θ, ˆ waarbij θ ˆ de maximum likelihood schatter Voor grote datasets geldt θ 1 −1 d ˆ ˜ = A. Dit betekent is. Stel dat θ ∈ R , dan θ ' Nd (θ, i ) en niθ ≈ −g 00 (θ) ˜ + loge p(X) ≈ g(θ)

n θ

dat |A| = det(niθ ) ≈ nd |iθ | 25

met i de verwachte Fisher informatie voor ´e´en waarneming. Nu geldt ˆ + log p(θ) ˆ + d log 2π − d log n − 1 log |i| (17) loge p(X) ≈ loge p(X|θ) e e e e 2 2 2 De eerste term aan de rechterkant van vergelijking (17) is van orde O(n) en de vierde term van orde O(loge n). De rest van de termen zijn van orde O(1) of kleiner en worden verwaarloosd. De vergelijking die overblijft is ˆ − loge p(X) ≈ loge p(X|θ)

d loge n 2

Met de benadering van (18) kan de Bayesfactor B21 = naderd, waarbij dus θ 1 ∈

Rd1

en θ 2 ∈

(18) p(X|M2 ) p(X|M1 )

worden be-

Rd2 .

p(X|M2 ) p(X|M1 ) = 2 loge p(X|M2 ) − 2 loge p(X|M1 ) ˆ 2 , M2 ) − d2 log n] − 2[log p(X|θ ˆ 1 , M1 ) − d1 log n] ≈ 2[loge p(X|θ e e e 2 2 ˆ 2 , M2 ) p(X|θ − (d2 − d1 ) loge n = 2 loge ˆ 1 , M1 ) p(X|θ

2 loge B21 = 2 loge

= −BIC21 De BIC is dus BIC21 ≈ −2 loge

ˆ 2 , M2 ) p(X|θ − (d2 − d1 ) loge n ˆ 1 , M1 ) p(X|θ

(19)

en de verhouding tussen BIC21 en B21 is gelijk aan B21 ≈ e− 2 BIC21 1

5.4

(20)

Voorbeeld

Om de theorie in paragraaf 5.2 en 5.3 beter te begrijpen en deze te koppelen aan de dataset van TNT, is er in R een klein onderzoek gedaan aan de hand van een voorbeeld. De volgende gegevens zijn nodig in dit voorbeeld: • Random datapunten X = (X1 , ..., Xn ) • Model 1 (M1 ): Xi ∼ P oisson(λ), a priori verdeling op λ is Exp(1), λ>0 • Model 2 (M2 ): Xi ∼ Bin(10, θ), a priori verdeling op θ is Beta(α, β), θ ∈ [0, 1]

26

In R worden verschillende datasets gegenereerd, waarbij n datapunten getrokken worden uit een P ois(λ)-verdeling met λ = 0.5, 1, 2 en n = 10, 100, 1000. Van deze datasets wordt de Bayesfactor B21 uitgerekend en dit wordt 1000 keer herhaald. Van deze 1000 B21 -waarden wordt een histogram gemaakt. Het script dat gebruikt is in R, is terug te vinden in Appendix A.4. Aangezien R moeite heeft met grote getallen die in de Gamma- en Beta-verdeling kunnen voorkomen, is gebruik gemaakt van log(B21 ) en dus de log(Beta) en log(Gamma). Uitendelijk is weer de e-macht van deze log(B21 ) genomen om de juiste waarde voor B21 te vinden. Om de 2 modellen te vergelijken en te bepalen welk model beter de data fit, wordt er gebruikt gemaakt van vergelijking (15) uit paragraaf 5.2. Hierbij wordt aangenomen dat P (M1 ) = P (M2 ) = 12 en dus geldt p(M2 |X) p(X|M2 ) = = B21 p(M1 |X) p(X|M1 ) waarbij B21 geschreven kan worden als R1 0 p(X|M2 , θ)p(θ)dθ B21 = R ∞ 0 p(X|M1 , λ)p(λ)dλ met p(X|M1 , λ) =

n Y i=1

en p(X|M2 , θ) =

n   Y 10

xi

i=1

e−λ

(21)

λ xi xi !

θxi (1 − θ)10−xi

Nu worden de integralen berekend, waarbij λ en θ een verdeling hebben van respectievelijk Exp(1) en Beta(α, β) zoals eerder gegeven, dus de kansdichtheid p(λ) en p(θ) worden hierdoor p(λ) = e−λ , λ ≥ 0 Γ(α + β) α−1 p(θ) = θ (1 − θ)β−1 , Γ(α)Γ(β)

θ ∈ [0, 1]

De noemer (21) wordt nu: Z

∞

∞

λΣxi −λ e−λn Qn e dλ 0 i=1 xi ! Z ∞ = c(x) e−λ(n+1) λΣxi dλ Z

p(X|M1 , λ)p(λ)dλ = 0

0

met

27

c(x) = Qn

1

i=1 xi !

(22)

= =

∞

tΣxi 1 dt Σxi n + 1 (n + 1) 0 Z ∞ c(x) e−t tΣxi +1−1 dt (n + 1)Σxi +1 0 1 Q Γ(Σxi + 1) Σx i (n + 1) +1 ni=1 xi ! Z

= c(x)

e−t

(23)

De derde regel van bovenstaande vergelijking wordt verkregen door substit tutie van t = λ(n + 1) → λ = n+1 . De laatste regel wordt verkregen door de definitie van de Gamma-functie, Z ∞ Γ(z) = e−t tz−1 dt met <(z) > 0 0

De teller van vergelijking (21) kan herschreven worden naar: Z 1 Z 1Y n   10 xi Γ(α + β) α−1 θ (1 − θ)β−1 dθ p(X|M2 , θ)p(θ)dθ = θ (1 − θ)10−xi Γ(α)Γ(β) 0 0 i=1 xi " n  # Z Y 10 Γ(α + β) 1 Σxi +α−1 = θ (1 − θ)10n−Σxi +β−1 dθ xi Γ(α)Γ(β) 0 i=1 " n  # Y 10 Γ(α + β) Γ(Σxi + α)Γ(10n − Σxi + β) · = Γ(α)Γ(β) Γ(Σxi + α + 10n − Σxi + β) xi i=1 Z 1 Γ(Σxi + α + 10n − Σxi + β) Σxi +α−1 θ (1 − θ)10n−Σxi +β−1 dθ · Γ(Σx α)Γ(10n − Σx + β) i i 0 Aangezien de Beta-functie geschreven kan worden als Beta(α, β) =

Γ(α)Γ(β) Γ(α + β)

is de term in de integraal de kansdichtheid van de Beta(Σxi +α, 10n−Σxi +β) verdeling en is de integraal 1. Hierdoor kan de laatste regel vereenvoudigd worden tot " n  # Z 1 Y 10 Beta(Σxi + α, 10n − Σxi + β) p(X|M2 , θ)p(θ)dθ = (24) xi Beta(α, β) 0 i=1

Vervolgens kan met de verkregen vergelijkingen (23) en (24) de Bayesfactor B21 berekend worden. R1 0 p(X|M2 , θ)p(θ)dθ B21 = R ∞ p(X|M1 , λ)p(λ)dλ "0n  # Q Y 10 Beta(Σxi + α, 10n − Σxi + β) (n + 1)Σxi +1 ni=1 xi ! = · xi Beta(α, β) Γ(Σxi + 1) i=1 hQ i n 10! Σxi +1 i=1 (10−xi )! (n + 1) Beta(Σxi + α, 10n − Σxi + β) = · (25) Γ(Σxi + 1) Beta(α, β) 28

Waarbij gebruik wordt gemaakt van n   Y 10 i=1

xi

=

n Y i=1

10! (10 − xi )!xi !

Nu kan er met behulp van vergelijking (25) gekeken worden welke model het beste de data fit. Als B21 < 1 dan is model 1 beter, bij B21 > 1 is model 2 beter. Er wordt nu gekeken in hoeverre dit klopt voor dit voorbeeld. De 9 histogrammen zijn te zien in figuur 2. In tabel 7 is te zien voor de verschillende n en λ wat de fractie van B21 < 1. Zowel in de histogrammen als in de tabel is af te lezen dat voor bijna elke combinatie van n = 10, 100, 1000 en λ = 0.5, 1, 2 het grootste gedeelte van de waarden B21 < 1, wat betekent dat model 1 de data het beste fit. Voor de P ois(2)verdeling met n = 10 en n = 100 neigt de data naar model 1, maar niet heel duidelijk. Dit komt door de kleine steekproefgrootte en dus door een sterke invloed van de prior.

n = 10 n = 100 n = 1000

λ = 0.5 1 1 1

λ=1 0.998 1 0.971

λ=2 0.756 0.854 0.998

Tabel 7: Ratio B21 -waarden < 1

29

Figuur 2: Histogrammen voor de B21 -waarden van P ois(λ), met n = 10, 100, 1000 en λ = 0.5, 1, 2

30

Nu wordt met behulp van de vergelijking (19) uit de vorige paragraaf de BIC berekend voor deze twee modellen. Voor model 1 is de parameter λ en dus is d1 = 1 en geldt ˆ = p(X|M1 , λ)

n Y i=1

ˆ ˆλ

e −λ

xi

xi !

ˆ=X ¯ n de maximum likelihood schatter (MLE). Voor model 2 is Hierbij is λ de parameter θ en dus is d2 = 1 en geldt n   Y 10 ˆxi ˆ ˆ 10−xi p(X|M2 , θ) = θ (1 − θ) xi i=1

¯ n /10 de MLE. Voor de berekening van de twee MLE’s, zie Hierbij is θˆ = X Appendix A.4.1. Invullen in vergelijking (19) wordt ˆ M2 ) p(X|θ, − (d2 − d1 ) loge n ˆ M1 ) p(X|λ, Qn 10
ˆxi ˆ 10−xi i=1 xi θ (1 − θ) −2 loge − (1 − 1) loge (n) Qn −λˆ λˆ xi i=1 e xi ! Qn 10
X¯ n xi ¯ n 10−xi X i=1 xi ( 10 ) (1 − 10 ) −2 loge Qn −X¯ X¯ nxi n i=1 e xi !   ¯ nxi 10−X ¯ ¯ n −xi X X 10! n Y (10−xi )!xi ! 10xi ( 10 n )10 ( 10− ) 10   −2 loge ¯ n ¯ xi 1 − X e X n xi ! i=1   10 (10−X ¯ ¯ n )−xi (10− X ) n 10! 1 n Y x −x 10 i 10 10 i  (10−xi )! 10  −2 loge ¯n −X e i=1 n  Y ¯ n )10−xi ¯  (10 − X 10! −2 loge eXn (10 − xi )! 1010

BIC21 ≈ −2 loge =

=

=

=

=

(26)

i=1

Aangezien de verhouding tussen B21 en BIC21 gegeven wordt door vergelijBIC king (20), geldt voor B21 BIC ≈ e− 2 BIC21 B21 1

− 21 −2 loge

= e n  Y = i=1

Qn

i=1



¯ n )10−xi (10−X 10! (10−xi )! 1010

¯

e Xn

¯ n )10−xi ¯ 10! (10 − X eXn (10 − xi )! 1010



 (27)

Met deze vergelijking (27) wordt nu in R voor verschillende waarden van BIC uitgerekend. De resultaten zijn weergegeven in de vorm van n en λ de B21 31

BIC < 1 in tabel 8. Het script histogrammen in figuur 3 en de ratio van B21 dat gebruikt is, is te vinden in Appendix A.4. Zoals te zien is in de tabel komt model 1 voor n = 10 en n = 100 niet als best fittende model naar voren. Dit komt door de kleine steekproefgrootte en dus door een sterke invloed van de prior. Voor n = 1000 komt wel model 1 als best fittende model naar voren.

BIC -waarden van P ois(λ), met n = Figuur 3: Histogrammen voor de B21 10, 100, 1000 en λ = 0.5, 1, 2 In dit voorbeeld is de gegenereerde dataset zo gekozen, dat er uit de resultaten naar voren moest komen dat model 1 het best fittende model was. Zoals in de resultaten te zien is komt dit ook naar voren in de verschillende grafieken. Wel moet gezegd worden dat hierbij een grote dataset het beste 32

n = 10 n = 100 n = 1000

λ = 0.5 0.302 0.49 0.713

λ=1 0.383 0.595 0.888

λ=2 0.433 0.745 0.997

BIC -waarden < 1 Tabel 8: Ratio B21 resultaat zal bieden bij modelselectie via BIC. Dit blijkt ook uit ons voorbeeld, aangezien slechts voor de ratio’s met n = 1000 voor het grootste deel BIC < 1. de B21

33

6

Onderzoek dataset TNT

In dit hoofdstuk zal met behulp van de package Bayesian Model Averaging (BMA) analyse gedaan worden op de dataset van TNT. In deze analyse wordt gebruik maakt van het Bayesiaans hypothese testen met behulp van het Bayesian Information Criterion (BIC). Het script dat hiervoor gebruikt is, staat in Appendix B.

6.1

Analyse in Bayesian Model Averaging

Bayesian Model Averaging (BMA) is een manier van hypothese testen waarbij rekening wordt gehouden met de onzekerheid van een model. De BMA posterior verdeling is een gewogen gemiddelde van zijn posterior verdelingen van de modellen die beschouwd worden. Hierbij is het gewicht van een model gelijk aan de posterior kans dat het model correct is, gegeven dat 1 van de modellen die wordt beschouwd correct is. Als de geobserveerde data wordt gemodelleerd door gegeneraliseerde lineaire modellen kan de model onzekerheid groot zijn. Er zijn dus twee moeilijke problemen: Het berekenen van de ge¨ıntegreerde likelihoods van alle modellen en het middelen over alle modellen, waarbij het aantal modellen heel groot kan zijn. De R package BMA gebruikt de BIC benadering om de ge¨ıntegreerde likelihoods te berekenen. Verder wordt gebruik gemaakt van een ’leaps and bound’ algoritme en worden a posteriori de modellen die het minst waarschijnlijk zijn eruit gehaald. Dit is een exhaustive search (uitgebreid onderzoek waarbij alle mogelijkheden worden bekeken) en deze vindt uiteindelijk het globaal optimale model.

6.2

Opzet onderzoek

Zoals eerder is vermeld, bestaat de dataset uit een aantal verklarende variabelen en ongeveer 20000 datapunten. De variabelen die gebruikt zijn voor het onderzoek zijn besproken in paragraaf 2.2. Deze variabelen komen allemaal voor in het eerste model, namelijk: Optijd ∼ β0 + β1 VerzendenDvdW + β2 Gewicht + β3 Formaat +β4 Schrift + β5 InterofIntra + β6 Streepjescode +β7 Poststroom + β8 InfoPRIC

(28)

De 8 verklarende variabelen in model (28) hebben catagorische subvariabelen, wat betekent dat ze per variabele meerdere mogelijkheden hebben. Voor de variable VerzendenDvdW heeft bijvoorbeeld de mogelijkheden Rest, Vrijdag of Zondag met respectievelijk de β’s β1rest , β1vrijdag en β1zondag . Dit geldt ook voor de overige variabelen. In R wordt deze codering met behulp van dummy variabelen gedaan.

34

6.3

Resultaat

Het resultaat voor de analyse op model 1 is te vinden in tabel 9 en figuren 4 en 5. Hierbij geeft tabel 9 een samenvatting van welke factoren wel of niet invloed hebben op het op tijd aankomen van de post. Hierbij is de kolom met ’p! = 0’ de posterior kans dat de variabele in het model zit in %, ’EV’ is de BMA posterior mean voor elke variabele en ’SD’ de BMA posterior standaard deviatie voor elke variabele. De kolommen daarna bevatten de parameter schatters per model per variabele. Onderaan in de tabel staan bij ’nVar’ het aantal variabelen in het model, bij ’BIC’ de BIC-waarden per model en bij ’post prob’ de bijbehorende posterior kansen per model. In figuur 4 zijn de posterior kansen van de verschillende variabelen te zien. De verticale lijn bij 0 is de posterior kans dat de variabele niet in het model zit. De kromme is de gewogen model posterior kans van de co¨effici¨ent, gegeven dat de variabele in het model zit. In figuur 5 is te zien welke verklarende variabelen invloed hebben op het op tijd verzenden van de post. Blauw betekent een negatieve parameter schatter van de variabele en rood een positieve parameter schatter. Zoals te zien is in de tabel en de figuren, komen er twee modellen uit de analyse. Model 1 met variabelen VerzendenDvdW, Streepjescode en InfoPric met een posterior kans van 0.741 en model 2 met de VerzendenDvdW, InterofIntra, Streepjescode en InfoPric met een posterior kans van 0.259. Dit betekent dat de variabelen VerzendenDvdW, InterofIntra, Streepjescode en InfoPric 100% van invloed zijn op het overkomstpercentage van de post en de variabele InterofIntra 25.9%.

35

p!=0

EV

SD

model 1

model 2

100

3.555965

0.06501

3.577e+00

3.497e+00

.Vrijdag

-0.006104

0.09332

-6.023e-03

-6.335e-03

.Zondag

-0.691765

0.09312

-6.916e-01

-6.923e-01

0.000000

0.00000

.

.

.C5

0.000000

0.00000

.

.

.C6

0.000000

0.00000

.

.

0.000000

0.00000

.

.

0.000000

0.00000

.

.

0.049605

0.09098

.

1.918e-01

-0.485287

0.10213

-4.848e-01

-4.867e-01

.Groot

0.000000

0.00000

.

.

.Klein

0.000000

0.00000

.

.

-1.448984

0.18386

-1.448e+00

-1.453e+00

Intercept VerzendenDvdW

100.0

Gewicht

0.0

Formaat

0.0

.EA5 Schrift

0.0 .Machine

InterofIntra

25.9 .Intra

Streepjescode

100.0 .nee

Poststroom

0.0

InfoPRIC

100.0 .aanwonjuist .Geen

-0.793888

0.09478

-7.933e-01

-7.957e-01

.nietgemeten

-19.583226

213.01181

-1.958e+01

-1.959e+01

nVar

3

4

BIC

-2.035e+05

-2.035e+05

0.741

0.259

post prob

Tabel 9: Resultaten na Bayesiaanse analyse met InfoPRIC

36

Figuur 4: Resultaat posterior kansen per variabele met InfoPRIC

Figuur 5: Resultaat modellen geselecteerd door BMA met InfoPRIC 37

Naar aanleiding van dit resultaat is een tweede model geanalyseerd, zonder de variabele InfoPRIC. Deze is uit het model gehaald aangezien het bleek dat er van veel datapunten geen Pric-code was of de Pric-code niet gemeten (NA) is. De resultaten zijn weergegeven in tabel 10 en figuren 6 en 7. Aangezien tabel 10 vrij groot is, is deze voor de duidelijkheid geroteerd afgebeeld.

Figuur 6: Resultaat posterior kansen per variabele zonder InfoPRIC

Figuur 7: Resultaat modellen geselecteerd door BMA zonder InfoPRIC

38

39

nVar BIC post prob

Poststroom

Streepjescode

InterofIntra

Schrift

Gewicht Formaat

Intercept VerzendenDvdW

.Groot .Klein

.nee

.Intra

.Machine

.C5 .C6 .EA5

.Vrijdag .Zondag

-0.15389 -0.10177

-1.44560

0.01990

0.03955

0.00000 0.00000 0.00000

-0.04464 -0.68747 0.00000

EV 3.47345

0.26139 0.18122

0.07091

0.05838

0.08035

0.00000 0.00000 0.00000

0.08670 0.08913 0.00000

SD 0.19999

. . 2 -2.029e+05 0.477

-1.436e+00

.

.

. . .

-4.458e-02 -6.902e-01 .

model 1 3.392e+00

-5.631e-01 -3.727e-01 3 -2.029e+05 0.182

-1.472e+00

.

.

. . .

-4.482e-02 -6.819e-01 .

model 2 3.788e+00

. . 3 -2.029e+05 0.160

-1.433e+00

.

1.796e-01

. . .

-6.881e-01 . .

model 3 3.310e+00

Tabel 10: Resultaten na Bayasiaanse analyse zonder InfoPRIC

27.4

100.0

12.0

22.0

0.0 0.0

p!=0 100 100.0

. . 3 -2.029e+05 0.089

-1.439e+00

1.661e-01

.

. . .

-6.901e-01 . .

model 4 3.322e+00

-5.622e-01 -3.707e-01 4 -2.029e+05 0.060

-1.469e+00

.

1.793e-01

. . .

-6.800e-01 . .

model 5 3.705e+00

In de figuren hierboven zijn de catagorische variabelen als geheel genomen en is te zien welke variabelen als geheel invloed hebben op het overkomstpercentage van de post, namelijk met een posterior kans van 100% VerzendenDvdW en Streepjescode en in mindere mate de variabelen Schrift (22%), InterofIntra (12%) en Poststroom (27.4%). In de figuren 8 en 9 is hetzelfde model gebruikt als hierboven, maar nu wordt gekeken of de aparte levels binnen een catagorische variabele in het model zit. Uit de figuren blijkt dat VerzendenDvdWZondag en Streepjescodenee in alle modellen invloed heeft en PoststroomGroot en FormaatC6 in bijna alle modellen voorkomt.

Figuur 8: Resultaat posterior kansen per variabele zonder InfoPRIC, aparte levels

40

Figuur 9: Resultaat modellen geselecteerd door BMA zonder InfoPRIC, aparte levels

41

6.4

Testen data 2008

Aan de hand van de analyse met BMA op de data uit 2007 is gebleken dat de factoren VerzendenDvdW en Streepjescode voornamelijk invloed hebben op het overkomstpercentage. Dit is ook terug te zien als model (29) op dezelfde manier wordt geanalyseerd als gedaan is in het vooronderzoek in paragraaf 4.2.1 oftewel met behulp van analyse met p-waarden. Het model dat gebruikt is Optijd ∼ β0 + β1 VerzendenDvdW + β2 Gewicht + β3 Formaat +β4 Schrift + β5 InterofIntra + β6 Streepjescode +β7 Poststroom

(29)

De resultaten die hierbij horen staan in tabel 11. Hieruit blijkt dat de

Intercept VerzendenDvdWVrijdag VerzendenDvdWZondag Gewicht FormaatC5 FormaatC6 FormaatEA5 SchriftMachine InterofIntraIntra Streepjescodenee PoststroomGroot PoststroomKlein

Estimate 3.828e+00 -4.487e-02 -6.807e-01 -5.985e-05 -2.362e-01 -5.035e-01 -1.906e-01 1.780e-01 1.638e-01 -1.478e+00 -7.227e-01 -1.464e-01

Std. Error 2.443e-01 8.681e-02 8.933e-02 1.203e-04 2.564e-01 2.728e-01 3.262e-01 6.472e-02 6.507e-02 7.636e-02 2.214e-01 1.848e-01

z value 15.668 -0.517 -7.620 -0.498 -0.921 -1.846 -0.584 2.750 2.518 -19.355 -3.263 -0.792

Pr(> |z|) < 2e-16 0.60529 2.54e-14 0.61881 0.35702 0.06494 0.55894 0.00596 0.01181 < 2e-16 0.00110 0.42816

Tabel 11: Resultaten data 2007, model (29) factoren VerzendenDvdWZondag en Streepjescodenee van grote invloed zijn op het overkomstpercentage van de post, aangezien beide een grote negatieve waarde hebben (Estimate) en de p-waarden (Pr(> |z|)) voor deze variabelen heel laag zijn. Dit is ook gebleken uit analyse met BMA. In mindere mate heeft PoststroomGroot en SchriftMachine invloed. Hetzelfde is gedaan voor de data van 2008. De resultaten hiervan staan in tabel 12. Ook uit deze tabel blijkt dat VerzendenDvdWZondag en Streepjescodenee van invloed zijn op het overkomstpercentage van de post. Het script dat gebruikt is voor deze resultaten is ook te vinden in Appendix B.

6.5

Conclusie onderzoek TNT

Als er gekeken wordt naar de factoren op zich, dan blijkt uit de tabellen en figuren in paragraaf 6.3 dat de factoren VerzendenDvdW en Streepjescode 42

Intercept VerzendenDvdWVrijdag VerzendenDvdWZondag Gewicht FormaatC5 FormaatC6 FormaatEA5 SchriftMachine InterofIntraIntra Streepjescodenee PoststroomGroot PoststroomKlein

Estimate 3.5263085 -0.1652666 -0.7607952 0.0001329 0.2997801 0.1724856 0.9135383 0.1428450 0.2567374 -1.3676346 -0.3147914 -0.4936951

Std. Error 0.2448625 0.0875812 0.0933721 0.0001323 0.2663927 0.2818344 0.3668452 0.0675161 0.0696719 0.0832544 0.2198468 0.2015441

z value 14.401 -1.887 -8.148 1.004 1.125 0.612 2.490 2.116 3.685 -16.427 -1.432 -2.450

Pr(> |z|) < 2e-16 0.059159 3.7e-16 0.315262 0.260449 0.540531 0.012765 0.034369 0.000229 < 2e-16 0.152182 0.014303

Tabel 12: Resultaten data 2008, model (29) met een posterior kans van 100% invloed hebben op het overkomstpercentage van de post. Als er gekeken wordt naar de verschillende levels van een variabele, dan komen de variabelen VerzendenDvdWZondag en Streepjescodenee beide met posterior kans 100% uit het model naar voren. Zoals te zien is in de tabel 11 en 12 geldt dit ook als er wordt gekeken naar de p-waarden van de variabelen. Dit geldt voor de data van 2007 en voor de data van 2008.

43

7

Conclusie

In dit verslag is een aantal manieren bekeken om kleine en grote datasets met catagorische variabelen te analyseren. Er is gebruikt gemaakt van gegeneraliseerde lineaire modellen, Bayesiaanse aanpak, AIC, p-waarden en BIC. Naast de theorie die is uitgelegd, zijn er ook voorbeelden gegeven waarin de theorie wordt toegepast. In deze voorbeelden is gebleken dat met een grote dataset met catagorische variabelen met de Bayesiaanse aanpak het beste resultaat verkregen kan worden. Er is dus daarom gekozen om met behulp van de Bayesian MOdel Averaging (BMA) analyse te doen op de TNT dataset. Uit dit onderzoek blijkt dat de factoren VerzendenDvdW en Streepjescode met posterior kans 100% de grootste invloed hebben op het overkomstpercentage van de poststukken, in het bijzonder de subvariabelen VerzendenDvdWZondag en Streepjescodenee. Ook blijkt dat als er geen Pric-code is gemeten of deze er niet opstaat, deze factor ook een grote invloed heeft. Dit betekent dat er door TNT extra aandacht kan worden besteed aan het aanbrengen van een Streepjescode en een Pric-code op de poststukken. Ook betekent dit dat het verzenden van de post op zondag of het ontvangen van de post op maandag minder goed verloopt dan op de andere dagen.

44

A

Script vooronderzoek, datasets en resultaten

Hieronder staat het script dat in R is uitgevoerd voor het toetsen van de drie modellen voor de verschillende kleine datasets, het Freedman experiment, de nabootsing van de TNT dataset en het Bayes voorbeeld. Ook is het hele onderzoek voor de kleine datasets beschreven met datasets en resulaten.

A.1

Kleine dataset

In deze paragraaf wordt het gehele onderzoek met de kleine datasets weergegeven. Er zijn drie datasets gegenereerd waarop telkens op dezelfde manier modelselectie wordt gedaan tussen drie modellen. De data bestaan uit een responsvariabele y die 0 of 1 is en twee catagorische variabelen f1 en f2. De drie modellen die gebruikt zijn, zijn Model 1: y ∼ β0 + β1 f 1 + β2 f 2 Model 2: y ∼ β0 + β1 f 1 + β2 f 2 + β3 (f 1 : f 2) Model 3: y ∼ β1 f 1 + β2 f 2 + β3 (f 1 : f 2) Hierbij is • y de responsvariabele die gelijk is aan 0 of 1 • f1 een variabele die gelijk is aan ’G’, ’B’ of ’R’ • f2 een variabele die gelijk is aan ’Hoog’ of ’Laag’ • f1:f2 de interactie tussen de variabele f1 en f2 In de resultaten komen een paar termen voor die hieronder verder zijn uitgelegd. Estimate: De geschatte waarde van de parameters β die horen bij de intercept of variabelen f1, f2 en f1:f2 Std Error: Standaard fout van de parameters β die horen bij de intercept of variabelen f1, f2 en f1:f2 z value: De toetsingsgrootheid z. Pr(< |z|): p-waarde van de parameter. Deze is uitgelegd in paragraaf 4.1 AIC: De AIC waarde van het model. Deze is uitgelegd in paragraaf 4.1

45

A.1.1

Test 1

De data die gebruikt is in deze test staat in tabel 13. De data is niet random gegenereerd, maar met de hand geconstrueerd. In deze test komt elke combinatie van variabelen ´e´en keer voor. Het script dat gebruikt is in R om de verschillende modellen te toetsen is: x2<-data.frame(y=rep(c(0,1),rep(6,2)), f1=rep(c("G","R","B"),rep(2,3)), f2=rep(c("Hoog","Laag"),rep(1,2))) x2.glm1<-glm(y~f1+f2,family=binomial, data=x2) summary(x2.glm1) x2.glm2<-glm(y~f1*f2,family=binomial, data=x2) summary(x2.glm2) x2.glm3<-glm(y~-1+f1:f2,family=binomial, data=x2) summary(x2.glm3) y 0 0 0 0 0 0 1 1 1 1 1 1

f1 G G R R B B G G R R B B

f2 Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag

Tabel 13: Data van test 1 Resultaten voor model 1: Waarden: zie tabel 14 AIC: 24.636 Resultaten voor model 2: Waarden: zie tabel 15 AIC: 28.636 Resultaten voor model 3:

46

Model 1 Intercept f1G f1R f2Laag

Estimate 0 0 0 0

Std Error 1.155 1.414 1.414 1.155

z value 0 0 0 0

Pr(> |z|) 1 1 1 1

Tabel 14: Resulaten test 1 van model 1 Model 2 Intercept f1G f1R f2Laag f1G:f2Laag f1R:f2Laag

Estimate 0 0 0 0 0 0

Std Error 1.414 2.000 2.000 2.000 2.828 2.828

z value 0 0 0 0 0 0

Pr(> |z|) 1 1 1 1 1 1

Tabel 15: Resulaten test 1 van model 2 Waarden: zie tabel 16 AIC: 28.636 Zoals te zien is in de resultaten van de drie modellen, blijkt uit de waarden dat er geen combinatie van factoren eruit spingt als meest voorkomende combinatie. Dit is ook logisch aangezien de dataset in tabel 13 zo geconstrueerd is dat alle mogelijke combinaties even vaak voorkomen. Uit de analyse met behulp van AIC blijkt dat model 1 het beste is voor deze test, aangezien deze de laagste AIC-waarde heeft. A.1.2

Test 2

De data die gebruik is in deze test staat in tabel 17. De data lijkt op de data van test 1, behalve dat nu blijkt uit de gegevens dat G=0, R=1 en B=1. Het script dat gebruikt is in R om de verschillende modellen te toetsen is: x3<-data.frame(y=c(0,0,1,1,1,1,0,0,1,1,1,1), f1=rep(c("G","R","B"),rep(2,3)), f2=rep(c("Hoog","Laag"),rep(1,2))) x3.glm1<-glm(y~f1+f2,family=binomial, data=x3) summary(x3.glm1) x3.glm2<-glm(y~f1*f2,family=binomial, data=x3) summary(x3.glm2) x3.glm3<-glm(y~-1+f1:f2,family=binomial, data=x3) summary(x3.glm3) Resultaten voor model 1: 47

Model 3 f1B:f2Hoog f1G:f2Hoog f1R:f2Hoog f1B:f2Laag f1G:f2Laag f1R:f2Laag

Estimate 0 0 0 0 0 0

Std Error 1.414 1.414 1.414 1.414 1.414 1.414

z value 0 0 0 0 0 0

Pr(> |z|) 1 1 1 1 1 1

Tabel 16: Resulaten test 1 van model 3 y 0 0 1 1 1 1 0 0 1 1 1 1

f1 G G R R B B G G R R B B

f2 Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag

Tabel 17: Data van test 2 Waarden: zie tabel 18 AIC: 8 Model 1 Intercept f1G f1R f2Laag

Estimate 2.457e01 -4.913e01 0 0

Std Error 7.564e04 9.264e04 9.264e04 7.564e04

z value 3.25e-04 -0.001 0 0

Tabel 18: Resulaten test 2 van model 1 Resultaten voor model 2: Waarden: zie tabel 19 AIC: 12 Resultaten voor model 3:

48

Pr(> |z|) 1 1 1 1

Model 2 Intercept f1G f1R f2Laag f1G:f2Laag f1R:f2Laag

Estimate 2.457e01 -4.913e+01 0 0 0 0

Std Error 9.264e04 1.310e+05 1.310e+05 1.310e+05 1.853e+05 1.853e+05

z value 2.65e-04 -3.75e-04 0 0 0 0

Pr(> |z|) 1 1 1 1 1 1

Tabel 19: Resulaten test 2 van model 2 Waarden: zie tabel 20 AIC: 12 Model 3 f1B:f2Hoog f1G:f2Hoog f1R:f2Hoog f1B:f2Laag f1G:f2Laag f1R:f2Laag

Estimate 24.57 -24.57 24.57 24.57 -24.57 24.57

Std Error 92638.56 92638.55 92638.56 92638.56 92638.55 92638.61

z value 0.000265 -0.000265 0.000265 0.000265 -0.000265 0.000265

Pr(> |z|) 1 1 1 1 1 1

Tabel 20: Resulaten test 2 van model 3 Nu is er wel verschil te zien in de waarden in de tabellen van de drie modellen. Duidelijk wordt dat een combinatie met f1G een negatieve invloed heeft. Deze is gelijk voor zowel f1G:f2Hoog als f1G:f2Laag. Dit komt ook terug in de dataset aangezien deze zo geconstrueerd is dat G=0, R=1 en B=1. Volgens de AIC-waarden is wederom model 1 het beste model voor deze test. A.1.3

Test 3

De data die gebruik is in deze test staat in tabel 21. Deze data is ook weer met de hand geconstrueerd, maar mist de combinatie G=Laag en R=0. Het script dat gebruikt is in R om de verschillende modellen te toetsen is: x4<-data.frame(y=c(0,1,1,1,1,1,1,1,0,1,1,0), f1=c("G","R","B","R","B","B","G","R","B","R","R","B"), f2=rep(c("Hoog","Laag"),rep(1,2))) x4.glm1<-glm(y~f1+f2,family=binomial, data=x4) summary(x4.glm1) x4.glm2<-glm(y~f1*f2,family=binomial, data=x4) summary(x4.glm2) x4.glm3<-glm(y~-1+f1:f2,family=binomial, data=x4) 49

summary(x4.glm3) x4.glm4<-glm(y~f1,family=binomial, data=x4) summary(x4.glm4) x4.glm5<-glm(y~f2,family=binomial, data=x4) summary(x4.glm5) y 0 1 1 1 1 1 1 1 0 1 1 0

f1 G R B R B B G R B R R B

f2 Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag Hoog Laag

Tabel 21: Data van test 3 De resultaten voor model 1: Waarden: zie tabel 22 AIC: 17.364 Model 1 Intercept f1G f1R f2Laag

Estimate 0.6931 -0.6931 19.4469 -0.6931

Std Error 1.2247 1.8708 4776.3212 1.8708

z value 0.566 -0.371 0.004 -0.371

Tabel 22: Resulaten test 3 van model 1 Resultaten voor model 2: Waarden: zie tabel 23 AIC: 19.364 Resultaten voor model 3: Waarden: zie tabel 24 AIC: 19.364 50

Pr(> |z|) 0.571 0.711 0.997 0.711

Model 2 Intercept f1G f1R f2Laag f1G:f2Laag f1R:f2Laag

Estimate 0.6931 -0.6931 18.8729 -0.6931 NA 0.6931

Std Error 1.2247 1.8708 10754.0130 1.8708 NA 12023.3521

z value 0.566 -0.371 0.002 -0.371 NA 5.77e-05

Pr(> |z|) 0.571 0.711 0.999 0.711 NA 1.000

Tabel 23: Resulaten test 3 van model 2 Model 3 f1B:f2Hoog f1G:f2Hoog f1R:f2Hoog f1B:f2Laag f1G:f2Laag f1R:f2Laag

Estimate 6.931e-01 0 1.957e+01 0 NA 1.957e+01

Std Error 1.225 1.414 1.075e+04 1.414 NA 5.377e+03

z value 0.566 0 0.002 0 NA 0.004

Pr(> |z|) 0.571 1.000 0.999 1.000 NA 0.997

Tabel 24: Resulaten test 3 van model 3 Na aanleiding van bovenstaande resultaten zijn er voor deze data nog 2 modellen onderzocht: Model 4: y ∼ β0 + β1 f 1 Model 5: y ∼ β0 + β1 f 2 De resultaten voor model 4: Waarden: zie tabel 25 AIC: 15.503 Model 4 Intercept f1G f1R

Estimate 0.4055 -0.4055 19.1606

Std Error 0.9129 1.6833 4809.3409

z value 0.444 -0.241 0.004

Tabel 25: Resulaten test 3 van model 4 De resultaten voor model 5: Waarden: zie tabel 26 AIC: 17.045

51

Pr(> |z|) 0.657 0.810 0.997

Model 5 Intercept f2Laag

Estimate 0.6931 0.9163

Std Error 0.8660 1.3964

z value 0.800 0.656

Pr(> |z|) 0.423 0.512

Tabel 26: Resulaten test 3 van model 5 Ook nu zijn er weer verschillen te zien tussen de waarden in de tabellen van de drie modellen. Duidelijk is te zien dat de waarde f1G:f2Laag niet voorkomt in de dataset door de NA-waarde in de tabellen. In de laatste kolom van de tabel, de kolom met de p-waarden, zijn nu ook verschillen te zien. Dit komt doordat in de dataset de verschillende combinaties niet meer even vaak voorkomen. Volgens de AIC-waarden blijkt in eerste instantie model 1 weer het beste model. Na toevoeging van model 4 en 5, blijkt model 4 het beste te zijn voor deze test.

A.2

Freedman Experiment

set.seed(12) #Verander het getal om verschillende datasets te krijgen #Waarden gebruikt: 12, 46, 65, 32, 87, 29, 20, 21, 49, 13 #dataset genereren n<-100 X <- matrix(rnorm(n*50),n,50) y <- matrix(rnorm(n),n,1) freed <- data.frame(X,y) freed.lm <-lm(y~X, data=freed) summary(freed.lm) #9 van de 50 variabelen zijn significant op 0.25 #Nogmaals het experiment, maar nu met deze 9 variabelen. XX <- X[,c(17,20,28,33,34,35,38,44,50)] freed2 <- data.frame(XX,y) freed2.lm <-lm(y~XX, data=freed2) summary(freed2.lm) # 8 van de 9 variabelen significant op 0.25

A.3

Nabootsing TNT dataset

set.seed(140) #Verander het getal om verschillende datasets te krijgen #Waarden gebruikt: 14, 46, 154, 12, 86, 99, 23, 209, 6, 140 n<-1000 y<-sample(c(0,1),n,replace=TRUE) x1<-sample(c("A", "B", "C", "D", "E", "F", "G", "H", "I",

52

"J"),n,replace=TRUE) x2<-sample(c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"),n,replace=TRUE) x3<-sample(c("A", "B", "C"),n,replace=TRUE) x4<-sample(c("A", "B", "C", "D", "E", "F"),n,replace=TRUE) x5<-sample(c("A", "B", "C", "D", "E", "F"),n,replace=TRUE) x6<-sample(c("A", "B", "C", "D"),n,replace=TRUE) x7<-sample(c("A", "B", "C", "D"),n,replace=TRUE) x8<-sample(c("A", "B"),n,replace=TRUE) x9<-sample(c("A", "B"),n,replace=TRUE) x10<-sample(c("A", "B", "C"),n,replace=TRUE) x11<-sample(c("A", "B", "C"),n,replace=TRUE) x12<-sample(c("A", "B", "C"),n,replace=TRUE) x13<-sample(c("A", "B", "C"),n,replace=TRUE) d <- data.frame(y,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11, x12,x13) #Volledige model mod <-glm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13, family=binomial, data=d) anova(mod, test="Chisq") #Variabelen met p-waarden groter dan 0.25 eruit halen mod1<-glm(y~x1+x3+x9+x11, family=binomial, data=d) anova(mod1, test="Chisq")

A.4

Voorbeeld Bayes

In deze appendix is de berekening van de maximum likelihoods en het script van het voorbeeld in paragraaf 5.4 te vinden. A.4.1

Berekening MLE

Berekening van de maximum likelihood schatter (MLE) van model 1 en model 2. Deze wordt voor model 1 verkregen door de afgeleide naar λ te berekenen van de kansdichtheidsfunctie van model 1 en deze gelijk te stellen aan nul. Dus ! n d Y −λ λxi 0 = e dλ xi ! i=1   1 d = Qn e−nλ λΣxi i=1 xi ! dλ

53

=

1

h

Qn

i=1 xi !

−ne−nλ λΣxi + e−nλ Σxi λΣxi −1

i

Hieruit blijkt dat e−nλ Σxi λΣxi −1 Qn i=1 xi !

ne−nλ λΣxi Qn i=1 xi !

=

ne−nλ λΣxi

= e−nλ Σxi λΣxi λ−1

n = Σxi λ−1 Σxi λ = n ˆ = Σxi = X ¯ n . Voor model 2 geldt eenzelfde Dus de MLE voor model 1 is λ n berekening. De MLE wordt verkregen door de afgeleide naar θ te berekenen van de kansdichtheidsfunctie van model 2 en deze gelijk te stellen aan nul. Dus ! n   d Y 10 xi 0 = θ (1 − θ)10−xi dθ xi i=1 n   Y
10 d Σxi = θ (1 − θ)10n−Σxi xi dθ i=1   n Y  10  = Σxi θΣxi −1 (1 − θ)10n−Σxi + θΣxi − 1(10n − Σxi )(1 − θ)10n−Σxi −1 xi i=1

Hieruit blijkt dat n   Y 10 i=1

xi

Σxi θ

Σxi θ

Σxi −1

Σxi −1

θ

n   Y 10 Σxi = θ (10n − Σxi )(1 − θ)10n−Σxi −1 xi

10n−Σxi

(1 − θ)

i=1 Σxi

10n−Σxi

(1 − θ)

= θ

(10n − Σxi )(1 − θ)10n−Σxi (1 − θ)−1

Σxi θ−1 = (10n − Σxi )(1 − θ)−1 1−θ 10n − Σxi = θ Σxi Σxi θ = 10n Dus de MLE voor model 2 is θˆ = A.4.2

Σxi 10n

¯ n /10. =X

Script in R alleen B21

rm(list = ls()); graphics.off(); library(BMA) set.seed(12) n10<-10; n100<-100; n1000<-1000 lambda05<-0.5; lambda1<-1; lambda2<-2 54

B21.10.05 <- numeric(1000);B21.10.1 <- numeric(1000) B21.10.2 <- numeric(1000);B21.100.05 <- numeric(1000) B21.100.1 <- numeric(1000);B21.100.2 <- numeric(1000) B21.1000.05 <- numeric(1000);B21.1000.1 <- numeric(1000) B21.1000.2 <- numeric(1000) count10.05<-0;count10.1<-0;count10.2<-0 count100.05<-0;count100.1<-0;count100.2<-0 count1000.05<-0;count1000.1<-0;count1000.2<-0 for(i in 1:1000){ data.pois10.05<-rpois(n10,lambda05) data.pois10.1<-rpois(n10,lambda1) data.pois10.2<-rpois(n10,lambda2) data.pois100.05<-rpois(n100,lambda05) data.pois100.1<-rpois(n100,lambda1) data.pois100.2<-rpois(n100,lambda2) data.pois1000.05<-rpois(n1000,lambda05) data.pois1000.1<-rpois(n1000,lambda1) data.pois1000.2<-rpois(n1000,lambda2) #Model1: Pois(lambda) verdeling, met lambda = Exp(1) #Model2: Bin(10,theta) verdeling, met theta=Beta(alpha, beta) a <- 2; b <- 2 somX10.05 <- sum(data.pois10.05) logB21.10.05 <- sum(log(factorial(10)/factorial(10data.pois10.05)))+(somX10.05+1)*(log(n10+1))+lbeta(somX10.05 +a,10*n10-somX10.05+b)-lbeta(a,b)-lgamma(somX10.05+1) B21.10.05a <- exp(logB21.10.05) B21.10.05b <- ifelse(is.nan(B21.10.05a), 0, B21.10.05a) count10.05<- ifelse(B21.10.05b<1, count10.05+1,count10.05) B21.10.05[i]<- B21.10.05b somX10.1 <- sum(data.pois10.1) logB21.10.1 <- sum(log(factorial(10)/factorial(10data.pois10.1)))+(somX10.1+1)*(log(n10+1))+lbeta(somX10.1 +a, 10*n10-somX10.1+b)-lbeta(a,b)-lgamma(somX10.1+1) B21.10.1a <- exp(logB21.10.1) B21.10.1b <- ifelse(is.nan(B21.10.1a), 0, B21.10.1a) count10.1<- ifelse(B21.10.1b<1, count10.1+1,count10.1) B21.10.1[i] <- B21.10.1b somX10.2 <- sum(data.pois10.2) logB21.10.2 <- sum(log(factorial(10)/factorial(10data.pois10.2)))+(somX10.2+1)*(log(n10+1))+lbeta(somX10.2 +a, 10*n10-somX10.2+b)-lbeta(a,b)-lgamma(somX10.2+1) 55

B21.10.2a <- exp(logB21.10.2) B21.10.2b <- ifelse(is.nan(B21.10.2a), 0, B21.10.2a) count10.2 <- ifelse(B21.10.2b<1, count10.2+1,count10.2) B21.10.2[i] <- B21.10.2b somX100.05 <- sum(data.pois100.05) logB21.100.05 <- sum(log(factorial(10)/factorial(10data.pois100.05)))+(somX100.05+1)*(log(n100+1))+lbeta( somX100.05+a, 10*n100-somX100.05+b)-lbeta(a,b)-lgamma( somX100.05+1) B21.100.05a <- exp(logB21.100.05) B21.100.05b <- ifelse(is.nan(B21.100.05a), 0, B21.100.05a) count100.05<- ifelse(B21.100.05b<1, count100.05+1, count100.05) B21.100.05[i] <- B21.100.05b somX100.1 <- sum(data.pois100.1) logB21.100.1 <- sum(log(factorial(10)/factorial(10data.pois100.1)))+(somX100.1+1)*(log(n100+1))+lbeta( somX100.1+a, 10*n100-somX100.1+b)-lbeta(a,b)-lgamma( somX100.1+1) B21.100.1a <- exp(logB21.100.1) B21.100.1b <- ifelse(is.nan(B21.100.1a), 0, B21.100.1a) count100.1<- ifelse(B21.100.1b<1, count100.1+1,count100.1) B21.100.1[i] <- B21.100.1b somX100.2 <- sum(data.pois100.2) logB21.100.2 <- sum(log(factorial(10)/factorial(10data.pois100.2)))+(somX100.2+1)*(log(n100+1))+lbeta( somX100.2+a, 10*n100-somX100.2+b)-lbeta(a,b)-lgamma( somX100.2+1) B21.100.2a <- exp(logB21.100.2) B21.100.2b <- ifelse(is.nan(B21.100.2a), 0, B21.100.2a) count100.2<- ifelse(B21.100.2b<1, count100.2+1,count100.2) B21.100.2[i] <- B21.100.2b somX1000.05 <- sum(data.pois1000.05) logB21.1000.05<-sum(log(factorial(10)/factorial(10data.pois1000.05)))+(somX1000.05+1)*(log(n1000+1))+lbeta( somX1000.05+a, 10*n1000-somX1000.05+b)-lbeta(a,b)-lgamma( somX1000.05+1) B21.1000.05a <- exp(logB21.1000.05) B21.1000.05b <- ifelse(is.nan(B21.1000.05a), 0, B21.1000.05a) count1000.05<- ifelse(B21.1000.05b<1, count1000.05+1, count1000.05) B21.1000.05[i] <- B21.1000.05b 56

somX1000.1 <- sum(data.pois1000.1) logB21.1000.1 <- sum(log(factorial(10)/factorial(10data.pois1000.1)))+(somX1000.1+1)*(log(n1000+1))+lbeta( somX1000.1+a, 10*n1000-somX1000.1+b)-lbeta(a,b)-lgamma( somX1000.1+1) B21.1000.1a <- exp(logB21.1000.1) B21.1000.1b <- ifelse(is.nan(B21.1000.1a), 0, B21.1000.1a) count1000.1<- ifelse(B21.1000.1b<1, count1000.1+1,count1000.1) B21.1000.1[i] <- B21.1000.1b somX1000.2 <- sum(data.pois1000.2) logB21.1000.2 <- sum(log(factorial(10)/factorial(10data.pois1000.2)))+(somX1000.2+1)*(log(n1000+1))+lbeta( somX1000.2+a, 10*n1000-somX1000.2+b)-lbeta(a,b)-lgamma( somX1000.2+1) B21.1000.2a <- exp(logB21.1000.2) B21.1000.2b <- ifelse(is.nan(B21.1000.2a), 0, B21.1000.2a) count1000.2<- ifelse(B21.1000.2b<1, count1000.2+1,count1000.2) B21.1000.2[i] <- B21.1000.2b} par(mfrow=c(3,3)) #Histogram van B21 h10.05<-hist(B21.10.05, plot=F, breaks=50) plot(h10.05$mids, h10.05$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=10, Pois(0.5)", type="h") h10.1<-hist(B21.10.1, plot=F, breaks=50) plot(h10.1$mids, h10.1$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=10, Pois(1)", type="h") h10.2<-hist(B21.10.2, plot=F, breaks=50) plot(h10.2$mids, h10.2$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=10, Pois(2)", type="h") h100.05<-hist(B21.100.05, plot=F, breaks=50) plot(h100.05$mids, h100.05$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=100, Pois(0.5)", type="h") h100.1<-hist(B21.100.1, plot=F, breaks=50) plot(h100.1$mids, h100.1$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=100, Pois(1)", type="h") h100.2<-hist(B21.100.2, plot=F, breaks=50) plot(h100.2$mids, h100.2$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=100, Pois(2)", type="h") h1000.05<-hist(B21.1000.05, plot=F, breaks=50) plot(h1000.05$mids, h1000.05$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=1000, Pois(0.5)", type="h") 57

h1000.1<-hist(B21.1000.1, plot=F, breaks=50) plot(h1000.1$mids, h1000.1$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=1000, Pois(1)", type="h") h1000.2<-hist(B21.1000.2, plot=F, breaks=50) plot(h1000.2$mids, h1000.2$counts, log="x", pch=20, xlab="B21", ylab="Aantal", main="n=1000, Pois(2)", type="h") ratio10.05<-count10.05/1000; ratio10.05 ratio10.1<-count10.1/1000; ratio10.1 ratio10.2<-count10.2/1000; ratio10.2 ratio100.05<-count100.05/1000; ratio100.05 ratio100.1<-count100.1/1000; ratio100.1 ratio100.2<-count100.2/1000; ratio100.2 ratio1000.05<-count1000.05/1000; ratio1000.05 ratio1000.1<-count1000.1/1000; ratio1000.1 ratio1000.2<-count1000.2/1000; ratio1000.2 A.4.3

Script in R, B21 met BIC

rm(list = ls()); graphics.off();library(BMA) set.seed(12) n10<-10; n100<-100; n1000<-1000 lambda05<-0.5; lambda1<-1; lambda2<-2 appB21.10.05 <- numeric(1000);appB21.10.1 <- numeric(1000) appB21.10.2 <- numeric(1000);appB21.100.05 <- numeric(1000) appB21.100.1 <- numeric(1000);appB21.100.2 <- numeric(1000) appB21.1000.05 <- numeric(1000);appB21.1000.1 <- numeric(1000) appB21.1000.2 <- numeric(1000) count10.05<-0;count10.1<-0;count10.2<-0 count100.05<-0;count100.1<-0;count100.2<-0 count1000.05<-0;count1000.1<-0;count1000.2<-0 for(i in 1:1000){ data.pois10.05<-rpois(n10,lambda05) data.pois10.1<-rpois(n10,lambda1) data.pois10.2<-rpois(n10,lambda2) data.pois100.05<-rpois(n100,lambda05) data.pois100.1<-rpois(n100,lambda1) data.pois100.2<-rpois(n100,lambda2) data.pois1000.05<-rpois(n1000,lambda05) data.pois1000.1<-rpois(n1000,lambda1) data.pois1000.2<-rpois(n1000,lambda2) somX10.05 <- sum(data.pois10.05) 58

A10.5 <- factorial(10)/factorial(10-data.pois10.05) B10.5 <- (10-somX10.05/n10)^(10-data.pois10.05) / (10^10) logB21.10.05 <- prod(A10.5*B10.5*exp(somX10.05/n10)) appB21.10.05a <- ifelse(is.nan(logB21.10.05), 0, logB21.10.05) count10.05<- ifelse(appB21.10.05a<1, count10.05+1, count10.05) appB21.10.05[i]<- appB21.10.05a somX10.1 <- sum(data.pois10.1) A10.1 <- factorial(10)/factorial(10-data.pois10.1) B10.1 <- (10-somX10.1/n10)^(10-data.pois10.1) / (10^10) logB21.10.1 <- prod(A10.1*B10.1*exp(somX10.1/n10)) appB21.10.1a <- ifelse(is.nan(logB21.10.1), 0, logB21.10.1) count10.1<- ifelse(appB21.10.1a<1, count10.1+1, count10.1) appB21.10.1[i]<- appB21.10.1a somX10.2 <- sum(data.pois10.2) A10.2 <- factorial(10)/factorial(10-data.pois10.2) B10.2 <- (10-somX10.2/n10)^(10-data.pois10.2) / (10^10) logB21.10.2 <- prod(A10.2*B10.2*exp(somX10.2/n10)) appB21.10.2a <- ifelse(is.nan(logB21.10.2), 0, logB21.10.2) count10.2<- ifelse(appB21.10.2a<1, count10.2+1, count10.2) appB21.10.2[i]<- appB21.10.2a somX100.05 <- sum(data.pois100.05) A100.05 <- factorial(10)/factorial(10-data.pois100.05) B100.05 <- (10-somX100.05/n100)^(10-data.pois100.05)/(10^10) logB21.100.05 <- prod(A100.05*B100.05*exp(somX100.05/n100)) appB21.100.05a <- ifelse(is.nan(logB21.100.05), 0, logB21.100.05) count100.05<- ifelse(appB21.100.05a<1, count100.05+1, count100.05) appB21.100.05[i]<- appB21.100.05a somX100.1 <- sum(data.pois100.1) A100.1 <- factorial(10)/factorial(10-data.pois100.1) B100.1 <- (10-somX100.1/n100)^(10-data.pois100.1)/(10^10) logB21.100.1 <- prod(A100.1*B100.1*exp(somX100.1/n100)) appB21.100.1a <- ifelse(is.nan(logB21.100.1), 0, logB21.100.1) count100.1<- ifelse(appB21.100.1a<1, count100.1+1, count100.1) 59

appB21.100.1[i]<- appB21.100.1a somX100.2 <- sum(data.pois100.2) A100.2 <- factorial(10)/factorial(10-data.pois100.2) B100.2 <- (10-somX100.2/n100)^(10-data.pois100.2)/(10^10) logB21.100.2 <- prod(A100.2*B100.2*exp(somX100.2/n100)) appB21.100.2a <- ifelse(is.nan(logB21.100.2), 0, logB21.100.2) count100.2<- ifelse(appB21.100.2a<1, count100.2+1, count100.2) appB21.100.2[i]<- appB21.100.2a somX1000.05 <- sum(data.pois1000.05) A1000.05 <- factorial(10)/factorial(10-data.pois1000.05) B1000.05 <- (10-somX1000.05/n1000)^(10-data.pois1000.05) /(10^10) logB21.1000.05 <- prod(A1000.05*B1000.05*exp( somX1000.05/n1000)) appB21.1000.05a <- ifelse(is.nan(logB21.1000.05), 0, logB21.1000.05) count1000.05<- ifelse(appB21.1000.05a<1, count1000.05+1, count1000.05) appB21.1000.05[i]<- appB21.1000.05a somX1000.1 <- sum(data.pois1000.1) A1000.1 <- factorial(10)/factorial(10-data.pois1000.1) B1000.1 <- (10-somX1000.1/n1000)^(10-data.pois1000.1) /(10^10) logB21.1000.1 <- prod(A1000.1*B1000.1*exp( somX1000.1/n1000)) appB21.1000.1a <- ifelse(is.nan(logB21.1000.1), 0, logB21.1000.1) count1000.1<- ifelse(appB21.1000.1a<1, count1000.1+1, count1000.1) appB21.1000.1[i]<- appB21.1000.1a somX1000.2 <- sum(data.pois1000.2) A1000.2 <- factorial(10)/factorial(10-data.pois1000.2) B1000.2 <- (10-somX1000.2/n1000)^(10-data.pois1000.2) /(10^10) logB21.1000.2 <- prod(A1000.2*B1000.2*exp( somX1000.2/n1000)) appB21.1000.2a <- ifelse(is.nan(logB21.1000.2), 0, logB21.1000.2) count1000.2<- ifelse(appB21.1000.2a<1, count1000.2+1, count1000.2) appB21.1000.2[i]<- appB21.1000.2a 60

} par(mfrow=c(3,3)) #Histogram van B21 met BIC h10.05<-hist(appB21.10.05, plot=F, breaks=50) plot(h10.05$mids, h10.05$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=10, Pois(0.5)", type="h") h10.1<-hist(appB21.10.1, plot=F, breaks=50) plot(h10.1$mids, h10.1$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=10, Pois(1)", type="h") h10.2<-hist(appB21.10.2, plot=F, breaks=50) plot(h10.2$mids, h10.2$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=10, Pois(2)", type="h") h100.05<-hist(appB21.100.05, plot=F, breaks=50) plot(h100.05$mids, h100.05$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=100, Pois(0.5)", type="h") h100.1<-hist(appB21.100.1, plot=F, breaks=50) plot(h100.1$mids, h100.1$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=100, Pois(1)", type="h") h100.2<-hist(appB21.100.2, plot=F, breaks=50) plot(h100.2$mids, h100.2$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=100, Pois(2)", type="h") h1000.05<-hist(appB21.1000.05, plot=F, breaks=50) plot(h1000.05$mids, h1000.05$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=1000, Pois(0.5)", type="h") h1000.1<-hist(appB21.1000.1, plot=F, breaks=50) plot(h1000.1$mids, h1000.1$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=1000, Pois(1)", type="h") h1000.2<-hist(appB21.1000.2, plot=F, breaks=50) plot(h1000.2$mids, h1000.2$counts, log="x", pch=20, xlab="B21^BIC", ylab="Aantal", main="n=1000, Pois(2)", type="h") ratio10.05<-count10.05/1000; ratio10.05 61

ratio10.1<-count10.1/1000; ratio10.1 ratio10.2<-count10.2/1000; ratio10.2 ratio100.05<-count100.05/1000; ratio100.05 ratio100.1<-count100.1/1000; ratio100.1 ratio100.2<-count100.2/1000; ratio100.2 ratio1000.05<-count1000.05/1000; ratio1000.05 ratio1000.1<-count1000.1/1000; ratio1000.1 ratio1000.2<-count1000.2/1000; ratio1000.2

62

B

Script onderzoek

rm(list = ls()); graphics.off(); library(BMA) setwd("G:/BachelorProject/") d2007<-read.csv2("data2007bep.csv", header=TRUE) d2007<-as.data.frame(d2007,stringsAsFactors=TRUE) d2008<-read.csv2("data2008bep.csv", header=TRUE) d2008<-as.data.frame(d2008,stringsAsFactors=TRUE) #alles met NA eruit halen dd2007<-d2007[,-c(15,16,17)] dd2008<-d2008[,-c(15,16,17)] #formule met InfoPRIC f<-formula(Optijd ~ VerzendenDvdW + Gewicht + Formaat + Schrift + InterofIntra + Streepjescode + Poststroom + InfoPRIC) glm.out.T<-bic.glm(f,data=dd2007,glm.family="binomial",factor.type=TRUE) summary(glm.out.T) plot(glm.out.T, mfrow=c(4,4)) layout(matrix(c(0,1), 1, 2, byrow = TRUE), widths=c(1,7)) imageplot.bma(glm.out.T, color = c("red", "blue", "#FFFFD5")) #formule zonder InfoPRIC f2<-formula(Optijd ~ VerzendenDvdW + Gewicht + Formaat + Schrift + InterofIntra + Streepjescode + Poststroom) glm.out.T2<-bic.glm(f2,data=dd2007,glm.family="binomial",factor.type=TRUE) summary(glm.out.T2) plot(glm.out.T2, mfrow=c(3,4)) layout(matrix(c(0,1), 1, 2, byrow = TRUE), widths=c(1,7)) imageplot.bma(glm.out.T2) #nu kunnen ook aparte levels van een factor uit of in het model zitten glm.out.F2<-bic.glm(f2,data=dd2007,glm.family="binomial",factor.type=FALSE) summary(glm.out.F2) plot(glm.out.F2, mfrow=c(3,4)) layout(matrix(c(0,1), 1, 2, byrow = TRUE), widths=c(1,7)) imageplot.bma(glm.out.F2) #Testen data model.2007<-glm(f2,family="binomial", data=dd2007) summary(model.2007) model.2008<-glm(f2,family="binomial", data=dd2008) summary(model.2008)

63

Referenties [1] Adrian E. Raftery, Ian S. Painter and Christopher T. Volinsky, BMA: AN R package for Baysian Model Averaging, R News Vol 5/2, november 2005, 2-8 [2] Adrian E. Raftery, Baysian Model Selection in Social Research, Sociological Methodology Vol. 25, 1995, 11-163 [3] S.N. Wood, (2006) Generalized Additive Models: An Introduction with R, Chapman & Hall/CRC, 2006, 1-119 [4] Fred. C. Pampel, Logistic Regression, Sage Publications, 2000 [5] John Fox, Applied Regression Analysis and Generalized Linear Models, Sage Publications, 2008, 124-128 & 611-626 [6] TNT, Codeerinformatie op de brief

64