Technische Universität M ünchen - mediaTUM · Technische Universität M ünchen Zentrum...

Technische Universitat Munchen

Zentrum Mathematik

Vergleich nicht-genesteter

Regressionsmodelle

fur Zahldaten

mit Hilfe von Bayes-Faktoren

Diplomarbeit

von

Anna Katharina Fendt

Themensteller: Prof. Dr. C. Czado, Dr. G. Sußmann (VKB)

Betreuer: Prof. Dr. C. Czado, Dr. G. Sußmann (VKB)

Abgabetermin: 15. April 2004

Hiermit erklare ich, dass ich die Diplomarbeit selbststandig angefertigt und nur die angegebenen

Quellen verwendet habe.

Munchen, 15. April 2004

Inhaltsverzeichnis

1 Einleitung 1

2 Theoretische Grundlagen 3

2.1 Die Poissonverteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 Generalisierte Lineare Modelle (GLMs) . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.1 Die Exponentielle Familie . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2.2 Die Komponenten eines GLM . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2.3 Offset in der Poissonregression . . . . . . . . . . . . . . . . . . . . . . . . 8

2.3 Maximum-Likelihood (ML)-Schatzung . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3.1 Beschreibung der Maximum-Likelihood-Schatzung fur den Fall

unabhangiger Beobachtungen . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3.2 Iterative Algorithmen zur numerischen Bestimmung des ML-Schatzers . . 12

2.4 Asymptotische Hypothesentests in GLMs . . . . . . . . . . . . . . . . . . . . . . 21

2.5 Goodness-of-Fit-Maße . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.5.1 Devianzanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.5.2 Residualanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.6 Die Negativ-Binomial-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.6.1 Definition und Eigenschaften des NB-Verteilung . . . . . . . . . . . . . . 28

2.6.2 Die NB-Verteilung mit bekanntem Parameter a als GLM . . . . . . . . . 33

2.6.3 Offset in der NB-Regression . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.6.4 Unbeobachtete Heterogenitat und gemischte Modelle . . . . . . . . . . . . 35

2.6.5 Fisher-Informationsmatrix und Asymptotik der Schatzer . . . . . . . . . . 36

3 Bayesianische Modelle und Bayes-Faktoren zur Modellwahl 40

3.1 Bayesianischer Ansatz und Definition der Bayes-Faktoren . . . . . . . . . . . . . 42

3.2 Approximation der Bayes-Faktoren . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.2.1 Die Laplace-Methode fur Integrale . . . . . . . . . . . . . . . . . . . . . . 47

3.2.2 Drei weitere Approximationen . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.3 Anwendung auf Generalisierte Lineare Modelle . . . . . . . . . . . . . . . . . . . 56

i

Inhaltsverzeichnis 1

3.4 Wahl der Gestalt der a priori Verteilung fur die Parameter in einem GLM . . . . 58

3.4.1 A Priori Verteilungen in Linearen Modellen . . . . . . . . . . . . . . . . . 58

3.4.2 A Priori Verteilungen in gewichteten Linearen Modellen . . . . . . . . . . 61

3.4.3 A Priori Verteilungen in GLMs . . . . . . . . . . . . . . . . . . . . . . . . 62

3.4.4 Wahl der Hyper-Parameter in der a priori Verteilung der Regressionspa-

rameter in GLMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.5 Unbekannter Dispersionsparameter und Uberdispersion . . . . . . . . . . . . . . 71

4 Anwendungsbeispiele 73

4.1 Patent-Daten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4.1.1 Explorative Datenanalyse der Patent-Daten . . . . . . . . . . . . . . . . . 73

4.1.2 Modellvergleiche mit Hilfe der Bayes-Faktoren . . . . . . . . . . . . . . . 74

5 Zusammenfassung und Ausblick 88

Anhang 90

Literaturverzeichnis 100

Abbildungsverzeichnis 102

Tabellenverzeichnis 103

Kapitel 1

Einleitung

Versicherungsunternehmen versuchen mit ihren Tarifen die Realitat moglichst genau abzubil-

den. Bei der Entwicklung eines neuen Tarifes ist man speziell in der Kfz-Versicherung daran

interessiert, wieviele Schaden im Durchschnitt fur das nachste Jahr zu erwarten sind und wie

sich Gruppen von Fahrzeughaltern, die bestimmte Merkmale aufweisen, hinsichtlich ihrer Scha-

denhaufigkeit unterscheiden. Grundsatzlich stellt sich immer die Frage, woran man einen Versi-

cherungsnehmer mit hohem oder niedrigem Unfallrisiko erkennt. Deshalb werden verschiedene

Großen wie das Alter des Fahrers, das Geschlecht des Fahrers, die Haltedauer, usw. erhoben. In

Abhangigkeit von diesen Merkmalen wird dann versucht, die interessierende Große Schadenan-

zahl zu erklaren. Dazu muss zunachst ein statistisches Modell erstellt werden. Dieses ist so zu

wahlen, dass es die beobachteten Daten in ihrer wesentlichen Struktur wiedergibt. Da es sich bei

der zu analysierenden Große Schadenanzahl um eine Zahlgroße handelt, dient ublicherweise das

Poissonmodell als Grundlage. Wegen seiner einem Zahldatenmodell entsprechenden Eigenschaf-

ten, eignet es sich sehr gut zur Modellierung der Schadenanzahl. Zu Beginn des zweiten Kapitels

werden deshalb kurz die charakteristischen Merkmale der Poissonverteilung zusammengefasst.

Mit die wichtigste Eigenschaft der Poissonverteilung ist die Gleichheit von Erwartungswert

und Varianz. Nun kann es aber vorkommen, dass genau diese Modellannahme eine gute An-

passung an die Daten verhindert. In der Praxis streuen die Zahldaten, im Beispiel der Kfz-

Versicherung die Schadenanzahl, meist in großerem Maße um ihren Erwartungswert, als es das

Poissonmodell erlaubt. In diesem Fall liegt eine Modellverletzung durch Uberdispersion in den

Daten vor. Weisen die betrachteten Daten Uberdispersion auf, so besteht die Moglichkeit, von

der Poissonverteilung auf die Negativ-Binomial-Verteilung uberzugehen. Diese Verteilung hat

im Vergleich zur Poissonverteilung einen Parameter mehr und enthalt die Poissonverteilung als

Spezialfall. Die Negativ-Binomial-Verteilung stellt keine so starke Forderung bzgl. der Beziehung

zwischen Erwartungswert und Varianz an die Daten. Nach der Einfuhrung der Generalisierten

Linearen Modelle (GLM), auf deren Basis die Anpassung der Daten geschehen soll, wird deshalb

auf die charakteristischen Eigenschaften der NB-Verteilung eingegangen.

1

2 Kapitel 1. Einleitung

In vielen Fallen und speziell im Rahmen der Tarifierung genugt es aber nicht, ein einziges

Poissonmodell an die Zahldaten anzupassen und bei eventuell auftretender Uberdispersion auf

ein NB-Modell uberzugehen. Beispielsweise fließen in die Berechnung eines neuen Tarifs oft auch

neu zu testende Merkmale ein. Dann ist es notig, Modelle mit herkommlichen Merkmalen mit

Modellen, in denen ein herkommliches Merkmal durch ein neu erhobenes ersetzt wurde, hinsicht-

lich ihrer Anpassungsgute zu vergleichen. Dadurch soll festgestellt werden, ob das neue Merkmal

mehr zu einer guten Modellanpassung beitragt und deshalb das herkommliche Merkmal ersetzen

sollte. Dies fuhrt dazu, dass die betrachteten Modelle nicht genestet sind. Zwei Modelle bezeich-

net man als genestet, wenn ein Modell als Spezialfall im anderen enthalten ist. Nicht-genestete

Modelle lassen sich mit herkommlichen auf p-Werten basierenden Signifikanztests nicht hinsicht-

lich ihrer Anpassungsgute vergleichen. Diese Tatsache motiviert die vorliegende Diplomarbeit.

Die Berechnung sogenannter Bayes-Faktoren zum Vergleich verschiedener Regressionsmo-

delle ermoglicht es, Aussagen daruber zu treffen, welches Modell die Daten besser anpasst,

unabhangig davon, ob die betrachteten Modelle genestet sind oder nicht. Basierend auf dem

Artikel Approximate Bayes factors and accounting for model uncertainty in generalised line-

ar models von A.E. Raftery [Raft 96] werden im dritten Kapitel die Bayes-Faktoren fur GLMs

definiert und verschiedene Approximationen hergeleitet. Eine analytische Berechnung der Bayes-

Faktoren ist nur in einigen elementaren Fallen moglich. Die Berechnung der Approximationen

vereinfacht sich, wenn man a priori, also vor der Datenerhebung, davon ausgeht, dass die Regres-

sionsparameter Normal-verteilt sind. Auf die genaue Wahl der Gestalt der a priori Verteilung

und deren Parameter wird im Anschluss an die Herleitung der Approximationen eingegangen.

Im Anwendungskapitel werden dann fur konkrete Datenbeispiele approximierte Bayes-Faktoren

zum Vergleich genesteter und vor allem auch nicht-genesteter Regressionsmodelle berechnet und

dazu benutzt, aus einer Gruppe von Regressionsmodellen dasjenige zu bestimmen, das die beste

Anpassung an die Daten liefert.

Kapitel 2

Theoretische Grundlagen

2.1 Die Poissonverteilung

Zur Analyse von Zahldaten wie z.B. der Anzahl von Schadensmeldungen bei einer Sachversi-

cherung bietet sich die Poissonverteilung an. In diesem Abschnitt werden deshalb die Wahr-

scheinlichkeitsfunktion, die momenterzeugende Funktion und einige wichtige Eigenschaften der

Poissonverteilung vorgestellt. Die Poissonverteilung gehort zur Klasse der diskreten Verteilun-

gen.

Definition 2.1 (Poissonverteilung) Sei Y eine Zufallsvariable mit einer diskreten Vertei-

lung, die definiert ist auf N ∪ 0 = 0, 1, 2, ....Y ist dann Poisson-verteilt mit Parameter (oder Rate) λ, kurz Y ∼ Poi(λ), falls die Wahrschein-

lichkeitsfunktion folgende Form hat:

P (Y = y) =e−λλy

y!fur λ ∈ R

+ und y = 0,1,2, . . .

Beispiel 2.1 (Schadensfalle bei einer Versicherung) (vgl. [Fahr 99],Seite 261)

Eine Versicherung will die Pramien fur Versicherungen gegen Großunfalle kalkulieren. Erfah-

rungswerte fuhren zu der Annahme, dass die Zufallsvariable

Y = ”Anzahl der Großunfalle im Winterhalbjahr (Oktober bis Marz)“

Poisson-verteilt ist mit Rate λ = 3.

Damit lasst sich nun die Wahrscheinlichkeit fur genau zwei Großunfalle in einem Winter be-

stimmen:

P (Y = 2) = e−3 32

2!≈ 0.22

3

4 Kapitel 2. Theoretische Grundlagen

Die Poissonverteilung zeichnet eine spezielle Eigenschaft aus:

E(Y ) = V ar(Y ) = λ

Die Gleichheit von Erwartungswert und Varianz bezeichnet man als Aquidispersion. Ist der Er-

wartungswert kleiner als die Varianz, spricht man von Uberdispersion. Und im umgekehrten Fall,

wenn also die Varianz kleiner als der Erwartungswert ist, liegt Unterdispersion vor.

Die Grafiken in Abbildung 2.1 sollen die Wahrscheinlichkeitsfunktion der Poissonverteilung fur

verschiedene λ-Werte veranschaulichen. Je kleiner λ ist, desto linkssteiler wird die Wahrschein-

lambda = 0.5

y

P(Y

=y)

0 2 4 6 8 10

0.0

0.1

0.2

0.3

0.4

0.5

0.6

lambda = 1

y

P(Y

=y)

0 2 4 6 8 10

0.0

0.1

0.2

0.3

lambda = 5

y

P(Y

=y)

0 5 10 15

0.0

0.05

0.10

0.15

lambda=10

y

P(Y

=y)

0 5 10 15 20

0.0

0.02

0.04

0.06

0.08

0.10

0.12

Abbildung 2.1: Darstellung der Poisson-Wahrscheinlichkeitsfunktion fur verschiedene λ-Werte

lichkeitsfunktion, und desto großer werden die Wahrscheinlichkeiten fur kleine Werte auf der

x-Achse. Fur großeres λ wird die Verteilung annahernd symmetrisch und lasst sich durch eine

Normalverteilungsdichte approximieren. (siehe [Fahr 99], Abschnitt 7.2)

2.1. Die Poissonverteilung 5

Satz 2.1 (Momenterzeugende Funktion der Poissonverteilung) Die momenterzeugende

Funktion der Poissonverteilung ist gegeben durch:

MY (t) = E[etY ] = eλ(et−1)

Beweis:

MY (t) = E[etY ] =∞∑y=0

etye−λλy

y!= e−λ

∞∑y=0

etyλy

y!= e−λ

∞∑y=0

(etλ)y

y!= e−λee

tλ = eλ(et−1).

Damit sind die Momente der Poissonverteilung definiert durch:

E[Y k] = M (k)(0) =∂kM(0)∂tk

, k = 0, 1, 2, . . . (2.1)

Der Erwartungswert lasst sich durch die Berechnung des ersten Moments bestatigen:

E[Y ] = M ′(0) = eλ(et−1)λet |t=0= λ

Fur die Varianz gilt nach [Rade 97] (Seite 420):

V ar(Y ) = M ′′(0) − [M ′(0)]2 (2.2)

Diese Formel ermoglicht uns nun auch eine Bestatigung der Varianz:

V ar(y) = M ′′(0) − [M ′(0)]2 = eλ(et−1)λetλet + eλ(et−1)λet|t=0 − [λ]2 = λ+ λ2 − [λ]2 = λ

Die Poissonverteilung besitzt ausserdem noch die Eigenschaft der Additivitat.

Satz 2.2 (Additivitat der Poissonverteilung) Es seien Yi ∼ Poi(λi), i = 1, 2, . . ., unabhangi-

ge Zufallsvariablen. Weiter gelte∑∞

i=1 λi <∞. Dann gilt:

Z =∞∑i=1

Yi ∼ Poi

( ∞∑i=1

λi

)

Beweis: Die Yi sind nach Voraussetzung unabhangig. Deshalb gilt fur die momenterzeugenden

Funktionen

MZ(t) = M∑∞i=1 Yi

(t) =∞∏i=1

MYi(t).

Ebenfalls nach Voraussetzung sind die Yi alle Poisson-verteilt. Demnach gilt:

MZ(t) =∞∏i=1

MYi(t) =∞∏i=1

eλi(et−1) = e∑∞

i=1 λi(et−1)

Dies stellt aber gerade die momenterzeugende Funktion einer Poisson-verteilten Zufallsvariable

mit Parameter∑∞

i=1 λi dar.


2.2 Generalisierte Lineare Modelle (GLMs)

Die Generalisierten Linearen Modelle bieten uns mehr Spielraum als die klassischen Linearen

Modelle. Wahrend die klassischen Linearen Modelle auf der Normalverteilung basieren, konnen

den GLMs viele verschiedene Verteilungen zugrundeliegen. Einzige Voraussetzung ist, dass die

gewahlte Verteilung zur linearen Exponentiellen Familie gehort. Zusatzlich ist in den GLMs

im Gegensatz zu den klassischen Linearen Modellen die Transformation des Erwartungswertes

erlaubt. Dies aussert sich in der Wahl der Linkfunktion. Eine ausfuhrliche Einfuhrung der GLMs

lasst sich z.B. in [McCu 89] finden.

2.2.1 Die Exponentielle Familie

Um zu entscheiden, ob eine Verteilung einem GLM zugrundegelegt werden darf, muss ihre Zu-

gehorigkeit zur linearen Exponentiellen Familie uberpruft werden.

Definition 2.2 (Lineare Exponentielle Familie) Eine Dichte f(y; θ, φ) gehort zur linearen

exponentiellen Familie mit naturlichem oder kanonischem Parameter θ und Stor- oder Skalen-

parameter φ, falls sie sich in der folgenden Form schreiben lasst:

f(y; θ, φ) = expθy − b(θ)a(φ)

+ c(y, φ).

Dabei beschreibt c(y, φ) eine normalisierende Konstante, und die Funktionen a(.) und b(.) be-

stimmen die jeweilige Verteilung.

Damit lautet die Log-Likelihood-Funktion

L(y; θ, φ) = ln f(y; θ, φ) =yθ − b(θ)a(φ)

+ c(y;φ).

Mit Hilfe der bekannten Gleichungen (siehe [Case 90], Seite 309ff)

E

(∂L∂θ

)= 0

und

E

(∂2L∂θ2

)+ E

(∂L∂θ

)2

= 0

erhalt man folgende Formeln fur den Erwartungswert und die Varianz:

E(Y ) = µ = b′(θ) (2.3)

V ar(Y ) = b′′(θ)a(φ) (2.4)

Dabei hangt die Funktion b′′(θ) nur vom kanonischen Parameter θ ab und wird auch als Vari-

anzfunktion v bezeichnet. Die Funktion a(φ) hingegen hangt nur vom Dispersionsparameter φ

ab. Eine genaue Herleitung dieser Formeln findet man in [Siko 02] (Seite 10f).

2.2. Generalisierte Lineare Modelle (GLMs) 7

Beispiel 2.2 (Die Poissonverteilung als Mitglied der Exponentiellen Familie) Da sich

die Dichte der Poissonverteilung in der Form

fpoi(y; θ, φ) = exp−λ+ y lnλ− ln(y!)

mit a(φ) = 1, θ = lnλ, b(θ) = eθ = elnλ = λ und c(y, φ) = − ln(y!) schreiben lasst, gehort

sie zur exponentiellen Familie. Aus den Formeln fur den Erwartungswert (2.3) und die Varianz

(2.4) ergibt sich

E(Y ) = µ = b′(θ) = eθ = λ

und

V ar(Y ) = b′′(θ)a(φ) = eθ = λ.

Damit bestatigt sich die Eigenschaft der Gleichheit von Erwartungswert und Varianz bei der

Poissonverteilung.

2.2.2 Die Komponenten eines GLM

Ein GLM besteht aus drei Komponenten, der Zufallskomponente, der systematischen Kompo-

nente und dem Link, der die Verbindung zwischen der Zufalls- und der systematischen Kompo-

nente darstellt. Um also ein GLM zu definieren, trifft man folgende Annahmen:

(i) Zufallskomponente bzw. Verteilungsannahme:

Die Zufallsvariablen Yi, i = 1, . . . , n, seien gegeben die Kovariablenvektoren xi, i = 1, . . . , p,

bedingt unabhangig verteilt nach einer Verteilung aus der exponentiellen Familie. Dies

stellt eine Erweiterung gegenuber den klassischen Linearen Modellen dar, in welchen aus-

schließlich die Normalverteilung angenommen wird.

(ii) Systematische Komponente bzw. Strukturelle Annahme: Im Folgenden sie i der

Index der n verschiedenen Beobachtungen und j der Index der p relevanten Kovariablen-

vektoren x1, . . . , xp. Der sogenannte Lineare Pradiktor η mit η = (η1, . . . , ηn)t ist gegeben

durch:

η = Xβ mit X = (x1, . . . ,xp) ∈ Rn×p und β = (β1, . . . , βp)t

In Komponentenschreibweise: ηi = xtiβ mit xi = (xi1, . . . , xip)t, i= 1,. . . ,n.

Der lineare Pradiktor fuhrt somit die Regressionsparameter β =(β1, . . . , βp)t in das Modell

ein.

(iii) Verbindung zwischen zufalliger und systematischer Komponente:

Nun mussen wir noch die zufallige mit der systematischen Komponente verbinden. Der

Erwartungswert µi wird uber eine Responsefunktion h mit dem linearen Pradiktor ηi ver-

knupft, d.h. µi = h(ηi) = h(xti β), i=1,. . . ,n.


Unter Verwendung der inversen Responsefunktion g := h−1(µi) erhalten wir:

g(µi) = ηi = xtiβ bzw. µi = g−1(ηi), i = 1, . . . , n

Die Funktion g: R → R wird im Allgemeinen als monoton und differenzierbar vorausgesetzt

und als Linkfunktion bezeichnet. Die Linkfunktion fungiert damit als Verbindung zwischen

zufalliger und systematischer Komponente.

Durch die Bedingung ηi = θi, i = 1, . . . , n, wird die kanonische Linkfunktion definiert. Fur diese

kanonische Linkfunktion existiert eine suffiziente Statistik.

Definition 2.3 (Suffiziente Statistiken) Seien X = (X1, . . . ,Xn)t Daten mit gemeinsamer

Dichte f(x, θ). Dann ist die Statistik T (X) suffizient fur den Parameter θ genau dann, wenn die

bedingte Verteilung von X gegeben T (X) = t unabhangig von θ ist, d.h. nachdem die Statistik T

bekannt ist, enthalten die Daten X keine weiteren Informationen uber den Parameter θ.

Proposition 2.1 (Suffiziente Statistik in einem GLM mit kanonischer Linkfunktion)

In einem GLM mit einer kanonischen Linkfunktion ist die Statistik (∑n

i=1 xi1yi, . . . ,∑n

i=1 xipyi)t

suffizient fur β = (β1, . . . , βp)t.

Beweis: mittels folgendem Faktorisierungstheorem

Satz 2.3 (Faktorisierungstheorem) Die Statistk T (X) mit Bild I ist suffizient fur den Pa-

rameter θ genau dann, wenn es Funktionen g(t, θ) definiert fur t ∈ I, θ ∈ Θ und h definiert auf

Rn gibt, so dass p(x, θ) = g(T (x), θ) · h(x) gilt.

Beweis: siehe [Bick 77], Seite 65.

Beispiel 2.3 (Die kanonische Linkfunktion der Poissonverteilung) Es gelte

(Yi|xi) ∼ Poi(µi), i=1,. . . ,n. Alle Beobachtungen seien unabhangig. Wie wir bereits aus Beispiel

2.2 wissen, gilt b(θ) = exp(θ). Daraus folgt µ = b′(θ) = exp(θ). D.h. die Bedingung η = θ ist

erfullt, wenn η = lnµ = ln(exp(θ)) = θ. Damit ist die Log-Linkfunktion g(µ) = ln(µ) die

kanonische Linkfunktion fur ein Regressionsmodell mit Poissonverteilung.

2.2.3 Offset in der Poissonregression

Bisher sind wir davon ausgegangen, dass die Zeitraume, in denen unsere Ereignisse auftreten,

einheitlich sind. Tatsachlich aber variieren in der Praxis die Beobachtungszeitraume sehr oft.

Deshalb fuhren wir eine neue Variable ti ein, die den jeweiligen Zeitraum der i-ten Beobachtung

wiedergibt, um die Zeiteinheit zu standardisieren. Somit ist unsere Responsevariable Yi Poisson-

verteilt mit Parameter (Rate) λ = tiµ∗i , d.h.:

P (Yi = yi) =e−tiµ∗i (tiµ∗i )

yi

yi!, i = 1, . . . , n, y = 0, 1, 2, . . .

2.3. Maximum-Likelihood (ML)-Schatzung 9

Hierbei gilt: µ∗i = exp(xtiβ).

Nun schreiben wir diese Dichte in der Form einer exponentiellen Familie (siehe Definition (2.2)):

fpoi(yi; θi, φ) = exp−tiµ∗i + yi ln(tiµ∗i ) − ln(yi!), i = 1, . . . , n.

Dabei gilt a(φ) = 1 = φ,

b(θi) = tiµ∗i ,

θi = ln(tiµ∗i ) und

c(y, φ)= − ln(yi!).

Daraus folgt: b(θi) = eθi mit θi = ln(tiµ∗i ), i = 1, . . . , n

Als Erwartungswert ergibt sich damit:

µi = E(Yi) = b′(θi) = eθi = tiµ∗i = ti exp(xt

iβ) = eln ti+xtiβ , i = 1, . . . , n

Dementsprechend gilt fur den linearen Pradiktor:

ηi = xtiβ = ln(µi) − ln(ti) = θi − ln(ti)︸︷︷︸

offset

, i = 1, . . . , n

Der letzte Term ln(ti) wird als offset bezeichnet und als bekannt vorausgesetzt. Der Offset ist

eine Regressionsvariable mit konstantem Koeffizienten 1 fur jede Beobachtung. Somit erhalt man

fur ti = 1 analog zu Beispiel 2.3 die kanonische Linkfunktion g(µi) = ln(µi), und im Falle ti = 1

die kanonische Linkfunktion

g(µi) = ln(µi) − ln(ti), i = 1, . . . , n.

2.3 Maximum-Likelihood (ML)-Schatzung

Im Anschluss an die Wahl eines geeigneten Modells mussen die Regressionsparameter geschatzt

werden. Dies erfolgt in den GLMs durch die Maximum-Likelihood-Methode. Um mit dieser Me-

thode die sogenannten ML-Schatzer fur die unbekannten Parameter zu bestimmen, ist es notig

vorauszusetzen, dass die zugrundeliegende Verteilung vollstandig bekannt ist und die Designma-

trix X der Kovariablen vollen Spaltenrang p besitzt.

Es sei allgemein f(y | θ) die Wahrscheinlichkeits- bzw. Dichtefunktion eines Zufallsvariablen-

vektors Y= (Y1, . . . , Yn)t, wobei der Parametervektor θ= (θ1, . . . , θm)t im zulassigen Parameter-

raum Θ liegt. Fur feste Realisierungen y nennt man f(y | θ) als Funktion von θ die ”Likelihood-

Funktion“ und schreibt L(θ).


2.3.1 Beschreibung der Maximum-Likelihood-Schatzung fur den Fall

unabhangiger Beobachtungen

Sind die Zufallsvariablen Y1, . . . , Yn gegeben die Kovariablenvektoren xi,i= 1,. . . ,n, unabhangig

diskret oder stetig verteilt mit Wahrscheinlichkeitsfunktion bzw. Dichte f(yi | xi,θ ), so lautet

wegen der Unabhangigkeit der Yi die Likelihood-Funktion

L(θ) =n∏i=1

f(yi | xi,θ ).

Bei der Maximum-Likelihood-Schatzung sucht man durch Maximieren der Likelihood-Funktion

den Wert θ des unbekannten Paramtervektors θ, fur den das Eintreten der beobachteten Stich-

probe maximale Wahrscheinlichkeit besitzt. Es muss also gelten:

L(θML) ≥ L(θ) ∀ θ ∈ Θ

Dabei bezeichnet Θ den Parameterraum.

Oft ist es allerdings leichter, die ”Log-Likelihood-Funktion“ L(θ)= lnL(θ) zu maximieren. Dies

ist erlaubt, da der Logarithmus eine monoton wachsende Funktion ist und sich damit durch das

Logarithmieren der gesuchte Maximalwert θML nicht verandert. Das Produkt in der Likelihood-

Funktion laßt sich in der Log-Likelihood-Funktion als Summe schreiben und erleichtert uns damit

das Differenzieren.

L(θ) = ln L(θ) = lnn∏i=1

f(yi | xi,θ) =n∑i=1

ln(f(yi|xi,θ))

Definition 2.4 (Maximum-Likelihood-Schatzer (MLE)) Die Zufallsvariablen Y1, . . . , Yn

seien unabhangig verteilt mit Wahrscheinlichkeitsfunktion bzw. Dichte f(yi | xi,θ), i = 1,. . . ,n.

Mit L(θ) =∏ni=1 f(yi|xi,θ) sei ihre Likelihood-Funktion, mit L(θ) = lnL(θ) ihre Log-Likelihood-

Funktion bezeichnet. Jede Maximalstelle

θML = maxL(θ); θ ∈ Θ = maxL(θ); θ ∈ Θ

heisst dann Maximum-Likelihood-Schatzer bzw. kurz ML-Schatzer fur θ. Der ML-Schatzer θML

muss nicht immer existieren und im Falle der Existenz nicht immer eindeutig sein.

Um das Maximum θML zu finden, ist es also notig, die Log-Likelihood-Funktion L(θ) nach θ zu

differenzieren und nullzusetzen. Die entstehenden Gleichungen sind im Allgemeinen nicht-linear

in θ, weshalb keine analytische Losung des Gleichungssystems moglich ist. In diesem Fall mussen

wir auf iterative Losungsverfahren zuruckgreifen. Passende Verfahren sind die Newton-Raphson-

Methode und das Fisher-Scoring-Verfahren. Um die Einfuhrung dieser Verfahren zu erleichtern,

werden zuerst noch einige notwendige Begriffe definiert.


Definition 2.5 (Score-Funktion) Die Score-Funktion bzw. der Score-Vektor s(θ) ist der Vek-

tor der Ableitungen der Log-Likelihood-Funktion nach θ:

s(θ) :=∂L(θ)∂θ

=n∑i=1

∂ ln fi∂θ

∈ Rp.

Die Maximum-Likelihood-Gleichungen lauten dann s(θ) = 0. Ihre Losung ist der ML-Schatzer.

Definition 2.6 (Fisher-Informationsmatrix) Die beobachtete Fisher-Informationsmatrix

wird durch

FIobs(θ) := −(∂2L(θ)∂θ∂θt

)definiert, wobei ∂θ∂θt die vektorwertige Differentiation abkurzt.

Die erwartete Fisher-Informationsmatrix dagegen ist durch die Kovarianzmatrix des Score-

Vektors definiert:

FI(θ) := Cov(s(θ)) = E(s(θ)s(θ)t).

Satz 2.4 (Fisher-Informationsmatrix) Unter Regularitatsbedingungen lasst sich die erwar-

tete Fisher-Informationsmatrix durch folgende Beziehung direkt aus der beobachteten Fisher-

Informationsmatrix berechnen:

FI(θ) = E(FIobs(θ)) = −E(∂2L(θ)∂θ∂θt

)

Beweis: Die erwartete Fisher-Informationsmatrix ist als Kovarianzmatrix des Score-Vektors

definiert. Deshalb ist zu zeigen, dass E(FIobs(θ)) = Cov(s(θ)) gilt.

Zunachst gilt:

∂2L(θ)∂θ∂θt

=∂2

∂θ∂θtlnL(θ) =

∂

∂θt

[∂∂θL(θ)

L(θ)

]=

∂2

∂θ∂θtL(θ)L(θ) − ∂∂θL(θ) ∂

∂θtL(θ)

(L(θ))2

=∂2

∂θ∂θtL(θ)

L(θ)−

∂∂θL(θ) ∂

∂θtL(θ)

L(θ)L(θ)︸︷︷︸∂

∂θlnL(θ) ∂

∂θt lnL(θ)

.

Daraus folgt:

E(FIobs(θ)) = −∫ [

∂2

∂θ∂θtlnL(θ)

]L(θ) dy

= −∫ [ ∂2

∂θ∂θtL(θ)

L(θ)− ∂

∂θlnL(θ)

∂

∂θtlnL(θ)

]L(θ) dy

= −∫

∂2

∂θ∂θtL(θ) dy +

∫∂

∂θlnL(θ)︸︷︷︸s(θ)

∂

∂θtlnL(θ)︸︷︷︸s(θ)t

L(θ) dy


(∗)= 0 +

∫s(θ)s(θ)tL(θ) dy

= E(s(θ)s(θ)t)

= Cov(s(θ)), da E(s(θ)) = 0.

Dabei gilt die Gleichheit in (*), da unter den Regularitatsbedingungen (siehe Ende diese Ab-

schnitts) beim ersten Integral Integration und Differentiation vertauscht werden durfen. Man

erhalt:

−∫

∂2

∂θ∂θtL(θ) dy = − ∂2

∂θ∂θt

∫L(θ) dy︸︷︷︸

=1, da L(θ)Dichte

= − ∂2

∂θ∂θt1 = 0

2.3.2 Iterative Algorithmen zur numerischen Bestimmung des ML-Schatzers

Das Newton-Raphson-Verfahren (allgemein)

Wir suchen das θ= (θ1, . . . , θr)t, das

f(θ) =

⎡⎢⎢⎣f1(θ1, . . . , θr)

...

fr(θ1, . . . , θr)

⎤⎥⎥⎦ =

⎡⎢⎢⎣

0...

0

⎤⎥⎥⎦

lost. Die Newton-Raphson-Methode liefert uns ein iteratives Verfahren, um diese Nullstellen-

gleichung approximativ zu losen. Es sei θ die Losung und θ0 nahe bei θ. Die Taylorentwicklung

erster Ordnung um θ0 lautet:

0 = f(θ) ≈ f(θ0) +Df(θ0)(θ − θ0) (2.5)

Dabei gilt:

Df(θ0) :=

⎡⎢⎢⎣

df1dθ1

· · · df1dθr

.... . .

...dfr

dθ1· · · dfr

dθr

⎤⎥⎥⎦

an der Stelle θ = θ0.

Die Approximation (2.5) kann nun unter der Voraussetzung, dass die Jacobimatrix Df(θ0) nicht

singular ist, nach θ aufgelost werden:

θ ≈ θ0 − [Df(θ0)]−1f(θ0)

Ausgehend vom Startwert θ0 berechne man nun iterativ

θi+1 = θi − [Df(θi)]−1f(θi), i = 0, 1, 2, . . . ;θ0 gegeben (2.6)


Diese Gleichung nennt sich Newton-Iteration. Zur Berechnung des neuen Naherungswertes benoti-

gen wir also f(θ) und die Jacobimatrix Df(θ) an der Stelle θi. Falls || θi+1 − θi || klein genug

ist, setzt man θ = θi+1. Dabei bezeichnet ||.|| die euklidische Norm.

Fur die Konvergenz des Newton-Verfahrens ist die Bedingung Df(θ ) = 0 notwendig, d.h. die

Jacobimatrix darf nicht singular sein, und die Bedingung | f(θ)D2f(θ)

Df(θ)2|≤ K < 1 (K = const.)

hinreichend. Diese Bedingungen mussen in einer Umgebung von θ, die alle Punkte θi sowie θ

enthalt, erfullt sein. Im Falle der Konvergenz des Newton-Verfahrens, verdoppelt sich bei je-

dem Iterationsschritt in etwa die Anzahl der genauen Stellen. Man spricht von quadratischer

Konvergenz (siehe [Bron 99]).

Das Newton-Raphson-Verfahren zur approximativen Berechnung des MLE in GLMs

Entsprechend Definition (2.5) lauten die Maximum-Likelihood-Gleichungen in GLMs s(β) = 0.

Dabei stellt die Score-Funktion s(β) den Vektor der Ableitungen der Log-Likelihood-Funktion

nach β dar. Um den MLE β zu erhalten, mussen wir also die Maximum-Likelihood- bzw. die

Score- Gleichungen s(β) = 0 losen. Der gesuchte Vektor β kann allerdings nicht immer analy-

tisch, d.h. durch direktes Umformen der Maximum-Likelihood-Gleichungen, ermittelt werden,

denn oft sind diese Gleichungen nicht-linear in β. In diesem Fall steht uns das Newton-Raphson-

Verfahren zur Verfugung.

Wir suchen also die Losung der Score-Gleichungen

s(β) =∂L(β)∂β

= 0.

Es sei β die Losung und β0 nahe bei β. Die zur Linearisierung der Scorefunktion verwendete

Taylorentwicklung erster Ordnung um β0 lautet:

0 = s(β) ≈ s(β0) + s′(β0)(β − β0) (2.7)

Dabei gilt nach Definition (2.6)

s′(β0) =∂2L(β0)∂β0∂β0

t = −FIobs(β0).

Die Approximation (2.7) kann nun nach β aufgelost werden:

β ≈ β0 + FI−1obs(β0)s(β0)

Zur Berechnung des MLE benotigen wir also die Inverse der beobachteten Fisher-Informationsmatrix

FIobs an der Stelle β0 und die Score-Funktion s an der Stelle β0.

Es sei nun β0 der Startwert. Im i-ten Iterationsschritt berechne man:

βi+1 = βi + FI−1obs(βi)s(βi) (2.8)


Das Produkt aus inverser beobachteter Fisher-Informationsmatrix an der Stelle βi und Score-

Vektor an der Stelle βi stellt dabei einen Korrekturfaktor dar. Im Laufe der Iterationen ver-

schwindet der Korrekturfaktor, weil die Score-Funktion gegen den Nullvektor strebt, und das

Verfahren bricht ab. Wenn also ||βi+1 − βi|| klein genug ist, setzt man β = βi+1. Eine andere

Moglichkeit ist, solange zu iterieren, bis ein bestimmtes Abstandskriterium unter einer Schranke

ε bleibt. Ein solches Kriterium ware z.B.:

||βi+1 − βi||||βi||

< ε , ε > 0

Wieder bezeichnet ||.|| die euklidische Norm.

Haben wir nun mit Hilfe des Newton-Raphson-Verfahrens einen Schatzer β errechnet, mussen

wir noch uberprufen, ob es sich dabei wirklich um ein Maximum handelt.

Dazu benotigen wir die Hesse-Matrix H(β):

H(β) :=

⎡⎢⎢⎣

∂s1(β)∂β1

· · · ∂s1(β)∂βp

.... . .

...∂sp(β)∂β1

· · · ∂sp(β)∂βp

⎤⎥⎥⎦ =

⎡⎢⎢⎣

∂2L(β)∂β1∂β1

· · · ∂2L(β)∂βp∂β1

.... . .

...∂2L(β)∂β1∂βp

· · · ∂2L(β)∂βp∂βp

⎤⎥⎥⎦ =

∂2L(β)∂β∂βt

Ist die Hesse-Matrix an der Stelle β = β negativ-definit, ist unser erhaltener Schatzer β wirklich

ein Maximum. Meist sind die berechneten Schatzer β allerdings keine globalen Maxima der Log-

Likelihood-Funktion. In der Regel stellen sie nur Losungen der Score-Gleichungen dar. Wenn die

Hesse-Matrix negativ-definit ist, entspricht dies also einem lokalen Maximum. Ist jedoch die Log-

Likelihood-Funktion konkav, stimmen das globale und das lokale Maximum uberein. Folgender

Satz sagt sogar, dass das Maximum fur strikt konkave Log-Likelihood-Funktionen eindeutig ist.

Satz 2.5 Sei V eine nicht-leere konkave Menge V ⊂ Rm und f : V → R eine strikt konkave

Funktion auf V, so gibt es hochstens einen Maximalpunkt von f auf V.

Beweis: Fur alle Paare x1,x2 ∈ V mit x1 = x2 und jede Zahl λ ∈ (0, 1) ist eine strikt konkave

Funktion definiert durch die Ungleichung

f(λx1 + (1 − λ)x2) > λf(x1) + (1 − λ)f(x2).

Nehmen wir nun an, es gabe zwei verschiedene Maximalpunkte x1 und x2 ∈ V von f , dann gilt

fur die Verbindungsgerade von x1 und x2

λf(x1) + (1 − λ)f(x2) = f(x1) = f(x2).

Nach Definition der Konkavitat gilt aber fur jedes λ ∈ (0, 1) und jedes x = λx1 +(1−λ)x2 ∈ V ,

das zwischen x1 und x2 liegt:

f(x) = f(λx1 + (1 − λ)x2) > λf(x1) + (1 − λ)f(x2) = f(x1) = f(x2)

Dies ist jedoch ein Widerspruch dazu, dass x1 und x2 Maxima sind.


Ausserdem gilt nach [Krab 83] (Seite 59) der folgende Satz:

Satz 2.6 Eine Funktion f : V → R auf einer nicht-leeren offenen und konkaven Teilmenge des

Rn, die fur jedes x∈ V stetige zweite partielle Ableitungen besitzt, ist genau dann strikt konkav

auf V, wenn ihre Hesse-Matrix H(x) fur jedes x∈ V negativ-definit ist.

Beispiel 2.4 (ML-Schatzung unter Poisson-Verteilungsannahme) Angenommen die

Yi, i = 1, . . . , n, sind Realisierungen Poisson-verteilter Zufallsvariablen gegeben die zugehorigen

Regressoren xi1, . . . , xip, i = 1, . . . , n. Die Designmatrix X = (xt1, . . . ,x

tn)t besitze vollen Spal-

tenrang p. Die kanonische Linkfunktion sei gegeben durch µi = E(Yi|xi) = exp(xtiβ).

Ausgehend von der Wahrscheinlichkeitsfunktion der Poissonverteilung fur unabhangige Zufalls-

variablen Yi, i = 1, . . . , n, lautet die Log-Likelihood-Funktion:

L(β) =n∑i=1

yiθi − b(θi)

a(φ)+ c(yi, φ)

=n∑i=1

−µi + yi lnµi − ln(yi!)

=n∑i=1

− exp(xtiβ) + yix

tiβ − ln(yi!)

Um den ML-Schatzer β von β = (β1, . . . , βp)t zu bestimmen, mussen wir die Log-Likelihood-

Funktion L(β) nach β ableiten und nullsetzen:

∂L(β)∂β

=n∑i=1

−xi exp(xtiβ) + yixi

=n∑i=1

yi − exp(xt

iβ)

xi

=

⎛⎜⎜⎝

∑ni=1

yi − exp(xt

iβxi1

...∑ni=1

yi − exp(xt

iβxin

⎞⎟⎟⎠ = 0

Wie zuvor schon erwahnt wurde, benotigen wir aber noch die Hesse-Matrix um zu uberprufen, ob

tatsachlich ein Maximum vorliegt. Dazu berechnen wir die einzelnen Matrixeintrage, die zweiten

partiellen Ableitungen.

∂2L∂βs∂βr

=∂

∂βs

n∑i=1

yi − exp(xt

iβ)xir

= −n∑i=1

exp(xtiβ)xirxis, r, s = 1, . . . , p.

Damit lautet die vollstandige Hesse-Matrix:

H(β) =∂2L(β)∂β∂βt

= −n∑i=1

exp(xtiβ)xix

ti


Es wurde bereits in [Siko 02](Seite 20f) gezeigt, dass die Hesse-Matrix fur beliebige β negativ-

definit ist (ztH(β)z < 0 ∀ z = 0). Die Log-Likelihood-Funktion ist damit nach Satz 2.6

global strikt konkav, und nach Satz 2.5 sind die Bedingungen fur ein eindeutiges sprich globales

Maximum erfullt. Der ML-Schatzer β stimmt als lokales Maximum der Log-Likelihood-Funktion

L(β) mit dem globalen Maximum uberein. Der Newton-Raphson-Algorithmus kann nun benutzt

werden, um dieses Maximum zu bestimmen.

IWLS-Algorithmus und das Fisher-Scoring-Verfahren (siehe [McCu 89], Seite 40ff)

Die Maximum-Likelihood-Schatzer der Regressionsparameter β im linearen Pradiktor ηi = xitβ

konnen durch iterative gewichtete kleinste-Quadrate-Schatzung (IWLS) bestimmt werden. In

dieser Regression ist die abhangige Variable nicht y sondern z, eine linearisierte Form der Link-

funktion angewandt auf y. Die Gewichte sind Funktionen der gefitteten Werte µ. Der Prozess

ist iterativ, da sowohl die adjusted dependent Variablen z = (z1, . . . , zn)t als auch die Gewichte

wi, i = 1, . . . , n, von den gefitteten Werten µ abhangen, fur die aber nur momentane Schatzer

verfugbar sind. Der IWLS-Algorithmus kann wie folgt beschrieben werden:

IWLS-Algorithmus

Schritt 1:

Es sei b0 der momentane Schatzer des MLE β. Man bilde

η0i := xi

tb0 und µ0i = g−1(η0

i ) , i = 1, . . . , n.

Die adjusted dependent Variablen z0 = (z01 , . . . , z

0n)t werden durch

z0i := η0

i + (yi − µ0i )(dηidµi

)|ηi=η0i

definiert, wahrend die Gewichte durch

w0i :=

[dηidµi

|ηi=η0i

]−2 1v0i

bestimmt werden. Dabei bezeichnet v0i die Varianzfunktion an der Stelle µ0

i .

Schritt2:

Nun regressiere man z0i auf die Kovariablen xi1, . . . , xip mit Gewichtung w0

i , um neue Parame-

terschatzer b1 zu erhalten und gehe zu Schritt 1. Dieses Vorgehen wiederhole man so lange, bis

||b0 − b1|| klein genug ist. Dabei bezeichnet ||.|| die euklidische Norm.

Diese Anpassungsprozedur gilt es nun zu rechtfertigen. Dazu werden zunachst die Maximum-

Likelihood-Gleichungen fur die βj , j = 1, . . . , p, in einem GLM hergeleitet.


Der Log-Likelihood fur die Beobachtung yi lautet:

Li(β) = L(yi, µi, φ) =yiθi − b(θi)

a(φ)+ c(yi, φ),

wobei φ bekannt, g(µi) = ηi, µi = E(Yi) = h(θi) , θi der kanonische Parameter, ηi = xtiβ und

xi = (xi1, . . . , xip)t ist. Der Log-Likelihood aller Beobachtungen ist dann gegeben durch

L(β) =n∑i=1

L(yi, µi, φ). (2.9)

Zur Bestimmung des MLE sind die Score-Gleichungen

sj(β) =∂L(β)∂βj

=n∑i=1

L(yi, µi, φ)∂βj

= 0, j = 1, . . . , p (2.10)

zu losen. Wir betrachten zunachst

∂Li(β)∂βj

=∂Li(β)∂µi

∂µi∂ηi

∂ηi∂βj

, j = 1, . . . , p.

Da

ηi =p∑j=1

βjxij, gilt∂ηi∂βj

= xij , j = 1, . . . , p.

Ferner gilt

∂Li(β)∂θi

=∂Li(β)∂µi

dµidθi

⇒ ∂Li(β)∂µi

=∂Li(β)∂θi

dµidθi

=yi−b′(θi)a(φ)

dµidθi

.

Da µi = b′(θi) und b′′(θi) = vi, gilt dµi/dθi = vi. Dementsprechend folgt:

∂Li(β)∂µi

=yi−b′(θi)a(φ)

vi

Damit ergibt sich als Score-Gleichung in einem GLM fur die Beobachtung yi

∂Li(β)∂βj

=yi − b′(θi)a(φ)

1vi

dµidηi

xij =wia(φ)

(yi − µi)dηidµi

xij = 0

mit

wi := v−1i g′(µi)−2 = v−1

i

(dµidηi

)2

.

Die Score-Gleichungen in einem GLM fur n Beobachtung ergeben sich somit als:

sj(β) =∂L(y,β)∂βj

=n∑i=1

∂L(yi, µi, φ)∂βj

=n∑i=1

(yi − µi)a(φ)

widηidµi

xij = 0, j = 1, . . . , p (2.11)

Fur konstanten Dispersionsparameter (a(φ) = const.) verschwindet der Faktor a(φ). Wie wir be-

reits wissen, sind die Score-Gleichungen allerdings im Allgemeinen nicht-linear in β und konnen


somit nur approximativ numerisch gelost werden. Fur die Anwendung der Newton-Raphson-

Methode zum Losen der Score-Gleichungen s(β) = 0 fehlt uns allerdings noch die Hesse-Matrix

H(β) = ∂2L∂βs∂βr

, r, s = 1, . . . , p. Ein Eintrag der Hessematrix ergibt sich mit Hilfe der Produkt-

regel folgendermaßen:

H(β) =∂2L

∂βs∂βr=

∂

∂βs

[n∑i=1

1a(φ)

(yi − µi)widηidµi

xir

]

=n∑i=1

1a(φ)

(yi − µi)∂

∂βs

[wi∂ηi∂µi

xir

]+

n∑i=1

1a(φ)

widηidµi

xir

∂

∂βs(yi − µi)︸︷︷︸

− dµidηi

∂ηi∂βs

dηi∂βs

=xis

=n∑i=1

1a(φ)

(yi − µi)∂

∂βs

[wi∂ηi∂µi

xir

]−

n∑i=1

1a(φ)

wixirxis (2.12)

Da im Allgemeinen ∂2L∂βs∂βr

von den Daten abhangt, modifizieren wir die Newton-Raphson-

Methode und benutzen an Stelle der beobachteten Informationsmatrix die erwartete. Bei der

Bildung des Erwartungswertes verschwindet der erste Term in Gleichung (2.12). Insbesondere

gilt

−E(

∂2L∂βs∂βr

)= 0 − E

(−

n∑i=1

1a(φ)

wixirxis

)=

1a(φ)

n∑i=1

wixirxis.

Damit lasst sich die erwartete Fisher-Informationsmatrix wie folgt darstellen:

FI := −(E

(∂2L

∂βs∂βr

))s,r=1,...,p

=

(1

a(φ)

n∑i=1

wixirxis

)s,r=1,...,p

=1

a(φ)XtWX (2.13)

Dabei gilt:

W = diag(w1, . . . , wn) ∈ Rn×n mit wi =

1

vi

(dηidµi

)2 =(dµidηi

)2

v−1i

Der Ubergang von der beobachteten zur erwarteten Fisher-Informationsmatrix verlangt nach

einem neuen Algorithmus. Das Fisher-Scoring-Verfahren benutzt den Score-Vektor

s(β) = ∂L∂β

und die erwartete Fisher-Informationsmatrix FI. Ausgehend vom momentanen

Schatzer b0 fur den MLE β leiten wir eine Anpassung δb0 her, definiert als Losung von

FI δb0 =∂L∂b0

= s(b0).

Der neue Schatzer b1 = b0 + δb0 von β erfullt somit die Gleichung

FI b1 = FI b0 + FI δb0 = FI b0 + s(b0).


Wie aus Gleichung (2.11) ersichtlich, gilt fur die Komponenten von s(b0) unter Vernachlassigung

des Dispersionsparameters

sr(b0) =n∑i=1

w0i (yi − µ0

i )dηidµi

|ηi=η0ixir, r = 1, . . . , p, (2.14)

und entsprechend Gleichung (2.13)

FIrs =n∑i=1

w0i xirxis, s, r = 1, . . . , p. (2.15)

Es folgt:

(FIb0)r =p∑s=1

FIrs b0s =

p∑s=1

n∑i=1

w0i xirxisb

0s (2.16)

Der neue Schatzer b1 erfullt somit:

(FIb1)r = (FI b0)r + sr(b0) =n∑i=1

w0i xir

[p∑s=1

xisb0s + (yi − µ0

i )dηidµi

|ηi=η0i

]

=n∑i=1

w0i xir

[η0i + (yi − µ0

i )dηidµi

|ηi=η0i

](2.17)

Andererseits gilt:

(FIb1)r =p∑s=1

FIrsb1s

(2.15)=

p∑s=1

(n∑i=1

w0i xirxis

)b1s =

n∑i=1

w0i xir

p∑s=1

xisb1s =

n∑i=1

w0i xirη

1i (2.18)

Mit Gleichung (2.17) und Gleichung (2.18) folgt:

n∑i=1

w0i xir

[η0i + (yi − µ0

i )dηidµi

|ηi=η0i

]=

n∑i=1

w0i xirη

1i (2.19)

Diese Gleichungen entsprechen den Normalengleichungen einer gewichteten Kleinste-Quadrate-

Schatzung (Weighted Least Squares) mit Gewichten w0i = 1

v0i

[dηi

dµi|ηi=η0i

]−2, den Kovariablen-

vektoren x1, . . . , xp und der abhangigen Variable z0i := η0

i + (yi − µ0i )dηidµi

|ηi=η0i, i = 1, . . . , n.

Dabei bezeichnet man die Variablen z als adjusted independent variables.

Es sei noch angemerkt, dass fur kanonische Linkfunktionen eine Vereinfachung auftritt. In diesem

Fall stimmen namlich der erwartete und der beobachtete Wert der Hessematrix (siehe Gleichung

(2.12)) uberein, so dass die Fisher-Scoring-Methode und das Newton-Raphson-Verfahren densel-

ben Algorithums darstellen. Das liegt daran, dass die lineare Gewichtsfunktion W dη/dµ in den

Score-Gleichungen (2.11) in diesem Fall eine Konstante ist, wodurch der erste Term in (2.12)

identisch gleich Null ist.


Da fur den ML-Schatzer β in einem GLM im Allgemeinen keine analytischen Losungen exi-

stieren, konnen nur asymptotische Eigenschaften gefolgert werden. Diese Eigenschaften sind

allerdings nur unter bestimmten Voraussetzungen, den sogenannten Regularitatsbedingungen,

gegeben. Diese lassen sich nach [Came 98] wie folgt zusammenfassen:

(i) Die Dichte f(y|xi,β) ist global definiert, und es gilt f(y|xi,β1) = f(y|xi,β1) fur alle

β1 = β2.

(ii) Es gilt β ∈ Θ, wobei der Parameterraum Θ endlichdimensional, abgeschlossen und kom-

pakt ist.

(iii) Es existieren stetige und beschrankte Ableitungen von L(β) bis zur 3. Ordnung.

(iv) Die Reihenfolge von Differentiation und Integration der Likelihood-Funktion darf ver-

tauscht werden.

(v) Die Kovariablenvektoren xi, i = 1, . . . , n erfullen folgende Bedingungen:

• xtixi <∞

• E(ω2i )∑

i E(ω2i )

= 0 ∀i, wobei ωi ≡ xit ∂ ln f(yi|xi,β)

∂β

• limn→∞∑n

i=1 E(ω2i |Ωi−1)∑n

i=1 E(ω2i )

= 1, wobei Ωi−1 = (x1, x2, . . . , xi−1).

Die erste Bedingung stellt sicher, dass ein eindeutiges Maximum existiert. Die zweite Voraus-

setzung verhindert mogliche Probleme am Rand des Parameterraums Θ und kann, wenn die

Log-Likelihood-Funktion L z.B. global konkav ist, weggelassen werden. Die dritte Bedingung

kann oft entscharft werden, indem man die Existenz der Ableitungen nur bis zur Ordnung 2

fordert. Die vierte Bedingung schließt solche Dichten aus, bei denen der Wertebereich von yi

von θ abhangt. Die letzte Bedingung sorgt dafur, dass jede Beobachtung ausgeschlossen wird,

deren Anteil an der Likelihood-Funktion zu groß ist.

Eine der wichtigsten Eigenschaften des Maximum-Likelihood-Schatzers β in GLMs ist nach

[Fahr 94] die asymptotische Normalverteilung:

Die Verteilung des ML-Schatzers ist eine Normalverteilung fur n→ ∞, formal

βa∼ N(β, F I−1(β)). (2.20)

Dies ist gleichbedeutend mit

F 1/2(β)(β − β) d→ N(0, Ip). (2.21)

Dabei steht das d fur Konvergenz in Verteilung, und Ip bezeichnet eine p-dimensionale Einheits-

matrix.

2.4. Asymptotische Hypothesentests in GLMs 21

Der Parametervektor β ist somit asymptotisch Normal-verteilt, und die asymptotische Kovari-

anzmatrix von β kann durch

Cov(β)a≈ FI−1(β) (2.22)

approximiert werden. Dabei stellt FI−1(β) die Inverse der Fisher-Informationsmatrix ausgewer-

tet an der Stelle β = β dar.

2.4 Asymptotische Hypothesentests in GLMs

Wurden durch explorative Datenanalyse bereits erste geeignete und vor allem sinnvolle Modelle

vorgeschlagen, gilt es, diese Modelle mittels statistischer Hypothesentests zu uberprufen. Deshalb

werden im Folgenden die wichtigsten Hypothesentests in GLMs vorgestellt.

Meist handelt es sich bei Testproblemen fur den Parametervektor β um lineare Hypothesen der

Form

H : Cβ = ξ gegen K : Cβ = ξ. (2.23)

Dabei muss die Matrix C vollen Zeilenrang q ≤ p + 1 haben. Uns interessiert insbesondere der

wichtige Spezialfall

H : βr = 0 gegen K : βr = 0. (2.24)

Man bezeichnet dabei mit βr einen Teilvektor von β. Der Hypothesentest (2.24) dient also zur

Uberprufung der Signifikanz der zu βr gehorenden Effekte. Es wird das durch βr = 0 definierte

Teilmodell mit dem betrachteten vollen Modell verglichen.

Zum Testen der beschriebenen Hypothesen verwendet man meist die Likelihood-Quotienten-

Statistik (Likelihood-Ratio-Statistic) oder die Wald-Statistik.

(i) Likelihood-Quotienten-Teststatistik:

Die Likelihood-Quotienten-Statistik lautet

lq = −2L(β) − L(β). (2.25)

Sie vergleicht das unrestringierte Maximum L(β) des Log-Likelihoods mit dem Maximum

L(β). Dabei bezeichnet β den restringierten ML-Schatzer unter der Nullhypothese

H : Cβ = ξ.

(ii) Wald-Teststatistik:

Die Wald-Statistik lautet

w = (Cβ − ξ)t[CFI−1(β)Ct

]−1(Cβ − ξ) (2.26)


Diese Statistik misst die durch die asymptotische Kovarianzmatrix CFI−1(β)Ct von Cβ

gewichtete Distanz zwischen dem unrestringierten Schatzer Cβ von Cβ und dem hypo-

thetischen Wert ξ = Cβ unter der Nullhypothese.

Diese Teststatistiken sind nach [Fahr 94] (Seite 46) unter der Nullhypothese von (2.23) asympto-

tisch aquivalent und besitzen dieselbe χ2-Grenzverteilung mit q Freiheitsgraden, d.h. lq, w a∼ χ2q.

Analog weisen die beiden Teststatistiken unter der Hypothese von (2.24) die χ2-Grenzverteilung

mit r Freiheitsgraden auf. Einen ausfuhrlichen Beweis fur die asymptotische χ2-Verteilung der

Teststatistiken findet man in [Fran 01] (Seite 24f). Fur weitergehende Erlauterungen lese man

[Rao 73] (Seite 415ff). Dort werden die Hypothesentests fur den allgemeinen Fall beschrieben,

d.h. ohne die Annahme, dass die betrachteten Modelle zur Klasse der GLMs gehoren. Damit

lassen sich nun die Ablehnungsbereiche der betrachteten Tests angeben.

(i) Likelihood-Quotienten-Test

Die Hypothese H wird zum Signifikanzniveau α > 0 genau dann zugunsten der Alternative

K abgelehnt, wenn gilt:

lq > χ2q,1−α (2.27)

Dabei bezeichnet χ2q,1−α das (1 − α)-Quantil der χ2-Verteilung mit q Freiheitsgraden. Be-

nutzt man den Likelihood-Ratio-Test fur das Testproblem (2.24) erfordert dies die ML-

Schatzung im entsprechenden Teilmodell (βr = 0) zusatzlich zur Schatzung von β im

vollen Modell und damit weitere Iterationen im Fisher-Scoring-Verfahren. Fur die Berech-

nung des ML-Schatzers β im Testproblem (2.23) allerdings ist der numerische Aufwand

wegen der allgemeineren linearen Gleichungsrestriktionen noch erheblich großer.

(ii) Wald-Test (t-Test):

Die Hypothese H wird zum Signifikanzniveau α > 0 genau dann zugunsten der Alternative

K abgelehnt, wenn gilt:

w > χ2q,1−α (2.28)

Betrachten wir das Testproblem (2.24) und wollen nicht die Signifikanz mehrerer Effekte

uberprufen, sondern interessieren uns nur fur die Signifikanz des i-ten Effekts, lautet die

lineare Hypothese

H : βi = 0 gegen K : βi = 0 (i fest), (2.29)

und die Wald-Statistik vereinfacht sich zu

w =β2i

FI−1ii

. (2.30)

2.5. Goodness-of-Fit-Maße 23

Wir verwerfen die Hypothese H : βi = 0 gegen die Alternative K : βi = 0 zum Niveau

α > 0 genau dann, wenn gilt:

w > χ21,1−α (2.31)

In diesem Fall stimmt die Wald-Statistik mit dem Quadrat des sogenannten t-Wertes

ti =βi√FI−1

ii

, (2.32)

dem standardisierten Schatzer von βi, uberein. Mit FI−1ii sei dabei jeweils das i-te Dia-

gonalelement der geschatzten asymptotischen Kovarianzmatrix FI(β)−1 bezeichnet. Da

der Wald-Test die ML-Schatzer des entsprechenden Teilmodells nicht benotigt, ist der

numerische Aufwand im Vergleich zum Likelihood-Quotienten-Test wesentlich geringer.

2.5 Goodness-of-Fit-Maße

Nachdem wir nun wissen, wie wir die Parameter unseres Modells schatzen, benotigen wir ein

Maß, das uns sagt, wie gut die gefundenen Schatzer die Datenstruktur anpassen. Dazu mussen

wir wieder auf die Log-Likelihoods zuruckgreifen.

Liegen n Beobachtungen vor, konnen Modelle mit bis zu n Parametern angepasst werden. Im

einfachsten Modell, dem Nullmodell, mit nur einem Parameter β1 suchen wir den Schatzer β1, der

den gemeinsamen Schatzer fur alle Beobachtungen darstellt. Damit ordnet das Nullmodell die ge-

samte Variation zwischen den Beobachtungen der zufalligen Komponente zu. Im Gegensatz zum

Nullmodell steht das volle Modell mit n Parametern. Hier setzt man einfach yi = yi, i = 1, . . . , n.

Dieses Modell ordnet die gesamte Variation zwischen den Beobachtungen der systematischen

Komponente zu. In Bezug auf die praktische Anwendung ist das Nullmodell zu einfach und

das volle Modell nicht aussagekraftig genug, denn es fasst die Daten nicht in ihrer wesentlichen

Struktur zusammen, sondern gibt sie vollstandig wieder. Jedoch hilft uns das volle Modell, die

Diskrepanz zu einem Zwischenmodell mit p < n Parametern zu messen.

2.5.1 Devianzanalyse

Im Folgenden sei der Skalenparameter φ bekannt oder zumindest fest gewahlt. Den Log-Likelihood

des zu untersuchenden Modells mit den p < n Parametern bezeichnen wir mit L(µ,y), den Log-

Likelihood des vollen Modells mit L(y,y). Im vollen Modell ist y derjenige Schatzer, der den

Log-Likelihood maximiert. Damit konnen wir nun die Devianz als Goodness-of-Fit-Mass defi-

nieren, um die Abweichung zwischen den wahren und den gefitteten Werten zu messen.

Definition 2.7 (Devianz) Die skalierte Devianz D∗ ist definiert durch

D∗(y, µ, φ) = −2[L(µ,y) − L(y,y)].


Somit ist die skalierte Devianz eine Funktion von y und µ , die durch a(φ) geteilt wird. Wir

multiplizieren die skalierte Devianz D∗ mit a(φ) und erhalten die Devianz D:

D(y, µ) = −2a(φ)[L(µ,y) − L(y,y)]

Beispiel 2.5 (Die Devianz im Poissonmodell) Der Log-Likelihood im Poissonmodell wurde

bereits in Beispiel 2.4 bestimmt. Aus Beispiel 2.2 wissen wir, dass fur die Poissonverteilung

a(φ) = 1 gilt. Damit stimmen im Poissonmodell die skalierte Devianz D∗ und die Devianz D

uberein. Wir erhalten also:

DPoi(y, µ) = D∗(y, µ) = −2[L(µ,y) − L(y,y)]

= −2n∑i=1

(yi ln µi − µi − ln(yi!)) − (yi ln yi − yi − ln(yi!))

= 2

n∑i=1

yi ln(yi/µi) − (yi − µi)

Falls yi = 0 gilt, setze man yi ln yi = 0.

In der Definition der Devianz hangt L(y,y) nicht von den Parametern ab. Das heisst, die Maxi-

mierung von L(µ,y) ist aquivalent zur Minimierung von D∗(y, µ, φ) bezuglich µ. Die Konstante

2 dient lediglich der Normalisierung. Durch die 2 als Vorfaktor ist sichergestellt, dass im Falle

der Normalverteilung die Devianz mit der Residuenquadratsumme ubereinstimmt. Folglich hat

die Devianz bei Modellen mit Normalverteilung eine χ2 -Verteilung mit n − p Freiheitsgraden.

Es lasst sich zeigen, dass die Devianz unter anderen Verteilungen asymptotisch χ2n−p -verteilt

ist (siehe [Fahr 94], Seite 48).

Mit Hilfe der Devianz wollen wir zwei genestete1 Modelle vergleichen. Es interessiert uns, ob

die Hinzunahme eines Regressors die Anpassung eines Modells signifikant verbessert. Wir be-

trachten die Reduktion in der Devianz, die durch die Hinzunahme des neuen Regressors bedingt

ist:

D(y, µ0) −D(y, µΛ) = 2[L(µΛ) − L(µ0)] (2.33)

Dabei bezeichnet µ0 die angepassten Werte des urspruglichen Modells und µΛ die angepassten

Werte des erweiterten Modells mit dem zusatzlichen Parameter. Die Statistik (2.33) ist dann

asymptotisch χ21 -verteilt.

Mit der Devianz und ihrer asymptotischen χ2-Verteilung kann der sogenannte Residual-Deviance-

Test durchgefuhrt werden, um die Anpassung des gefundenen Modell zu uberprufen.

1Ein Modell Mi ist in Modell Ms genestet, wenn alle Kovariablen von Modell Mi im Modell Ms, s ≥ i + 1,

enthalten sind. M1 ⊂ M2 ⊂ . . . ⊂ Mr. M1 = Null-Modell, Mr =”saturated“ Modell


Man testet die allgemeine Hypothese

H: Modell passt gegen K: Modell passt nicht.

Die Nullhypothese verwirft man genau dann, wenn

Devianz > χ2n−p,1−α

gilt, wobei n die Anzahl der Beobachtungen, p die Anzahl der zu schatzenden Parameter, α das

Signifikanzniveau und χ2n−p,1−α das (1−α)-Quantil der χ2-Verteilung mit n− p Freiheitsgraden

darstellt. Muss man die Nullhypothese verwerfen, so bedeutet das, die Anpassung des erarbei-

teten Modells an die Datenstruktur ist noch nicht ausreichend.

Der Residual-Deviance-Test uberpruft die Anpassungsgute des gesamten Modells. Im Gegensatz

dazu misst der Partial-Deviance-Test die Signifikanz von nur einer oder mehreren Kovariablen.

Die zu vergleichenden Modelle M1 und M0 mussen dabei die Eigenschaft haben, dass Modell

M1 im Vergleich zu Modell M0 genau diejenigen Parameter zusatzlich enthalt, deren Einfluss

getestet werden soll. Es handelt sich also um genestete Modelle. Der Regressionsparametervek-

tor β wird aufgespalten in β = (β0,β1)t, wobei Modell M0 aus den Parametern besteht, die in

β0 enthalten sind, und Modell M1 zusatzlich die Parameter aus β1 aufweist. Der entsprechende

Test lautet

H: β1 = 0 gegen K: β1 = 0.

Die Nullhypothese wird verworfen, wenn

Devianz(M0) - Devianz(M1) > χ2p1−p0,1−α

gilt. Dabei bezeichnet p1 die Anzahl der Parameter im Modell M1 und p0 die Anzahl der Pa-

rameter im Modell M0. Wieder steht χ2p1−p0,1−α fur das (1 − α)-Quantil der χ2

p1−p0-Verteilung.

Wird die Nullhypothese abgelehnt, so gilt dies als Nachweis fur einen signifikanten Einfluss der

zu den Parametern aus β1 gehorenden Kovariablen.

Bezeichnet man mit V das volle Modell, ergibt sich beim Partial-Deviance-Test

Devianz(M0) − Devianz(M1) = −2a(φ)(L(M0) − L(V )) + 2a(φ)(L(M1) − L(V ))

= −2a(φ)(L(M0) − L(M1)).

Dies entspricht bis auf den Vorfaktor a(φ) der Testgroße des Likelihood-Quotienten-Tests, wenn

M0 das Modell unter der Nullhypothese und M1 das unrestringierte Modell abkurzt.

Ein weiteres Anpassungsmaß ist die Pearson-χ2 -Statistik.


Definition 2.8 (Pearson-Statistik) Die Pearson-Statistik fur unabhangige Beobachtungen lau-

tet:

P =n∑i=1

(yi − µi)2

V (µi)

Dabei ist µi Schatzer von µi und V (µi) die geschatzte Varianzfunktion.

Solange der Erwartungswert und die Varianz bekannt sind, eignet sich dieses Goodness-of-Fit-

Maß fur beliebige Regressionsmodelle. Die Devianz hingegen beschrankt sich auf parametrische

Regressionsmodelle, da der Likelihood definiert sein muss. Sind der Erwartungswert und die

Varianz korrekt spezifiziert, gilt:

E

(n∑i=1

(Yi − µi)2

V arYi

)= n , da E

((Yi − µi)2

V arYi

)= 1

Die Pearson-Statistik P besitzt bei Modellen mit Normalverteilung eine exakte χ2n−p -Verteilung,

da sich die Pearson-Statistik in diesem Fall aus der Summe quadrierter standardisierter Normal-

verteilter Zufallsvariablen ergibt. Wie die Devianz ist unter anderen Verteilungen auch die

Pearson-Statistik zumindest asymptotisch χ2-verteilt mit n−p Freiheitsgraden (siehe [Fahr 94],

Seite 48)).

Beispiel 2.6 (Pearson χ2 im Poissonmodell) Aus Beispiel 2.2 wissen wir, dass im

Poissonmodell V arYi = b′′(θi)a(φ) = µi gilt. Die Pearson-Statistik lautet somit:

PPoi =n∑i=1

(yi − µi)2

µi

Der Quotient auf der rechten Seite vergleicht die empirische Varianz mit dem empirischen Er-

wartungswert. Deshalb wird PPoi oft als Indikator fur die im Poissonmodell geforderte Gleich-

heit von Erwartungswert und Varianz benutzt. Der Dispersionsparameter φ lasst sich nach

[McCu 89](Seite 200) durch PPoin−p schatzen. Im Falle der Poissonverteilung ist φ = 1, wie wir

bereits aus Beispiel 2.2 wissen. Gilt also PPoi > n − p, schließen wir daraus, dass die Varianz

im Modell großer ist als der Erwartungswert. Gilt hingegen PPoi < n−p, interpretieren wir dies

als Hinweis dafur, dass die Varianz kleiner ist als der Erwartungswert.

Jedoch eignet sich die Pearson-Statistik nur beschrankt fur diese Diagnosen. Verallgemeinern

wir die Poissonverteilung, indem wir den Skalenparameter φ frei wahlen und setzen a(φ) = φ,

dann ist wegen Gleichung (2.4) und Beispiel 2.2 die Varianz ein Vielfaches des Erwartungswerts:

V arYi = φµi. Wahlen wir nun als Schatzer fur die Varianz V (µi) = φµi mit

φ = 1n−p

∑ni=1(yi−µi)2/µi, so ist die Pearson-Statistik immer gleich n−p und damit unabhangig

von φ. Sie lasst keine Ruckschlusse bzgl. der Gleichheit von Varianz und Erwartungswert mehr

zu.


Zur Beurteilung der Anpassungsgute eines oder mehrerer Modelle sollten allerdings nie nur die

Devianz und die Pearson-Statistik verwendet werden. Man benutzt diese Gutemaße, um eine

Vorentscheidung zu treffen, ob ein Modell geeignet ist oder nicht. Eine endgultige Aussage uber

die Anpassung des Modells erfordert eine detaillierte Residualanalyse zur genaueren Uberprufung

der Modellannahmen.

2.5.2 Residualanalyse

Residuen sind im klassischen Modell als Differenz zwischen dem beobachteten und dem ange-

passten Wert definiert:

ri = (yi − µi), i = 1, . . . , n

Dabei ist der angepasste Erwartungswert µi der bedingte Erwartungswert µi = µ(xtiβ) ausgewer-

tet an der Stelle β = β . Die Residuen dienen dazu, Modellmissspezifizierungen, Ausreisser oder

Beobachtungen mit schlechter Anpassung aufzuspuren. Mit Hilfe der Residuenanalyse konnen

moglicherweise die Art der Missspezifizierung bestimmt und Korrekturmoglichkeiten gefunden

werden. Auch zur Beurteilung des Ausmaßes der Missspezifizierung kann eine Residuenanalyse

hilfreich sein. In der klassischen linearen Regression mit normalverteilten, homoskedastischen

Fehlern gilt:

(yi − µi) ∼ N(0, σ2), i = 1, . . . , n

In großen Stichproben sind also die Residuen symmetrisch um Null verteilt mit konstanter Va-

rianz σ2. In dieser Diplomarbeit sollen allerdings Zahldaten untersucht werden. Fur diese ist

(yi − µi) heteroskedastisch und asymmetrisch. Es gibt kein Residuum, das einen Erwartungs-

wert von Null, konstante Varianz und symmetrische Verteilung gleichzeitig mitbringt. Um trotz-

dem einzelne Eigenschaften davon zu erhalten, stehen uns verschiedene Residuendefinitionen zur

Verfugung.

(i) Pearson-Residuum:

Das Pearson-Residuum rPi ist wie folgt definiert:

rPi =yi − µi√V (µi)

Hier wird zur Korrektur der Heteroskedastizitat das klassische Residuum mit der geschatz-

ten Standardabweichung skaliert. Somit hat das Pearson-Residuum in großen Stichproben

den Erwartungswert 0, ist homoskedastisch, allerdings asymmetrisch verteilt. Es gilt fol-

gende Beziehung zwischen der Pearson-Statistik und dem Pearsonresiduum, falls eine

Poissonverteilung vorliegt:n∑i=1

(rPi )2 = PPoi


(ii) Devianz-Residuum:

Verwendet man die Devianz als Anpassungsmaß in einem GLM, tragt jede Einheit einen

bestimmten Teil di zur Gesamtmasse bei:

n∑i=1

di = 2L(yi) − L(µi) = D

Das Devianz-Residuum lasst sich damit wie folgt definieren:

rDi = sign(yi − µi)√di mit sign(x) =

⎧⎪⎪⎨⎪⎪⎩

−1 fur x < 0

0 fur x = 0

1 fur x > 0

Dieses Maß wachst mit yi − µi, und es gilt∑n

i=1(rDi )2 = D . Fur den Fall der Poissonver-

teilung lassen sich die Residuen durch

rDi = sign(yi − µi)[2yi ln

yiµi

− (yi − µi)] 1

2

berechnen.

2.6 Die Negativ-Binomial-Verteilung (NB-Verteilung)

Die Modellierung von Zahldaten mit der Poissonverteilung stoßt leicht an ihre Grenzen. Die von

der Poissonverteilung geforderte Gleichheit von Erwartungswert und Varianz (Aquidispersion)

wird meist durch die Daten verletzt, da sie eine hohere Variabilitat als vom Poissonmodell erwar-

tet aufweisen. Diese Uberdispersion kann verschiedene Ursachen haben. Fehlende Regressoren

oder eine falsch gewahlte Linkfunktion zerstoren die Aquidispersion. Im Falle der Poissonver-

teilung kommt es aber insbesondere zu Problemen, wenn die zugrundeliegende Beobachtungs-

einheit wie z.B die Zeit oder das Volumen nicht fest sondern variabel ist. Um das Problem der

Uberdispersion in den Griff zu bekommen, bietet es sich an, eine Verteilung zu wahlen, die auf

die Forderung nach Aquidispersion verzichtet. Eine Moglichkeit, Uberdispersion zu modellieren,

bietet die Negativ-Binomial-Verteilung, deren Varianz eine quadratische Funktion des Erwar-

tungswertes ist und die damit bei der Modellierung der Varianz wesentlich mehr Flexibilitat

zulasst.

2.6.1 Definition und Eigenschaften der Negativ-Binomial-Verteilung

Die NB-Verteilung besitzt im Gegensatz zur Poissonverteilung einen zusatzlichen und damit

zwei Parameter. Zum besseren Verstandnis und zur Veranschaulichung wird die NB-Verteilung

vorerst als diskrete Verteilung mit ganzzahligem Parameter r eingefuhrt. Spater gehen wir dann

zur Erweiterung mit stetigem Parameter r uber.

2.6. Die Negativ-Binomial-Verteilung 29

Definition 2.9 (NB-Verteilung mit r ∈ N) Die Zufallsvariable Y heisst Negativ-Binomial-

verteilt mit den Parametern r und p, kurz Y ∼ NB(r, p), falls sich ihre Wahrscheinlichkeits-

funktion wie folgt darstellen lasst:

P (Y = y) =(y + r − 1

y

)pr(1 − p)y, y = 0, 1, 2, . . . , r ∈ N; 0 < p < 1.

Fur eine NB(r,p)-verteilte Zufallsvariable Y bezeichnet Y=y die Anzahl der Misserfolge vor dem

r-ten Erfolg in einer Reihe von unabhangigen Bernoulli-Versuchen mit Erfolgswahrscheinlichkeit

p.

Beispiel 2.7 Eine Straßenbahn verspatet sich bei ihrer Ankunft am Bahnhof in 5 von 6 Fallen.

Die Wahrscheinlichkeit fur zehn Verspatungen vor der dritten punktlichen Ankunft ist dann

gegeben durch:

P (Y = 10) =(

10 + 3 − 110

)(16

)3 (1 − 1

6

)1

≈ 0.25

Man sagt: ”Y ist NB(3, 16)-verteilt“.

Die Abbildungen 2.2 und 2.3 zeigen die Wahrscheinlichkeitsfunktionen von NB-Verteilungen

mit unterschiedlich gewahlten Parametern. In Abbildung 2.2 wurde die Anzahl der Erfolge auf

15 festgelegt und die Erfolgswahrscheinlichkeit variiert. Je kleiner die Erfolgswahrscheinlichkeit

p, desto weiter ist die Funktion nach rechts verschoben. Das heisst, es werden mehr Versuche

benotigt, bis wieder ein Erfolg eintritt.

0 10 20 30 40 50 60 70 80 90 1000

0.02

0.04

0.06

0.08

0.1

0.12NB(15,0.2)NB(15,0.4)NB(15,0.6)

Abbildung 2.2: NB-Verteilungsdichten fur verschiedene Erfolgswahrscheinlichkeiten p


In Abbildung 2.3 hingegen wurde die Erfolgswahrscheinlichkeit festgehalten und die Anzahl der

gewunschten Erfolge unterschiedlich gewahlt. Wie erwartet benotigt man bei gleicher Erfolgs-

wahrscheinlichkeit mehr Versuche, wenn man mehr Erfolge erzielen will.

0 10 20 30 40 50 60 70 80 90 1000

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

0.045NB(15,0.2)NB(10,0.2)NB(5,0.2)

Abbildung 2.3: NB-Verteilungsdichten fur verschiedene r-Werte

Satz 2.7 (Momenterzeugende Funktion der NB-Verteilung) Die momenterzeugende

Funktion der NB-Verteilung lautet:

MY (t) = E(etY ) =(

p

1 − (1 − p) et

)r, t < | ln(1 − p)|.

Beweis: siehe [John 93]

Mit Hilfe der Formeln (2.1) und (2.2) ergeben sich fur Negativ-Binomial-verteilte Zufallsva-

riablen folgender Erwartungswert und folgende Varianz (siehe [Siko 02], Seite 34f):

E(Y ) = r1 − p

pund V ar(Y ) = r

1 − p

p2=E(Y )p

(2.34)

Dies bestatigt, dass sich die NB-Verteilung besser als die Poissonverteilung zur Modellierung

von Daten mit Uberdispersion eignet, da die Varianz immer großer als der Erwartungswert ist.

Satz 2.8 (Additivitat der NB-Verteilung) Fur die Zufallsvariable Z = X + Y mit

X ∼ NB(r, p) und Y ∼ NB(s, p) und X, Y unabhangig gilt:

Z = X + Y ∼ NB(r + s, p)


Beweis: Nach Voraussetzung sind die Zufallsvariablen X und Y unabhangig. Daher gilt fur die

momenterzeugende Funktion

MZ(t) = MX+Y (t) = Mx(t)MY (t) =(

p

1 − (1 − p)et

)r ( p

1 − (1 − p)et

)s=

(p

1 − (1 − p)et

)r+s.

Dies ist aber gerade die momenterzeugende Funktion einer NB(r+s,p)-verteilten Zufallsvariable.

Nun wollen wir die Definition der NB-Verteilung auf stetige Parameter erweitern. Da die Y-

Variable trotz der stetigen Parameter nur diskrete Werte annimmt, gehort auch die erweiterte

Form der NB-Verteilung noch zur Klasse der diskreten Verteilungen.

Definition 2.10 (Negativ-Binomial-Verteilung: erweiterte Form ) Die Zufallsvariable Y

heisst Negativ-Binomial-verteilt, falls ihre Wahrscheinlichkeitsfunktion folgende Form hat:

P (Y = y) =Γ(y + a−1)Γ(a−1)y!

(aµ

1 + aµ

)y ( 11 + aµ

)a−1

=Γ(y + a−1)Γ(a−1)y!

(µ

a−1 + µ

)y ( a−1

a−1 + µ

)a−1

, y = 0, 1, 2, . . .

Der Parameter a ≥ 0 wird oft als Dispersionsindex bezeichnet. Man schreibt: Y ∼ NB(µ, a)

Als nachstes sind wir analog zum Fall der vereinfachten NB-Verteilung an der Darstellung des

Erwartungswertes und der Varianz interessiert. Dazu benotigen wir die momenterzeugende Funk-

tion fur die erweiterte Form.

Satz 2.9 (Momenterzeugende Funktion der erweiterten NB-Verteilung) Die moment-

erzeugende Funktion fur die erweiterte NB-Verteilung ist gegeben durch:

MY (t) =

⎛⎝ a−1

a−1+µ

1 −(1 − a−1

a−1+µ

)et

⎞⎠a−1

Beweis: siehe [John 93]

Mit Hilfe der bereits bekannten Formeln (2.1) und (2.2) lassen sich aus der momenterzeugenden

Funktion der Erwartungswert und die Varianz berechnen. Man erhalt

E(Y ) =dM(0)dt

= µ und (2.35)

V ar(Y ) =d2M(0)dt2

− dM(0)dt

= µ+ aµ2. (2.36)


Nimmt man an, dass der stetige Parameter a nur diskrete Werte annimmt, lasst sich durch

leichtes Umformen die Aquivalenz der zwei Darstellungen der NB-Verteilung zeigen. Ersetzt man

den Paramter r durch a−1 und den Parameter p durch a−1

a−1+µ, muss noch gezeigt werden, dass

der Binomialkoeffizient in Definition 2.9 dem Γ-Bruch in der erweiterten Form (siehe Definition

2.10) entspricht:

Γ(y + a−1)Γ(a−1)y!

=(y + a−1 − 1)!(a−1 − 1)!y!

=(y + a−1 − 1

y

), fur a−1 ∈ N

Im Vergleich zur Poissonverteilung lasst also die Negativ-Binomial-Verteilung mehr Flexibilitat

bei der Modellierung der Varianz zu und eignet sich damit besser zur Modellierung von Daten mit

Uberdispersion. Die Poissonverteilung mit ihrer Eigenschaft der Gleichheit von Erwartungswert

und Varianz ist in der NB-Verteilung als Spezialfall enthalten. Wahlt man in der erweiterten

Form der NB-Verteilung den Dispersionsparameter a = 0 um auszudrucken, dass keine Uberdis-

persion vorliegt, erhalt man die Poissonverteilung.

Satz 2.10 (Die Poissonverteilung als Grenzfall der NB-Verteilung) Es sei Y eine

Negativ-Binomial-verteilte Zufallsvariable mit den Parametern µ und a = 0. Dies entspricht

einer Poisson-verteilten Zufallsvariable mit Parameter µ.

Beweis: siehe [Came 98]

P (Y = y) =Γ(y + a−1)Γ(a−1)y!

(aµ

1 + aµ

)y ( 11 + aµ

)a−1

α=a−1

=Γ(y + α)Γ(α)y!

(1αµ

1 + 1αµ

)y (1

1 + 1αµ

)α

=Γ(y + α)Γ(α)y!

(µ

α+ µ

)y ( α

α+ µ

)αMit

Γ(y + α)Γ(α)

=(y + α− 1)!

(α− 1)!= (y + α− 1)(y + α− 2) . . . (y + α− y) =

y−1∏j=0

(j + α)

folgt:

P (Y = y) =

⎛⎝y−1∏j=0

(j + α)

⎞⎠ 1y!

(1

α+ µ

)yµy

(α

α+ µ

)α

Mit (1

α+ µ

)y=

y−1∏j=0

(1

α+ µ

)


folgt:

P (Y = y) =

⎛⎝y−1∏j=0

j + α

α+ µ

⎞⎠ 1y!µy

(α

α+ µ

)α

=

⎛⎝y−1∏j=0

jα + 11 + µ

α

⎞⎠(

11 + µ

α

)α µyy!

Und wegen⎛⎝y−1∏j=0

jα + 11 + µ

α

⎞⎠ α→∞→ 1 und

(1

1 + µα

)α=

(α

α+ µ

)α=

(1 +

µ

α

)−α α→∞→ e−µ

ergibt sich

lima→0

P (Y = y) = limα→∞P (Y = y) = e−µ

µy

y!.

Dies ist gerade die Wahrscheinlichkeitsfunktion der Poissonverteilung mit Parameter µ.

2.6.2 Die NB-Verteilung mit bekanntem Parameter a als GLM

In Beispiel 2.2 wurde bereits gezeigt, dass die Poissonverteilung zur exponentiellen Familie

gehort. Dies hat den Vorteil, dass sie als zugrundeliegende Verteilung eines generalisierten li-

nearen Modells gewahlt werden kann. Diese Eigenschaft hat in etwas eingeschranktem Maße

auch die NB-Verteilung. Die Wahrscheinlichkeitsfunktion der erweiterten NB-Darstellung aus

Definition 2.10 lasst sich namlich schreiben als:

fNB(y; θ, φ) = exp

ln(

Γ(y + a−1)Γ(a−1)y!

)+ y ln

(aµ

1 + aµ

)− a−1 ln(1 + aµ)

(2.37)

Ist der Parameter a bekannt, gehort die NB-Verteilung zur linearen exponentiellen Familie (De-

finition 2.2) mit θ = ln(

aµ1+aµ

), a(φ) = 1, c(y, φ) = ln

(Γ(y+a−1)Γ(a−1)y!

)und b(θ) = a−1 ln(1 + aµ).

Da b(θ) in Abhangigkeit von θ dargestellt werden sollte, losen wir θ = ln(

aµ1+aµ

)nach µ auf und

setzen dies in b(θ) ein:

µ =eθ

a(1 + eθ)

b(θ) = a−1 ln(1 + aµ) = a−1 ln(

1 + aeθ

a(1 − eθ)

)

= a−1 ln(

1 +eθ

1 − eθ

)= a−1 ln

(1

1 − eθ

)


Damit lassen sich nun wieder der Erwartungswert und die Varianz bestatigen. Formel (2.3)

zufolge gilt:

E(Y ) = b′(θ) = a−1(1 − eθ)eθ

(1 − eθ)2

= a−1 eθ

1 − eθ= µ

Zur Uberprufung der Varianz benutzen wir wieder die Formel (2.4):

V ar(Y ) = b′′(θ)a(φ) = a−1 eθ(1 − eθ) + eθeθ

(1 − eθ)2

= a−1 eθ

1 − eθ+ a−1 (eθ)2

(1 − eθ)2= µ+ aµ2

Ist auch der Parameter a unbekannt, so ist die NB-Verteilung kein GLM mehr. Nach [Came 98]

(Seite 73) allerdings gehort sie zumindest noch zur Klasse der linearen exponentiellen Familie

mit dem Storparameter φ. Diese Familienklasse ist den GLMs sehr ahnlich (siehe [Came 98],

Seite 33).

2.6.3 Offset in der NB-Regression

Analog zur Poissonregression (siehe Abschnitt 2.2.3) werden wir nun den Offset in der NB-

Regression herleiten. Wieder benutzen wir die zusatzliche Variable ti, um die Zeiteinheit zu

standardisieren. Die zu untersuchende Variable Yi ist damit NB-verteilt mit den Parametern a

und tiµ∗i :

P (Yi = yi) =Γ(yi + a−1)Γ(a−1)yi!

(atiµ

∗i

1 + atiµ∗i

)yi(

11 + atiµ

∗i

)a−1

, i = 1, . . . , n; y = 0, 1, 2, . . .

Hierbei gilt µ∗i = exp(xtiβ).

Dies lasst sich fur bekanntes a in der Form einer exponentiellen Familie schreiben:

fNB(yi; θ, φ) = exp

ln(

Γ(yi + a−1)Γ(a−1)yi!

)+ yi ln

(atiµ

∗i

1 + atiµ∗i

)− a−1 ln(1 + atiµ

∗i )

Dabei gilt entsprechend Definition 2.2:

θi = ln(

atiµ∗i

1 + atiµ∗i

), a(φ) = 1, ci(yi, φ) = ln

(Γ(yi + a−1)Γ(a−1)yi!

)und b(θi) = a−1 ln(1 + atiµ

∗i )

Da b(θi) in Abhangigkeit von θi dargestellt werden sollte, losen wir θi = ln(

atiµ∗i1+atiµ∗i

)nach tiµ

∗i

auf und setzen dies in b(θi) ein:

tiµ∗i =

eθia(1 − eθi )


b(θi) = a−1 ln(1 + atiµ∗i ) = a−1 ln

(1 + a

eθia(1 − eθi )

)

= a−1 ln(

1 +eθi

1 − eθi

)= a−1 ln

(1

1 − eθi

)

Daraus folgt mit Formel (2.3) fur den Erwartungwert:

E(Yi) = µi = b′(θi) = a−1(1 − eθi )eθi

(1 − eθi )2

= a−1 eθi1 − eθi

= tiµ∗i = ti exp(xt

iβ) = exp(ln(ti) + xtiβ), i = 1, . . . , n

Die Varianz ergibt sich mit Formel (2.4) als:

V ar(Yi) = b′′(θi)a(φ) = a−1 eθi (1 − eθi ) + eθi e

θi

(1 − eθi )2

= a−1 eθi1 − eθi

+ a−1 (eθi )2

(1 − eθi )2= µi + aµ2

i , i = 1, . . . , n

Mit Hilfe des Erwartungswertes lasst sich nun wieder nach dem linearen Pradiktor auflosen, um

die Linkfunktion zu erhalten:

ηi = xtiβ = ln(µi) − ln(ti) = θi − ln(ti)︸︷︷︸

offset

= g(µi), i = 1, . . . , n

Sowohl in der Poissonregression als auch in der NB-Regression muss man also im Falle uneinheit-

licher Beobachtungszeitraume von der kanonischen Linkfunktion g(µi) = ln(µi) zur kanonischen

Linkfunktion g(µi) = ln(µi) − ln(ti) mit offset ln(ti) ubergehen.

2.6.4 Unbeobachtete Heterogenitat und gemischte Modelle

Scheitert die Poissonregression bei der Modellierung von Zahldaten, so kann ein Grund dafur

die Verletzung der vom Poissonmodell geforderten Aquidispersion sein. Meist außert sich diese

Verletzung in Form von Uberdispersion. Uberdispersion in der Poissonregression liegt vor, wenn

V ar(Yi) > E(Yi) gilt. Stellt man sich die Frage nach den Ursachen fur diese Uberdispersion,

so stoßt man auf viele Moglichkeiten. Am haufigsten sind fehlende relevante Regressoren ver-

antwortlich, oder es liegt - typisch fur die Poissonverteilung - unbeobachtete Heterogenitat vor.

Um diese in den Griff zu bekommen, wahlt man den Parameter der Poissonverteilung nicht fest,

sondern zufallig. Dies fuhrt uns zur Familie der gemischten Modelle.

Beispiel 2.8 (Das Poisson-Gamma-Mischmodell) Es gelte: Yi|Zi ∼ Poisson(Zi), Zi ≥ 0

sei eine Zufallsvariable, Zi und Yi seien unabhangig mit E(Zi) = µi, µi = ti exp(xtiβ).

Wahlt man den Parameter der Poissonverteilung nun nicht fest, sondern lasst ihn nach der

Gammaverteilung variieren, erhalt man ein gemischtes Modell - das Poisson-Gamma-Modell. Im


Folgenden wird gezeigt, dass das Ergebnis dieser Mischung eine NB-Verteilung ist. Wir nehmen

an, Zi sei Gamma-verteilt mit Erwartungswert µi und Index φµi. Das bedeutet:

Zi ∼ Γ(φµi, φ) mit E(Zi) =φµiφ

= µi

Fur die Dichte von Zi gilt damit:

f(zi) =1

Γ(φµi)(φzi)φµiexp−φzi 1

zi

Die Verteilung des Mischmodells erhalten wir durch Integration uber die gemischte Dichte bzgl.

zi:

P (Yi = yi) =∫ ∞

0P (Yi = yi|Zi = zi)fzi(zi)dzi

=∫ ∞

0

ezizyii

yi!1

Γ(φµi)(φzi)φµiexp−φzi 1

zidzi

=φφµi

Γ(φµi)yi!

∫ ∞

0e−zi(1+φ)z

(yi+φµi−1)i dzi

−zi(1+φ)=:si=φφµi

Γ(φµi)yi!

∫ ∞

0esi

s(yi+φµi−1)i

(1 + φ)yi + φµi − 11

1 + φdsi

=φφµi

Γ(φµi)yi!Γ(yi + φµi)

(1 + φ)yi+φµi

=Γ(yi + φµi)Γ(φµi)yi!

(1

1/φ+ 1

)φµi(

1/φ1 + 1/φ

)yi

Vergleichen wir dies mit der erweiterten Form der NB-Verteilung (siehe Definition 2.10), so

folgt: Yi ist NB-verteilt mit a−1i = φµi. Wir schreiben

Yi ∼ NB(µi,1φµi

).

Mit Gleichung (2.35) und Gleichung (2.36) erhalten wir :

E(Yi) = µi

V ar(Yi) = µi + aµ2i = µi +

1φµi

µ2i = µi +

1φµi = µi(1 + 1/φ)

Damit gilt V ar(Yi) > E(Yi) , falls 1/φ = 0, d.h. es liegt Uberdispersion vor. Falls φ = ∞, liegt

keine Uberdispersion vor.

2.6.5 Fisher-Informationsmatrix und Asymptotik der Schatzer

Am Ende des Abschnitts 2.3.2 wurde bereits erlautert, dass die Poisson-ML-Schatzer asymp-

totisch Normal-verteilt sind mit asymptotischer Kovarianzmatrix FI−1(β). Nun interessiert


uns die Asymptotik der Parameter β und a der NB-Regression. Die dafur notwendige Fisher-

Informationsmatrix wird im Folgenden in groben Schritten hergeleitet. Ausfuhrlich lasst sich

dies in [Siko 02] (Seite 47ff) finden.

Als erstes gilt es, die Likelihood-Funktion zu berechnen. Dabei wahlen wir die Log-Linkfunktion

fur das NB-Modell. Damit gilt:

Yi ∼ NB(µi, a), i = 1, . . . , n, unabhangig mitµi = exp(xtiβ)

Die Likelihood-Funktion lautet:

L(β, a) =n∏i=1

Γ(yi + a−1)Γ(a−1)yi!

(aµi

1 + aµi

)yi(

11 + aµi

)a−1

Mit Γ(x+ 1) = xΓ(x) fur x > 0 gilt:

Γ(yi + a−1)Γ(a−1)

= (yi + a−1 − 1) . . . (yi + a−1 − yi) =y∗i∏j=0

(j + a−1)

Dabei ist y∗i = yi − 1, und∏y∗ij=0 wird Null fur y∗i < 0. Die Likelihood-Funktion lasst sich damit

umschreiben zu:

L(β, a) =n∏i=1

⎡⎣⎛⎝ y∗i∏j=0

(j + a−1)

⎞⎠ 1yi!

(aµi

1 + aµi

)yi(

11 + aµi

)a−1⎤⎦

=n∏i=1

⎡⎣⎛⎝ 1ayi

y∗i∏j=0

(aj + 1)

⎞⎠ 1yi!ayi

(µi

1 + aµi

)yi(

11 + aµi

)a−1⎤⎦

=n∏i=1

⎡⎣⎛⎝ y∗i∏j=0

(aj + 1)

⎞⎠ 1yi!

(µi

1 + aµi

)yi(

11 + aµi

)a−1⎤⎦

Durch Logarithmieren erhalt man die Log-Likelihood-Funktion:

L(β, a) = lnL(β, a) =n∑i=1

⎡⎣ y∗i∑j=0

ln(aj + 1) − ln(yi!) + yi lnµi − (yi + a−1) ln(1 + aµi)

⎤⎦ (2.38)

Zur Bestimmung der Fisher-Informationsmatrix werden zunachst die ersten Ableitungen der

Log-Likelihood-Funktion nach β1 bis βp und nach a benotigt. Ersetzt man in der Log-Likelihood-

Funktion µi durch µi = exp(xtiβ), lasst sich wie gewohnt nach βr, r = 1, . . . , p ableiten. Man

erhalt

∂L(β, a)∂βr

=n∑i=1

xir(yi − µi)1 + aµi

fur r = 1, . . . , p

und

∂L(β, a)∂a

=n∑i=1

⎡⎣ y∗i∑j=0

(j

aj + 1

)+

1a2

ln(1 + aµi) − (yi + a−1)µi

1 + aµi

⎤⎦ .


Die Fisher-Informationsmatrix ist nach Satz 2.4 der Erwartungswert der negativen Hesse-Matrix.

Deshalb benotigen wir die zweiten Ableitungen der Log-Likelihood-Funktion. Diese lauten

∂2L(β, a)∂βr∂βs

=n∑i=1

−xirxisµi(1 + ayi)(1 + aµi)2

, (2.39)

∂2L(β, a)∂βr∂a

=n∑i=1

−xir(yi − µi)µi(1 + aµi)2

(2.40)

und

∂2L(β, a)∂a2

= −n∑i=1

⎡⎣ y∗i∑j=0

(j

aj + 1

)2

+2a3

ln(1 + aµi) − 2a2

µi1 + aµi

− (yi + a−1)µ2i

(1 + aµi)2

⎤⎦ . (2.41)

Mittels dieser Ableitungen wollen wir nun die Eintrage der Fisher-Informationsmatrix bestim-

men. Es handelt sich dabei um eine ((p+1)× (p+1))-Matrix, deren erster Block aus den Erwar-

tungswerten der negativen zweiten Ableitungen nach β besteht und mit FIr,s(β, a) abgekurzt

wird (r, s=1, . . . ,p). Aufgrund der Symmetrie der Fisher-Informationsmatrix sind sowohl in die

letzte Spalte als auch in die letzte Zeile die Erwartungswerte der negativen gemischten Ablei-

tungen einzutragen. Diese werden mit FIr,p+1(β, a) bezeichnet. Den Eintrag (p+1)× (p+1) der

Matrix belegt der Erwartungswert der negativen zweiten Ableitung nach a, kurz FIp+1,p+1(β, a).

Mit Gleichung (2.39) erhalten wir

FIr,s(β, a) =n∑i=1

µixirxis1 + aµi

fur r, s = 1, . . . , p. (2.42)

Desweiteren folgt aus Gleichung (2.40)

FIr,p+1(β, a) = 0 fur r = 1, . . . , p. (2.43)

Und schließlich liefert der Erwartungswert der negativen Ableitung (2.41):

FIp+1,p+1(β, a) = a−4

⎡⎣ n∑i=1

⎡⎣ ∞∑j=0

1(j + a−1)2

PNB(Yi ≥ j) − aµiµi + a−1

⎤⎦⎤⎦ (2.44)

Nun konnen wir die Fisher-Informationsmatrix in Blockform angeben.

FI(β, a) =

[FIr,s(β, a) 0

0t FIp+1,p+1(β, a)

]fur r, s = 1, . . . , p. (2.45)

In dieser Darstellung ist FIr,s eine (p×p)-Blockmatrix, und FIp+1,p+1 ist ein Skalar. 0 bezeichnet

einen Nullvektor der Lange p.


Wie im Poissonfall (siehe Ende des Abschnitts 2.3.2) sind auch die ML-Schatzer im NB-

Fall asymptotisch Normal-verteilt. In [Lawl 87] wird ohne konkreten Beweis erlautert, dass fur

a > 0 und unter schwachen Bedingungen fur die Kovariablenvektoren xi, i=1,. . . ,n, die die

Approximation von n−1FI(β, a) gegen einen positiv definiten Grenzwert fur n → ∞ sichern,

die ML-Schatzer als asymptotisch Normal-verteilt behandelt werden konnen. Man kann also fur

große n zeigen, dass

√n(β − β, a− a) a∼ N(0, nFI(β, a)−1) (2.46)

gilt. Dabei ist die Kovarianzmatrix gegeben durch

nFI(β, a)−1 = n

[FI−1

r,s (β, a) 0

0t FI−1p+1,p+1(β, a)

]. (2.47)

Die Eintrage FI−1r,s (β, a) und FI−1

p+1,p+1(β, a) ergeben sich aus den Gleichungen (2.42) und (2.44)

ausgewertet an (β, a) = (β, a). Referenzen fur Beweisansatze entnehme man [Siko 02] (Seite 55).

Man beachte noch, dass die Kovarianzmatrix eine Diagonal-Blockmatrix ist, da nach Gleichung

(2.43) die Parameter a und β asymptotisch unkorreliert sind.

Kapitel 3

Bayesianische Modelle und

Bayes-Faktoren zur Modellwahl

Bis jetzt galt die Annahme, dass bzgl. der Parameter keine Informationen außer den Daten

vorliegen. In vielen Situation sind allerdings weitere Informationen vorhanden.

Beispiel 3.1 (Qualitatssicherung) Eine Lieferung mit N Teilen soll auf ihre Qualitat uber-

pruft werden. Es ist aber zu teuer, alle Teile zu testen. Angenommen es fanden schon viele

Lieferungen mit N Teilen statt, und es liegen Kundeninformationen uber die Anzahl der defek-

ten Teile vor. Dann laßt sich die empirische Haufigkeitsverteilung aufstellen:

πi =Anzahl der Lieferungen mit i defekten Teilen

K

K sei dabei die Anzahl der vorherigen Lieferungen. θ sei der Anteil der defekten Teile. Damit

ist θ also eine ZV mit P (θ = i/N) = πi. Es gilt:

θ ∼ π = P (θ), π(i/N) = P (θ = i/N) = πi

X sei die Anzahl der defekten Teile in einer Stichprobe vom Umfang n. Damit folgt fur den

Likelihood:

P (X = k|θ = i/N) =

( ik

)(N−in−k

)(Nn

)Dabei beschreibt

( ik

)die Anzahl der Moglichkeiten, in der Stichprobe k defekte Teile von Nθ = i

defekten Teilen zu ziehen.(N−in−k

)gibt die Anzahl der Moglichkeiten an, in der Stichprobe (n− k)

nicht-defekte Teile von N −Nθ = N − i nicht-defekten Teilen zu ziehen. Geteilt wird durch(Nn

),

die Anzahl der Moglichkeiten, die Stichprobe zu ziehen. Daraus folgt:

P (X = k, θ = i/N) = P (X = k|θ = i/N)P (θ = i/N) = πi

( ik

)(N−in−k

)(Nn

)40

41

Bayesianischer Ansatz:

Bayesianer nehmen an, dass der Wert eines Parameters θ eine Realisierung einer ZV θ mit

gegebener Verteilung (a priori Verteilung) ist. Die a priori Verteilung erfasst die Annahmen des

Forschers uber den wahren Wert des Parameters, bevor die Daten erhoben werden. Man spricht

von ”subjektiver Inferenz“. Bezeichnen wir mit π(θ) die a priori Verteilungsdichte von θ, so gilt

θ ∼ π(θ) = P (θ)

Mit Hilfe der Likelihood-Funktion lasst sich - wie bereits in Beispiel 3.1 gezeigt - die gemeinsame

Verteilungsdichte berechnen:

P (X = x,θ) = π(θ)P (x|θ)

Nachdem die Daten erhoben sind, mochte man die a priori Verteilung fur θ dem neuen Wissen

anpassen. Dies fuhrt zur a posteriori Verteilung. Die bedingte Dichte bzw. Wahrscheinlichkeits-

funktion (a posteriori Verteilungsdichte) von θ gegeben die Daten lasst sich wie folgt berechnen:

π(θ|X = x) =P (θ,X = x)P (X = x)

=

⎧⎪⎪⎪⎨⎪⎪⎪⎩

π(θ)P (x|θ)∑t π(t)P (x|t) , falls θ diskret

π(θ)P (x|θ)∫ ∞−∞ π(t)P (x|t)dt , falls θ stetig

(3.1)

Dabei stellt

P (X = x) =

∑t π(t)P (x|t)∫∞

−∞ π(t)P (x|t)dt

die Marginalverteilung von X dar. Vernachlassigt man die normierende Marginalverteilung, so

gilt:

π(θ|X) ∝ π(θ)P (x|θ)

Das heisst, die a posteriori Verteilungsdichte ist proportional zum Produkt von a priori Vertei-

lungsdichte und Likelihood:

a posteriori Verteilungsdichte ∝ a priori Verteilungsdichte × Likelihood (3.2)

Beispiel 3.2 (Berechnung der a posteriori Verteilung) Es sei θ der Anteil der Amerika-

ner/innen, die die Todesstrafe befurworten. A priori - also vor der Datenerhebung - nehmen wir

an, dass θ Beta-verteilt ist mit α = 1 und β = 2/3. Wir gehen also von einem Erwartungswert

von E(θ) = αα+β = 0.6 aus. Fur die a priori Verteilung gilt:

π = P (θ) = Beta(α, β) =Γ(α+ β)Γ(α)Γ(β)

θα−1(1 − θ)β−1 =Γ(5/3)

Γ(1)Γ(2/3)(1 − θ)−1/3

42 Kapitel 3. Bayesianische Modelle und Bayes-Faktoren zur Modellwahl

Die Likelihood-Funktion lautet in diesem Beispiel

P (y|θ) = Bin(n, θ) = θy(1 − θ)n−y.

Mittels der in einer Umfrage unter 1000 Amerikanern bzw. Amerikanerinnen erhobenen Daten

berechnen wir die a posteriori Verteilung. Mit (3.2) folgt:

P (θ|y) ∝ P (y|θ)P (θ) ∝ θy(1 − θ)n−yθα−1(1 − θ)β−1 = θy+α−1(1 − θ)n−y+β−1

∝ Beta(y + α, n − y + β)

Wenn sich bei der Umfrage herausgestellt hat, dass 65% bzw. 650 der 1000 befragten Menschen

fur die Todesstrafe sind, ergibt sich fur y = 650, n = 1000, α = 1 und β = 2/3 als a posteriori

Verteilung

P (θ|y) ∝ Beta(651, 35023).

Mit α = 651 und β = 35023 erhalten wir einen a posteriori Erwartungswert von

E(θ) = αα+β = 0.65. Die a priori Verteilung wurde somit dem neuen Wissen angepasst.

3.1 Bayesianischer Ansatz und Definition der Bayes-Faktoren

Dieses Kapitel basiert großtenteils auf dem 1996 in der Zeitschrift Biometrika veroffentlichen

Technischen Report 255 von Adrian E. Raftery [Raft 94]. Nach der Definition der Bayes-Faktoren

werden Approximationen fur Bayes-Faktoren in Generalisierten Linearen Modellen (GLMs) vor-

gestellt, da eine direkte Berechunung meist nicht moglich ist. Die erste Approximation basiert auf

der Laplace-Methode fur Integrale, benotigt aber Werte, die oft nicht im Output von Standard-

Software zu finden sind. Deshalb stellt A.E. Raftery drei weitere Approximationen vor, die nur

den Output einer Standard-Software benotigen. Es wird eine empfohlene Menge geeigneter a

priori Verteilungen entwickelt, um einerseits die Situation, wenn wenig a priori Information vor-

liegt, zu reprasentieren und andererseits die Sensitivitat der Ergebnisse in Bezug auf die a priori

Verteilung zu beurteilen. Die von A.E. Raftery vorgeschlagenen Methoden konnen bei bekann-

tem und unbekanntem Dispersionsparameter benutzt werden, wenn Uberdispersion vorliegt und

um Linkfunktionen, Fehlerverteilungen und Varianzfunktionen zu vergleichen.

Die Modellbildung in GLMs umfasst die Auswahl der unabhangigen Variablen, der Linkfunk-

tion und der Varianzfunktion. Jede mogliche Auswahlkombination definiert ein eigenes Modell.

Damit besteht der Prozess der Modellwahl aus dem Vergleich vieler konkurrierender Model-

le. Gewohnlich benutzt man dazu Folgen von Signifikanztests, die oft auf der approximativen

asymptotischen Verteilung der Devianz basieren.

3.1. Bayesianischer Ansatz und Definition der Bayes-Faktoren 43

Dabei treten allerdings verschiedene Probleme auf. So wurde die Theorie der p-Werte zum

Vergleich zweier genesteter 1 Modelle entwickelt. Im Allgemeinen handelt es sich aber meist um

hunderte Modelle, die großtenteils nicht genestet sind. Grundsatzlich sagt A.E. Raftery, dass

die auf p-Werten basierenden Signifikanztests die falsche Frage stellen. Sie wollen wissen, ob

das Null-Modell ”wahr“ ist. Dies lasst sich aber von vorne herein mit ”nein“ beantworten. Eine

wissenschaftlich relevantere Frage ist: ”Welches Modell sagt die Daten besser voraus?“ Oder

anders gesagt: ”Unter welchem Modell sind die Daten wahrscheinlicher beobachtet worden?“

(vgl. [Raft 93])

Das vielleicht großte Problem ist allerdings folgendes: Jeder Ansatz, der ein einzelnes Modell

auswahlt und dann Folgerungen bedingt unter diesem Modell macht, ignoriert die Modellunsi-

cherheit. Und diese vernachlassigte Modellunsicherheit kann einen großen Teil der Gesamtunsi-

cherheit bzgl. der interessierenden Großen ausmachen. Die Unterschatzung der Modellunsicher-

heit fuhrt dann oft zu Entscheidungen, die zu riskant sind. All diese Komplikationen konnen im

Prinzip durch die Anwendung des Bayesianischen Ansatzes vermieden werden. Dabei berechnet

man die a posteriori Verteilung einer interessierenden Große als einen gewichteten Durchschnitt

ihrer a posteriori Verteilungen unter den individuellen Modellen gewichtet mit den a posteriori

Modell-Wahrscheinlichkeiten. Der Bayesianischen Ansatz tendiert im Vergleich zum sequentiel-

len p-Wert-Ansatz dazu, einfachere und sparsamere Modelle zu bevorzugen. Nach A.E. Raftery’s

Ansicht liefert uns die Modellwahl im Rahmen des Bayes-Ansatzes das Beste aus beiden Wel-

ten: Einfachere und interpretierbare Modelle, die durch die Daten gut unterstutzt werden und

dennoch die tatsachlich zu untersuchende Hypothese widerspiegeln und auf ausreichend durch-

dachten Theorien beruhen (vgl. [Raft 93]).

In den folgenden Abschnitten werden zunachst die Grundideen der Bayes-Faktoren, der a

posteriori Wahrscheinlichkeiten und der Beachtung der Modellunsicherheit vorgestellt. Dann be-

nutzen wir die Laplace-Methode fur Integrale, um eine allgemeine Approximation fur die Bayes-

Faktoren zu erhalten. Daraufhin werden drei weitere sehr genaue Approximationen beschrieben,

die nur die MLEs der Parameter, die Devianz und die Informationsmatrix benutzen. Dement-

sprechend konnen diese Approximationen direkt aus dem Output einer Standard-Software zum

Schatzen generalisierter linearer Modelle berechnet werden.

Wir beginnen mit einem Datensatz D, von dem wir annehmen, dass er unter einem der

zwei Modelle M1und M0 entsprechend einer Wahrscheinlichkeitsdichte P (D|M1) oder P (D|M0)

entstanden ist. Sind die a priori Wahrscheinlichkeiten P (M1) und P (M0) = 1 − P (M1) gege-

ben, produziert der Datensatz die a posteriori Wahrscheinlichkeiten P (M1|D) und P (M0|D) =

1 − P (M1|D). Jede a priori Meinung wird durch Betrachten der Daten zu einer a posterio-

ri Meinung transformiert. Diese Transformation reprasentiert die Evidenz fur Modell 1 gegen

1Ein Modell Mi ist in Modell Ms genestet, wenn alle Kovariablen von Modell Mi im Modell Ms, s ≥ i + 1,

enthalten sind. M1 ⊂ M2 ⊂ . . . ⊂ Mr. M1 = Null-Modell, Mr =”saturated“ Modell


Modell 0, die die Daten liefern. Wenn wir zur sogenannten odds-Skala ubergehen (odds =

Wahrscheinlichkeit/(1 − Wahrscheinlichkeit)), bekommt die Transformation eine einfache

Form.

Die Transformation des a priori odds zum a posteriori odds bezeichnet man als Bayes-Faktor.

Um diesen herzuleiten, benotigen wir den Satz von Bayes:

Satz 3.1 (Satz von Bayes) Sei A1, . . . , Ak eine disjunkte Zerlegung von Ω, wobei fur minde-

stens ein i, i = 1, . . . , k, P (Ai) > 0 und P (B|Ai) > 0 erfullt ist. Dann gilt:

P (Aj |B) =P (B|Aj)P (Aj)∑ki=1 P (B|Ai)P (Ai)

=P (B|Aj)P (Aj)

P (B), j = 1, . . . , k

Beweis: Der Satz von Bayes folgt aus dem Produktsatz (Ps) und dem Satz von der totalen

Wahrscheinlichkeit (SvtW):

P (Aj |B) =P (B ∩Aj)P (B)

(Ps)=

P (B|Aj)P (Aj)P (B)

(SvtW )=

P (B|Aj)P (Aj)∑ki=1 P (B|Ai)P (Ai)

Mit Hilfe des Satzes von Bayes erhalten wir:

P (Mk|D) =P (D|Mk)P (Mk)

P (D|M1)P (M1) + P (D|M0)P (M0)(k = 0, 1) (3.3)

Daraus folgt

P (M1|D)P (M0|D)

=P (D|M1)P (D|M0)

P (M1)P (M0)

, (3.4)

d.h. die Transformation des a priori odds P (M1)P (M0) zum a posteriori odds P (M1|D)

P (M0|D) ist eine einfache

Multiplikation mit

B10 :=P (D|M1)P (D|M0)

(3.4)=

P (M1|D)P (M0|D)

P (M1)P (M0)

=P (M1|D)

1−P (M1|D)

P (M1)1−P (M1)

. (3.5)

B10 wird als Bayes-Faktor bezeichnet. Er stellt das Verhaltnis des a posteriori odds von Modell

M1 zu seinem a priori odds dar. Oft wahlt man fur den a priori odds den Wert 1, um neutrale a

priori Information zu reprasentieren, die kein Modell bevorzugt. Der Bayes-Faktor ist eine von

den Daten gelieferte Zusammenfassung der Evidenz fur Modell M1 gegen Modell M0.

Gleichung (3.4) in Worten:

a posteriori odds = Bayes-Faktor × a priori odds

3.1. Bayesianischer Ansatz und Definition der Bayes-Faktoren 45

In Gleichung (3.4) gilt fur den marginalen Likelihood :

P (D|Mk) =∫P (D|θk,Mk)P (θk|Mk)dθk (k = 0, 1), (3.6)

wobei θk der Parametervektor von Mk ist. P (θk|Mk) ist seine a priori Dichte, und P (D|θk,Mk)

der Likelihood. Der marginale Likelihood P (D|Mk) ist also die Vorhersagewahrscheinlichkeit

der Daten gegeben das Modell Mk und misst, wie gut das Modell die Daten vorhersagt. Nach

A.E. Raftery ist dies das richtige Kriterium zur Modellbewertung. Der Bayes-Faktor gibt an,

wie gut das Modell M1 die Daten im Vergleich zum Modell M2 prophezeit. Es sei noch erwahnt,

dass dieser Vergleich nicht von der Annahme abhangt, dass ein Modell ”wahr“ ist, wie es bei

Modellvergleichen der Fall ist, die auf p-Werten beruhen.

Oft kann es nutzlich sein, den zweifachen Logarithmus des Bayes-Faktors zu betrachten. Dann

bewegt man sich auf der Skala der gewohnlichen Devianz und der Likelihood-Ratio-Statistiken.

Den zweifachen Logarithmus des Bayes-Faktors (2 logB10) wollen wir in Zukunft als skalierten

Bayes-Faktor bezeichnen. Wir benutzen die gerundete Skala in Tabelle 3.1 zur Interpretation

des Bayes-Faktors bzw. des skalierten Bayes-Faktors. Diese Skala basiert auf der von H. Jeffreys

(1961) [Jeff 61], ist aber etwas grober und nicht so konservativ. Auf dieser Skala uberschreitet

B10 2 logB10 Evidenz fur Modell M1

< 1 < 0 Negativ (unterstutzt Modell M0)

1 bis 3 0 bis 2.2 Nicht mehr wert als bloße Erwahnung

3 bis 12 2.2 bis 5 Positiv

12 bis 150 5 bis 10 Stark

> 150 > 10 Entscheidend, maßgeblich

Tabelle 3.1: Tabelle fur die Interpretation der Bayes-Faktoren

der Bayes-Faktor B10 fur eine Reihe angemessener a priori Verteilungen selten mehr als eine

Intervallgrenze. Auf Jeffreys Originalskala hingegen nimmt der Bayes-Faktor oft gleich zwei

Hurden. Das erschwert die Interpretation.

Werden mehr als zwei Modelle betrachtet, erhalt man mit Hilfe der Bayes-Faktoren fur

alle Modelle folgendermaßen die a posteriori Wahrscheinlichkeiten. Wir nehmen an, es werden

(K + 1) Modelle M0,M1, . . . ,MK betrachtet. Vergleicht man die Modelle M1, . . . ,MK jeweils

mit M0, ergeben sich die Bayes-Faktoren B10, . . . , BK0. Entsprechend Gleichung (3.3) laßt sich

fur k = 0, . . . ,K die a posteriori Wahrscheinlichkeit von Mk wie folgt berechnen:

P (Mk|D) =P (D|Mk)P (Mk)∑Kr=0 P (D|Mr)P (Mr)

=P (D|Mk)P (D|M0)

P (Mk)P (M0)∑K

r=0P (D|Mr)P (D|M0)

P (Mr)P (M0)

(3.5)=

αkBk0∑Kr=0 αrBr0

, (3.7)

wobei αk := P (Mk)/P (M0) der a priori odds fur Mk gegen M0 ist (k = 0, . . . ,K). Ausserdem

gelte B00 = α0 = 1.


Nun wollen wir die a posteriori Verteilung einer interessierenden Große ∆, meist einer zukunf-

tigen Beobachtung, berechnen. In der Medizin interessiert z.B. die Wahrscheinlichkeit, dass

jemand, der eine bestimmte Krankheit jetzt nicht hat, diese spater bekommen wird. Nach Glei-

chung (3.6) lautet die a posteriori Verteilung gegeben die Daten fur jedes Modell Mk

P (∆|D,Mk) =∫P (∆|D,θk,Mk)P (θk|D,Mk) dθk.

Dabei gilt nach Gleichung (3.1), wenn θk der Parameter von Modell Mk ist:

P (θk|D,Mk) =P (D|θk,Mk)P (θk|Mk)∫P (D|θk,Mk)P (θk|Mk) dθk

(3.8)

Die a posteriori Verteilung unserer interessierenden Große ∆ unter den moglichen Modellen

M1, . . . ,Mk gegeben die Daten berechnen wir nun als gewichteten Durchschnitt der a posteriori

Verteilungen von ∆ unter jedem Modell P (∆|D,Mk) gewichtet durch ihre a posteriori Modell-

wahrscheinlichkeiten P (Mk|D) (Bayesian-mixing):

P (∆|D) :=K∑k=0

P (∆|D,Mk)P (Mk|D) (3.9)

Hierbei lasst sich P (Mk|D) mit Hilfe der Formel (3.7) berechnen.

Der gemeinsame a posteriori Erwartungswert und die Standardabweichung (bzw. der Punktschatzer

und der Standardfehler) konnen mit folgenden Formeln berechnet werden:

E[∆|D] =K∑k=0

E[∆|D,Mk]P (Mk|D) =K∑k=0

∆k P (Mk|D) (3.10)

V ar[∆|D] = E[(∆|D)2] − (E[∆|D])2

(3.9)=

K∑k=0

E[(∆|D,Mk)2]P (Mk|D) − (E[∆|D])2

=K∑k=0

[E[(∆|D,Mk)2] + (E[∆|D,Mk])2

−(E[∆|D,Mk])2]P (Mk|D) − (E[∆|D])2

=K∑k=0

(V ar[∆|D,Mk] + ∆2

k

)P (Mk|D) − (E[∆|D])2 (3.11)

Dabei gilt ∆k := E[∆|D,Mk]. Der Erwartungswert und die Varianz lassen sich aus dem Output

einer Standard-Maximum-Likelihood-Analyse berechnen, indem man ∆k durch den MLE und

V ar[∆|D,Mk] durch das Quadrat des entsprechenden Standardfehlers ersetzt (siehe [Raft 93]).

3.2. Approximation der Bayes-Faktoren 47

3.2 Approximation der Bayes-Faktoren

Nur in einigen elementaren Fallen kann das zur Berechnung der Bayes-Faktoren notige Integral

(3.6) analytisch bestimmt werden. Viel ofter ist es unlosbar. Deshalb muss auf numerische Me-

thoden zuruckgegriffen werden. Um die damit verbundenen Probleme zu umgehen, werden wir

die Integrale nur approximativ berechnen. Zu diesem Zweck werden im Folgenden vier unter-

schiedlich genaue Approximationen vorgestellt.

3.2.1 Die Laplace-Methode fur Integrale

Grundsatzlich basiert die Laplace-Methode fur Integrale (siehe [Brui 70], Kapitel 4.4) auf einer

Taylor-Reihen-Entwicklung der reellwertigen Funktion h(θ) des p-dimensionalen Vektors θ und

fuhrt zu folgender Approximation:∫eh(θ) dθ ≈ (2π)p/2|A| 12 exph(θ) (3.12)

Dabei ist θ der Wert von θ, an dem die Funktion h ein Maximum hat; A ist minus die Inverse

der Hesse-Matrix von h an der Stelle θ. Die Schreibweise |A| steht fur det(A) und bezeichnet die

Determinante von A. Diese Approximation wollen wir zur Berechnung des Integrals in Gleichung

(3.6) benutzen. Dabei gilt

h(θ) := logP (D|θ,M)P (θ|M). (3.13)

Herleitung der Laplace-Approximation in groben Schritten:

Vorerst beschranken wir uns auf den eindimensionalen Fall (p=1) und approximieren folgendes

Integral:

I :=∫ 1

0exph(θ)dθ (3.14)

Dabei ist h(θ) = logP (D|θ,M)P (θ|M) stetig, unimodal, zweimal differenzierbar und hat in

(0, 1) ein Maximum bei θ. Eine Taylorreihenentwicklung fur h(θ) um θ ergibt:

I ≈∫ 1

0exph(θ) + h′(θ)(θ − θ) +

12h′′(θ)(θ − θ)2dθ (3.15)

Nach Definition von θ gilt:

h′(θ) =∂

∂θ[logP (D|θ,M)P (θ|M)]θ=θ = 0

Die Approximation (3.15) vereinfacht sich zu:

I ≈ exph(θ)∫

exp−12(−h′′)(θ)(θ − θ)2dθ (3.16)


Die Grenzen des Integrals wurden heuristisch auf ∞ erweitert. Damit mussen wir also uber eine

nicht-normalisierte Dichte der Normalverteilung integrieren. Wie allgemein bekannt, gilt fur die

Dichte einer Normalverteilung N(µ, σ2) mit Erwartungswert µ und Varianz σ2∫1√

2πσ2exp− 1

2σ2(x− µ)2dx = 1. (3.17)

Daraus folgt ∫exp− 1

2σ2(x− µ)2dx =

√2πσ2. (3.18)

Entsprechend ergibt sich fur das Integral in Approximation (3.16) mit µ = θ und σ2 = 1−h′′(θ) :∫

exp−12(−h′′)(θ)(θ − θ)2dθ =

√2π

1−h′′(θ) (3.19)

Die eindimensionale Approximation des Integrals I lautet also:

I ≈ (2π)12

[1

−h′′(θ)

] 12

exph(θ) (3.20)

Somit ergibt sich fur die Approximation im p-dimensionalen Fall

I ≈ (2π)p2 |A| 12 exph(θ). (3.21)

Dabei ist A gleich minus die Inverse der Hessematrix, und |A| bezeichnet die Determinante der

Matrix A.

Diese Approximation wenden wir nun auf das Integral in Gleichung (3.6) an. Es ergibt sich:

P (D|Mk)(3.6)=

∫P (D|θk,Mk)P (θk|Mk)dθk

(3.12)≈ (2π)pk/2|Ψk|12P (D|θk,Mk)P (θk|Mk) (3.22)

Dabei ist pk die Dimension von θk, θk der a posteriori Modus von θk und Ψk minus die Inverse

der Hesse-Matrix von h(θk) = logP (D | θk,Mk)P (θk |Mk) an der Stelle θk = θk.

In regularen statistischen Modellen ist der relative Fehler in der Approximation (3.22) und

damit auch in der resultierenden Approximation fur den Bayes-Faktor B10 von der Ordnung

O(n−1) (vgl. [Tier 86]). Ein statistisches Modell bezeichnet man als regular, wenn entweder

alle Verteilungen stetig oder alle Verteilungen diskret sind. In einem statistischen Modell ist im

Vergleich zum Wahrscheinlichkeitsraum das Wahrscheinlichkeitsmaß nur bis auf den Parameter

θ bekannt. Unter Benutzung von (3.22) kann man den marginalen Likelihood P (D|Mk) (siehe

Gleichung (3.6)) in jedem regularen statistischen Modell approximieren. Jedoch liefert Standard-

Software gewohnlich nicht den a posteriori Modus θk und das Negative der inversen Hesse-Matrix


Ψk, sondern berechnet oft den Maximum-Likelihood-Schatzer (MLE) θk, die Devianz oder die

Likelihood-Ratio-Statistik und die beobachtete oder erwartete Fisher-Informationsmatrix Fk

oder ihre Inverse Vk.

Im Folgenden werden deshalb drei unterschiedlich genaue Approximationen hergeleitet, die

auf der Laplace-Approximation (3.22) basieren, aber nur die im Output von Standard-Software

verfugbaren Großen benotigen.

3.2.2 Drei weitere Approximationen

Angenommen fur die a priori Verteilung von θk gilt E[θk |Mk] =: ωk und V ar[θk |Mk] =: Wk.

Approximiert man den a posteriori Modus θk ausgehend vom MLE θk durch einen einzigen

Schritt des Newton-Raphson-Algorithmus und setzt das Ergebnis in Gleichung (3.22) ein, erhalt

man als erste Approximation fur den skalierten Bayes-Faktor:

2 logB10 ≈ χ2 + (E1 − E0) (3.23)

In dieser Gleichung gilt χ2 = 2L1(θ1) − L0(θ0) , wobei Lk(θk) = log(P (D|θk,Mk)) der Log-

Likelihood ist; χ2 ist die Standard-Likelihood-Ratio-Teststatistik, wenn M0 in M1 genestet ist.

Ferner gilt:

Ek = 2λk(θk) + λ′k(θk)T (Fk +Gk)−12I − Fk(Fk +Gk)−1λ′k(θk)

− log |Fk +Gk| + pk log(2π), (3.24)

wobei Gk := (V ar[θk|Mk])−1 = W−1k , λk(θk) := logP (θk|Mk) die Log-Prior-Dichte und λ′k(θk)

der pk-dimensionsionale Vektor der Ableitungen von λk(θk) bezuglich der Elemente von θk ist.

Mit Fk bezeichnen wir die beobachtete oder erwartete Fisher-Informationsmatrix (k = 0, 1).

Rechtfertigung der Approximation (3.23):

Mittels der Laplace-Approximation (3.22) folgt:

2 logB10 ≈ 2 log(P (D|M1)P (D|M0)

)

≈ 2 log

((2π)p1/2|Ψ1|1/2P (D|θ1,M1)P (θ1|M1)(2π)p0/2|Ψ0|1/2P (D|θ0,M0)P (θ0|M0)

)

Mit Lk(θk) = log(P (D|θk,Mk)) und λk(θk) := log(P (θk|Mk) erhalt man:

2 logB10 = 2(p1

2log(2π) +

12

log |Ψ1| + L1(θ1) + λ1(θ1)

−p0

2log(2π) − 1

2log |Ψ0| − L0(θ0) − λ0(θ0)) +O(n−1)

= 2L1(θ1) − L0(θ0)

+ 2

λ1(θ1) − λ0(θ0)

+ log |Ψ1| − log |Ψ0| + (p1 − p0) log(2π) +O(n−1) (3.25)


Die Indizes 0 und 1 werden ab jetzt weggelassen, da sie nicht zum besseren Verstandnis bei-

tragen. Ein Schritt der Newton-Raphson-Methode (siehe Kapitel 2, Seite 12) fuhrt zu folgender

Approximation des a posteriori Modus θ:

θ ≈ θ − h′′(θ)−1h′(θ) (3.26)

Mit h(θ) = −L(θ) − λ(θ) = − log(P (D|θ,M)) − log(P (θ|M)) folgt [Berg 85] (Seite 224)

h′′(θ) ≈ −(F +G). (3.27)

Beweis:

Wir zeigen zuerst, dass die a posteriori Dichte π(θ|D) ∝ P (D|θ,M)π(θ|M) = exp(h(θ)) durch

Np(θ, F I−1obs(θ)) approximiert werden kann mit FIobs die beobachtete Fisher-Informationsmatrix

mit FIobs(θ)i,j = −[ ∂2

∂θi∂θjlog(P (D|θ,M))]

θ=θ. Daraus leiten wir uber weitere Approximationen

die obige Behauptung her. Zur Vereinfachung betrachten wir nur den Fall p = 1. Eine Taylor-

Entwicklung des Log-Likelihood L(θ) = logP (D|θ,M) um θ ergibt fur θ nahe bei θ:

π(θ|D) =explog(P (D|θ,M)π(θ)∫explog(P (D|θ,M)π(θ)dθ

≈ explog(P (D|θ,M) − 12(θ − θ)2FIobs(θ)π(θ)∫

explog(P (D|θ,M) − 12(θ − θ)2FIobs(θ)π(θ)dθ

=explog(P (D|θ,M)π(θ) exp−1

2(θ − θ)2FIobs(θ)explog(P (D|θ,M)π(θ)

∫exp−1

2(θ − θ)2FIobs(θ)dθ

=exp−1

2(θ − θ)2FIobs(θ)[2πFIobs(θ)−1

] 12

Dabei beachte man, dass die a priori Dichte π(θ) fur θ nahe θ approximativ konstant ist, und die

erste Ableitung des Log-Likelihood L(θ) = log(P (D|θ,M) an der Stelle θ Null ist. Die a posteriori

Dichte kann also durch N(µ, σ2) approximiert werden mit µ = θ und σ2 =[FIobs(θ)

]−1.

Es gilt ausserdem θ ≈ θπ

mit θπ

der verallgemeinerte MLE fur θ, d.h. θπ

maximiert

L∗(θ) = P (D|θ,M)π(θ). Desweiteren gilt mit FI die erwartete Fisher-Informationsmatrix

[FI(θ)]−1 ≈ [FIobs(θ)]−1 ≈ [FI(θ)π]−1 =[−[

∂2

∂θ∂θlog(P (D|θ,M)π(θ))

]]−1

.

Die a posteriori Dichte π(θ|D) kann damit durch Np(θπ,[FI(θ)π

]−1) approximiert werden.

Es folgt

Ψ =[−h′′(θ)

]−1

θ=θ≈

[−[∂2

∂θ∂θlog(P (D|θ,M)) +

∂2

∂θ∂θπ(θ))

]]−1

θ=θ= (F +G)−1


mit G = W−1, W = V ar[θ|M ] die a priori Kovarianzmatrix von θ im Modell M und F gleich die

beobachtete oder erwartete Fisher-Informationsmatrix. Zudem gilt h′(θ) = L′(θ)+λ′(θ) = λ′(θ),

da die Ableitung der Log-Likelihood-Funktion an der Stelle θ nach der Definition des MLE Null

ist. Die Approximation (3.26) wird damit zu:

θ ≈ θ + (F +G)−1λ′(θ) ⇒ θ − θ ≈ (F +G)−1λ′(θ) (3.28)

Mittels multivariater Taylor-Reihen-Entwicklung fur den Log-Likelihood L(θ) folgt:

L(θ) ≈ L(θ) +12(θ − θ)tL′′(θ)(θ − θ)

(3.28)≈ L(θ) − 12λ′(θ)t(F +G)−1F (F +G)−1λ′(θ) (3.29)

Dabei haben wir benutzt, dass L′′(θ) ≈ −F gilt. Analog folgt:

λ(θ) ≈ λ(θ) + λ′(θ)t(θ − θ) (3.30)(3.28)≈ λ(θ) + λ′(θ)t(F +G)−1λ′(θ)

Desweiteren gilt wie soeben bewiesen Ψ ≈ (F +G)−1 (siehe [Berg 85], Seite 224). Dies und die

Approximationen (3.29) und (3.30) setzen wir nun in die Gleichung (3.25) ein:

2 logB10

(3.29)(3.30)≈ 2L1(θ1) − 1

2λ′1(θ1)t(F1 +G1)−1F1(F1 +G1)−1λ′1(θ1)

−L0(θ0) +12λ′0(θ0)t(F0 +G0)−1F0(F0 +G0)−1λ′0(θ0)

+2λ1(θ1) + λ′1(θ1)t(F1 +G1)−1λ′1(θ1)

−λ0(θ0) − λ′0(θ0)t(F0 +G0)−1λ′0(θ0)

+ log | (F1 +G1)−1 | − log | (F0 +G0)−1 | +(p1 − p0) log(2π)

= 2L1(θ1) − L0(θ0)

+2λ1(θ1) − λ′1(θ1)t(F1 +G1)−1F1(F1 +G1)−1λ′1(θ1)

+2λ′1(θ1)t(F1 +G1)−1λ′1(θ1) + log |(F1 +G1)−1| + p1 log(2π)

−2λ0(θ0) + λ′0(θ0)t(F0 +G0)−1F0(F0 +G0)−1λ′0(θ0)

−2λ′0(θ0)t(F0 +G0)−1λ′0(θ0) − log |(F0 +G0)−1| − p0 log(2π)(∗)= χ2

+2λ1(θ1) + λ′1(θ1)t(F1 +G1)−12I − F1(F1 +G1)−1λ′1(θ1)

− log |F1 +G1| + p1 log(2π)

−2λ0(θ0) − λ′0(θ0)t(F0 +G0)−12I − F0(F0 +G0)−1λ′0(θ0)

+ log |F0 +G0| − p0 log(2π)

= χ2 + (E1 −E0) (3.31)


mit χ2 = 2L1(θ1) − L0(θ0), wobei Lk(θk) = log(P (D|θk,Mk)) der Log-Likelihood ist; χ2

ist die Standard-Likelihood-Ratio-Teststatistik, wenn M0 in M1 genestet ist. Ek entnehme man

Gleichung (3.24). Die Gleichheit in (*) gilt, da log |(Fk+Gk)−1| = − log |Fk+Gk| wegen detA−1 =

(detA)−1 fur A ∈ GL(n;K). Dabei stellt K einen Korper dar, z.B. R, und die Menge GL(n;K) =

A ∈ M(n × n);K) : A invertierbar ist mit der Multiplikation von Matrizen eine Gruppe mit

neutralem Element En. Sie heisst allgemeine lineare Gruppe (general linear). (Fk +Gk) ist ein

Element der allgemeinen linearen Gruppe der Matrizen.

Diese Approximation ist naher an der Laplace-Approximation (3.22), wenn Fk die beobachte-

te Fisher- Informationsmatrix ist, als wenn man die erwartete Fisher-Information benutzt. Wenn

Fk die beobachtete Fisher-Information ist, belauft sich der relative Fehler auf O(n−1), wahrend

der Einsatz der erwarteten Fisher-Information den relativen Fehler auf O(n−1/2) erhoht (ver-

gleiche [Kass 92b]).

Nun nehmen wir an, dass die a priori Verteilung von θ die Normalverteilung ist, d.h. es gilt

θ ∼ N(ω,W ) mit W = V ar[θ|M ] und ω = E[θ|M ].

Die n-dimensionale Normalverteilung hat folgende Form:

f(θ) = P (θ|M) =1

(2π)p2

√|W | exp−12(θ − ω)tW−1(θ − ω)

Mit G = W−1 und λ(θ) = log(P (θ|M) folgt:

λ(θ) = log

1

(2π)p2

√|W | exp−12(θ − ω)tW−1(θ − ω)

=12

log |G| − 12(θ − ω)tG(θ − ω) − p

2log(2π) (3.32)

λ′(θ) = −G(θ − ω) (3.33)

Setzt man diese Werte in Gleichung (3.24) ein, ergibt sich:

Ek = 2(

12

log |Gk| − 12(θk − ωk)tGk(θk − ωk) − pk

2log(2π)

)+(θk − ωk)Gk(Fk +Gk)−1

2I − Fk(Fk +Gk)−1

Gk(θk − ωk)

− log |Fk +Gk| + pk log(2π)

= log |Gk| − (θk − ωk)tGk(θk − ωk)

+(θk − ωk)Gk(Fk +Gk)−12I − Fk(Fk +Gk)−1

Gk(θk − ωk)

− log |Fk +Gk| (3.34)


Mit Hk := (Fk +Gk)−1 erhalten wir:

Ek = log |Gk| − (θk − ωk)tGk(θk − ωk)

+(θk − ωk)GkHk 2I − FkHkGk(θk − ωk) − log |Fk +Gk|= log |Gk| − (θk − ωk)tGk I −Hk(2I − FkHk)Gk (θk − ωk)

− log |Fk +Gk| (3.35)

Zur ubersichtlicheren Darstellung setzen wir noch Ck := Gk I −Hk(2I − FkHk)Gk:

Ek = log |Gk| − (θk − ωk)tCk(θk − ωk) − log |Fk +Gk| (3.36)

Die Gleichung (3.36) erlaubt, Ek einfacher zu berechnen, wenn man a priori annimmt, dass θ

Normal-verteilt ist.

Wir wenden uns nun einer zweiten Approximation zu. Diese ist einfacher aber weniger genau.

Dabei nimmt man Gleichheit in den approximativen Relationen θk ≈ θk und Ψk ≈ F−1k an, so

dass gilt

2 logB10 ≈ χ2 + (E∗1 − E∗

0) (3.37)

mit

E∗k = − log |Fk| + 2λk(θk) + pk log(2π). (3.38)

Dies folgt analog zur Rechtfertigung der ersten Approximation aus

θ = θ

Ψk = F−1k

L(θ) = L(θ)

λ(θ) = λ(θ),

denn setzt man diese Werte wieder in (3.25) ein, folgt:

2logB10 ≈ 2L1(θ1) −L0(θ0)

+ 2

λ1(θ1) − λ0(θ0)

(3.39)

− log |F1| + log |F0| + p1 log(2π) − p0 log(2π) = χ2 + (E∗1 − E∗

0)

Wenn die a priori Verteilung von θ die Normalverteilung ist, vereinfacht sich Gleichung (3.38)

zu:

E∗k = − log |Fk| − (θk − ωk)tGk(θk − ωk) + log |Gk|, (3.40)

da λ(θ) = 12 log |G| − 1

2(θ − ω)tG(θ − ω) − p2 log(2π) (siehe Gleichung (3.32)).

Obwohl diese Approximation nicht so gut wie die vorhergehende ist, wird sie hier aufgefuhrt,


weil sie relativ gut arbeitet und fur analytische Arbeit den Vorteil hat, dass die Beitrage von a

priori Verteilung und Likelihood in separaten Termen stehen.

In Gleichung (3.38) ist jedes Element der Fisher-Informationsmatrix Fk von der Ordnung O(N),

wobei N die Gesamtstichprobengroße ist. Im Poisson-Fall ist N typischerweise die Summe der

Beobachtungen. Somit gilt:

|Fk| = O(Npk)

Weiter gilt:

log |Fk| = log(O(Npk)) = log(NpkO(1)) = log(Npk) + log(O(1)) = pk logN +O(1)

Dies fuhrt zur dritten Approximation, die die einfachste aber auch die ungenaueste ist:

2 logB10 ∼ χ2 − (p1 − p0) logN = χ2 + (E+1 − E+

0 ) (3.41)

mit

E+k = −pk logN (3.42)

In dieser Approximation belauft sich der relative Fehler auf O(1). Jedoch zeigen praktische

Experimente, dass diese Approximation angesichts ihres asymptotischen Fehlers dennoch uber-

raschend genau ist. In [Kass 92a] wird auf einen theoretischen Grund dafur hingewiesen. Es lasst

sich zeigen, dass beim Vergleich genesteter Modelle unter bestimmten Bedingungen der Fehler

in Gleichung (3.41) nur von der Ordnung O(n−12 ) ist, wenn die Menge an Information in der a

priori Verteilung gleich der in einer Beobachtung ist, und die Alternativhypothese ”nahe“ bei

Null ist.

Die Ergebnisse dieses Kapitels werden nun in Form zweier Tabellen zusammengefasst. In der

ersten Tabelle finden sich alle verwendeten Terme mit ihrer Beschreibung, in der zweiten Tabelle

erhalt man einen Uberblick uber die drei hergeleiteten Approximationen.


Term (k=0,1) Erklarung

h(θk) := logP (D|θk,Mk)P (θk|Mk) abkurzende Schreibweise

minus die Inverse der Hesse-Matrixψk := [−h′′(θk)]−1

von h(θk)

θk a posteriori Modus von θk

θk ML-Schatzer von θk

erwartete (beobachtete)Fk (siehe Definition 2.6)

Fisher-Informationsmatrix

Inverse derVk := F−1

k Fisher-Informationsmatrix

a priori Erwartungswertωk := E[θk|Mk] von θk im Modell Mk

a priori KovarianzmatrixWk := V ar[θk|Mk] von θk im Modell Mk

Inverse der a priori KovarianzmatrixGk := W−1

k von θk im Modell Mk

Lk(θk) Log-Likelihood fur Modell Mk

Likelihood-Ratio-Teststatistik,χ2 = 2L1(θ1) − L0(θ0) falls M0 in M1 genestet

Logarithmus der a priori Dichteλk(θk) := logP (θk|Mk) von θk im Modell Mk

p-dimensionaler Vektor der Ableitungenλ′k(θk) von λk(θk) bzgl. der Elemente von θk

Hk := (Fk +Gk)−1 abkurzende Schreibweise

Ck := GkI −Hk(2I − FkHk)Gk abkurzende Schreibweise

pk Dimension von θk

N Gesamtstichprobengroße

Tabelle 3.2: Fur die Approximation der skalierten Bayes-Faktoren notwendige Terme und deren

Erlauterung


Approximation Ek (k = 1, 0) Ek unter θk ∼ N(ωk,Wk)

2 logB10

(3.23)≈ Ek(3.24)= 2λk(θk) + λ′k(θk)

t(Fk +Gk)−1 Ek(3.36)= log |Gk|

χ2 + (E1 − E0) 2I − Fk(Fk +Gk)−1λ′k(θk) −(θk − ωk)tCk(θk − ωk)

− log |Fk +Gk| + pk log(2π) − log |Fk +Gk|

2 logB10

(3.37)≈ Ek(3.38)= 2λk(θk) − log |Fk| + pk log(2π) Ek

(3.40)= log |Gk|

χ2 + (E1 − E0) −(θk − ωk)tGk(θk − ωk)

− log |Fk|

2 logB10(3.41)∼ Ek

(3.42)= −pk logN

χ2 + (E1 − E0)

Tabelle 3.3: Approximationen zur Berechnung des skalierten Bayes-Faktors

Man beachte, dass in der ersten Spalte der Tabelle (3.3) zweimal die Notation ≈ und einmal die

Notation ∼ die Approximation beschreibt. Allgemein ausgedruckt schreiben wir an ≈ bn, wenn

limn→∞ | an−bn |= 0 gilt. Hingegen wahlen wir die Notation an ∼ bn, wenn limn→∞(an/bn) = 1

gilt.

3.3 Anwendung auf Generalisierte Lineare Modelle

Wir betrachten folgende generalisierte lineare Modelle: Es sei yi eine Zielvariable und

xi = (xi1, . . . , xip) ein zugehoriger Vektor unabhangiger Variablen fur i = 1, . . . , n. Das GLM

M1 wird definiert, indem wir P (yi|xi,β) so spezifizieren, dass

E[yi|xi] = µi und

V ar[yi|xi] = a(φ)b′′(θ) = σ2v(µi) (3.43)

gilt mit a(φ) = σ2, Varianzfunktion b′′(θ) = v(µi) und g(µi) = xitβ mit β = (β1, . . . , βp)t.

Dabei stellt g die Linkfunktion dar. Dadurch ergibt sich die Designmatrix X als eine (n × p)-

Matrix mit den Eintragen xij , i = 1, . . . , n, j = 1, . . . , p. Ausserdem nehmen wir an, dass xi1 = 1

fur i = 1, . . . , n. Das Modell M1 enthalt also einen Intercept. Der Dispersionsparameter σ2 sei

bekannt.

Wir wollen nun das Nullmodell M0 mit dem Modell M1 vergleichen. Dabei entsteht das

Nullmodell M0 aus M1, indem wir βj = 0 setzen fur alle j = 2, . . . , p. Um diesen Vergleich

mit Hilfe der in Abschnitt 3.2 beschriebenen Approximationen fur den skalierten Bayes-Faktor

3.3. Anwendung auf Generalisierte Lineare Modelle 57

durchfuhren zu konnen, ist es allerdings noch notig, sich die Berechnung der Bestandteile der

Approximationen genau vor Augen zu fuhren und vor allem, die a priori Verteilung der Regres-

sionsparameter genau zu spezifizieren. Wir wir bereits wissen, gilt

χ2 = 2[L1(θ1,y) − L0(θ0,y)

]. (3.44)

Leider finden sich die benotigten Log-Likelihoods nicht immer im Output von Standard-Software.

Speziell die von uns verwendete Software Splus gibt die Log-Likelihoods nur in Spezialfallen aus.

Unter Poisson-Verteilungsannahme ist es allerdings moglich, χ2 mit Hilfe der zu den Modellen

M0 und M1 korrespondierenden Devianzen zu erhalten. Nach Definition (2.7) berechnet sich die

Devianz von Modell Mk, k = 1, 0, durch

Dk(y, µk) = −2a(φ) [Lk(µk,y) − Lk(y,y)] .

Dabei bezeichnet Lk(y,y) den maximalen Log-Likelihood, der im zu Modell Mk korrespondie-

renden vollen Modell erreicht wird. Unter Poisson-Verteilungsannahme erhalten wir unabhangig

vom zugrundeliegenden ModellMk, k = 0, 1, als maximalen Log-Likelihood, der im vollen Modell

erreicht wird,

L(y,y) =n∑i=1

yi + yi ln yi − ln(yi!).

Dieser Ausdruck hangt nur von den Daten ab. Bei der Berechnung von χ2 konnen wir deshalb

auf die Devianz-Differenz ubergehen:

χ2 = 2[L1(θ1,y) − L0(θ0,y)

]= 2

[L(y,y) − L0(θ0,y)

]− 2

[L(y,y) − L1(θ1,y)

]=

D0(y, θ0)a(φ)

− D1(y, θ1)a(φ)

=D0(y, θ0) −D1(y, θ1)

σ2(3.45)

Unter NB-Verteilungsannahme hangt der maximale Log-Likelihood, der im vollen Modell er-

reicht wird, nicht nur von den Daten sondern auch vom Dispersionsparameter ak des zugrunde-

liegenden Modells Mk ab. Insbesondere gilt entsprechend Gleichung (2.38)

Lk(y,y, ak) =n∑i=1

⎡⎣ y∗i∑j=0

ln(ak j + 1) − ln(yi!) + yi ln yi − (yi + a−1k ) ln(1 + akyi)

⎤⎦ , k = 0, 1.

Sind die Dispersionsparameter der Modelle M1 und M0 unterschiedlich, muss zur Berechnung

von χ2

χ2 = 2[L1(θ1,y) − L0(θ0,y)

](3.46)


benutzt werden, da eine Umformung zur Devianz-Differenz analog zu Gleichung (3.45) nicht

moglich ist. Die erwartete Fisher-Informationsmatrix fur Modell M1 lautet mit Gleichung (2.13)

F1 = σ−2XtWX,

wobei W = diag(w1, . . . , wn) mit w−1i = g′(µ(1)

i )2v(µ(1)i ). Es gilt hier µ(1)

i = g−1(xtiβ(1)

), wobei

β(1)

der MLE bedingt unter Modell M1 ist.

Entsprechend folgt fur das Nullmodell M0 (βj = 0 fur j = 2, . . . , p) mit w−1i = g′(µ(0)

i )2v(µ(0)i )

F0 = σ−2(1, . . . , 1)W (1, . . . , 1)t = σ−2(w1 + w2 + . . .+ wn)

= σ−2nwi

= σ−2n[g′(µ(0)

i )2v(µ(0)i )

]−1,

da µ(0)i = g−1(β(0)), und damit wi unabhangig von i ist. Mit β(0) bezeichnen wir den MLE be-

dingt unter dem Nullmodell M0. Die beobachtete und die erwartete Fisher-Informationsmatrix

sind gleich, wenn g die kanonische Linkfunktion ist. Somit sind die Approximationen in diesem

Fall genauer. Die Werte fur χ2, β(0), β(1), F0 und F1 konnen nun direkt in die Approximationen

fur die skalierten Bayes-Faktoren (siehe Tabelle 3.3) eingesetzt werden.

Es fehlt uns aber noch die a priori Verteilung fur die Parameter in den zu vergleichenden GLMs.

3.4 Wahl der Gestalt der a priori Verteilung fur die Parameter

in einem GLM

Wenn spezifische Informationen verfugbar sind, sollte die a priori Verteilung das naturlich reflek-

tieren. Nehmen wir nun an, es existiere a priori wenig Information. Fur diese Situation schlagt

A.E. Raftery eine sinnvolle Menge an a priori Verteilungen fur die Regressionsparameter β in

einem GLM vor. Dazu wird zunachst der Fall der linearen Modelle betrachtet. Danach werden

die Uberlegungen zuerst auf gewichtete lineare Modelle und schließlich auf GLMs ubertragen.

3.4.1 A Priori Verteilungen in Linearen Modellen

In diesem Unterabschnitt wollen wir geeignete a priori Verteilungen fur die Regressionsparameter

β = (β1, β2, . . . , βp)t im linearen Modell M1

Y = Xβ + ε, i = 1, . . . , n (3.47)

herleiten. Dabei bezeichnet X die Designmatrix und ε den zugehorigen Fehler, d.h. εi := yi−xtiβ.

Die Linkfunktion sei die Identitat, d.h. es gilt g(µ) = µ = η. Der Erwartungswert wird also

nicht transformiert. Die Varianzfunktion sei gegeben durch v(µ) = 1. Mit (3.43) folgt

V ar(εi) = V ar(yi|xi) = σ2.

3.4. Wahl der Gestalt der a priori Verteilung fur die Parameter in einem GLM 59

Wie in einem linearen Modell verlangt, weisen die Fehler im Regressionsmodell M1 danach

homogene Varianzen auf.

Zur Vereinfachung soll das Modell M1 zunachst weiteren Annahmen unterliegen. Der Response y

und die Kovariablenvektoren x1, . . . , xp seien auf empirischen Erwartungswert 0 und empirische

Varianz 1 standardisiert. Um eine ubersichtlichere Modelldarstellung zu ermoglichen, definieren

wir den standardisierten Response und die standardisierten Kovariablen folgendermaßen:

ysi :=yi − y

s0

xsij :=xij − xjsj

, i = 1, . . . , n, j = 2, . . . , p (3.48)

Dabei bezeichnet y = 1n

∑ni=1 yi den Stichprobenerwartungswert von y, xj = 1

n

∑ni=1 xij den

Stichprobenerwartungswert von xj, s20 = 1n

∑ni=1(yi − y)2 die Stichprobenvarianz von y und

s2j = 1n

∑ni=1(xij− xj)2 die Stichprobenvarianz von xj. Die Regressionsparameter fur diesen Fall

seien mit γ = (γ1, . . . , γp) bezeichnet, wobei γ1 den Interceptparameter darstellt. Ausgehend

vom linearen Modell (3.47) gilt:⎛⎜⎜⎜⎜⎜⎝

y1

y2

...

yn

⎞⎟⎟⎟⎟⎟⎠ =

⎛⎜⎜⎜⎜⎜⎝

1 x12 · · · x1p

1 x22 · · · x2p

......

. . ....

1 xn2 · · · xnp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

β1

β2

...

βp

⎞⎟⎟⎟⎟⎟⎠ +

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠

⇔

⎛⎜⎜⎜⎜⎜⎝

y1 − y

y2 − y...

yn − y

⎞⎟⎟⎟⎟⎟⎠ =

⎛⎜⎜⎜⎜⎜⎝

1 x12 · · · x1p

1 x22 · · · x2p

......

. . ....

1 xn2 · · · xnp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

β1

β2

...

βp

⎞⎟⎟⎟⎟⎟⎠ +

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠−

⎛⎜⎜⎜⎜⎜⎝

y

y...

y

⎞⎟⎟⎟⎟⎟⎠

⇔

⎛⎜⎜⎜⎜⎜⎝

y1−ys0

y2−ys0...

yn−ys0

⎞⎟⎟⎟⎟⎟⎠ =

1s0

⎛⎜⎜⎜⎜⎜⎝

1 x12 · · · x1p

1 x22 · · · x2p

......

. . ....

1 xn2 · · · xnp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

β1

β2

...

βp

⎞⎟⎟⎟⎟⎟⎠ +

1s0

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠− 1

s0

⎛⎜⎜⎜⎜⎜⎝

y

y...

y

⎞⎟⎟⎟⎟⎟⎠

⇔

⎛⎜⎜⎜⎜⎜⎝

y1−ys0

y2−ys0...

yn−ys0

⎞⎟⎟⎟⎟⎟⎠ =

1s0

⎛⎜⎜⎜⎜⎜⎝

1 x12−x2s2

· · · x1p−xp

sp

1 x22−x2s2

· · · x2p−xp

sp

......

. . ....

1 xn2−x2s2

· · · xnp−xp

sp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

β1

β2s2...

βpsp

⎞⎟⎟⎟⎟⎟⎠ +

1s0

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠− 1

s0

⎛⎜⎜⎜⎜⎜⎝

y

y...

y

⎞⎟⎟⎟⎟⎟⎠

+1s0

⎛⎜⎜⎜⎜⎜⎝

0 x2 · · · xp

0 x2 · · · xp...

0 x2 · · · xp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

β1

β2

...

βp

⎞⎟⎟⎟⎟⎟⎠


Dies lasst sich vereinfachen und mit den Definitionen in (3.48) schreiben als:⎛⎜⎜⎜⎜⎜⎝

ys1

ys2...

ys2

⎞⎟⎟⎟⎟⎟⎠ =

⎛⎜⎜⎜⎜⎜⎝

1 xs12 · · · xs1p

1 xs22 · · · xs2p...

.... . .

...

1 xsn2 · · · xsnp

⎞⎟⎟⎟⎟⎟⎠

⎛⎜⎜⎜⎜⎜⎝

1s0

(β1 − y + β2x2 + · · · + βpxp)

β2s2s0...

βpsp

s0

⎞⎟⎟⎟⎟⎟⎠ +

1s0

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠

Da die Regressionsparameter unter den kanonischen Annahmen mit γ = (γ1, . . . , γp)t bezeichnet

werden sollen, gilt⎛⎜⎜⎜⎜⎜⎝

ys1

ys2...

ys2

⎞⎟⎟⎟⎟⎟⎠ =

⎛⎜⎜⎜⎜⎜⎝

1 xs12 · · · xs1p

1 xs22 · · · xs2p...

.... . .

...

1 xsn2 · · · xsnp

⎞⎟⎟⎟⎟⎟⎠⎛⎜⎜⎝

γ1

...

γp

⎞⎟⎟⎠ +

1s0

⎛⎜⎜⎜⎜⎜⎝

ε1

ε2...

εn

⎞⎟⎟⎟⎟⎟⎠ , mit

⎛⎜⎜⎝

γ1

...

γp

⎞⎟⎟⎠ :=

⎛⎜⎜⎜⎜⎜⎝

1s0

(β1 − y + β2x2 + · · · + βpxp)

β2s2s0...

βpsp

s0

⎞⎟⎟⎟⎟⎟⎠ . (3.49)

Unter den getroffenen Annahmen lasst sich das lineare Modell M1 damit schreiben als:

Y s = Xsγ + εs, mit V ar(εsi ) =1soV ar(εi) =

σ2

s20, i = 1, . . . , n (3.50)

Dabei bezeichnet Xs die Designmatrix der Kovariablenvektoren xjs und εs := 1

s0ε = ysi −(xsi )

t γ

den zugehorigen Fehler.

Wir nehmen weiter an, dass die a priori Verteilung von (γ|M1) eine Normalverteilung ist, und

dass die Regressionsparameter (γ2, . . . , γp) a priori unabhangig sind. Die einzelnen Kovariablen

seien also fur sich genommen von Interesse. Ausserdem gehen wir davon aus, dass die a priori

Verteilung der Regressionsparameter γ fur die Testsituation ”objektiv“ ist. Im Sinne von Berger

und Sellke (1987) [Berg 87] bedeutet das, dass die a priori Verteilung symmetrisch um den

Nullwert von γ, (γ1, 0, . . . , 0)t, ist und nicht-steigend, falls ein Parameter von seinem Nullwert

abruckt.

Diese Annahmen fuhren fur die Regressionsparameter γ = (γ1, . . . , γp)t unter Modell M1 zur

a priori Verteilung

(γ|M1) ∼ N(ν, U), (3.51)

wobei ν = (ν1, 0, . . . , 0) und U = diag(ψ2, φ2, . . . , φ2). Dabei bezeichnet N(µ,Σ) eine Normal-

verteilung mit Erwartungswertvektor µ und Kovarianzmatrix Σ. Die Kovarianzmatrix U ist eine

Diagonalmatrix, weil die Regressionsparameter γ = (γ1, γ2, . . . , γp)t nach Annahme a priori un-

abhangig sind, und Unabhangigkeit Unkorreliertheit bedingt.


Die a priori Verteilung fur γ unter Modell M0 ist einfach die a priori Verteilung von γ unter

Modell M1 gegeben γ2 = . . . = γp = 0, d.h.

(γ1 |M0) ∼ N(ν1, ψ2). (3.52)

Nun gilt es, dies in eine a priori Verteilung bezuglich der Originalparameter β = (β1, β2, . . . , βp)t

zu verwandeln. Der Vektor der Regressionsparameter γ (siehe (3.49)) lasst sich implizit durch

die Gleichung

β = v +Qγ (3.53)

definieren, wobei v := (y, 0, . . . , 0) und

Q := so

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣

1 −x2/s2 −x3/s3 · · · −xp/sp0 1/s2 0 · · · 0

0 0 1/s3 · · · 0...

......

. . . . . .

0 0 0 0 1/sp

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦. (3.54)

Fur die a priori Verteilung von β unter Modell M1 gilt also wegen (γ|M1) ∼ N(ν, U):

(β|M1)(3.53)= (v +Qγ|M1) ∼ N(ν + v, QUQt) (3.55)

Die a priori Verteilung von β unter Modell M0 ist entsprechend die Verteilung unter Modell M1

gegeben β2 = . . . = βp = 0. Es gilt

(β1|M0) ∼ N(ν1 + y, ψ2s20).

3.4.2 A Priori Verteilungen in gewichteten Linearen Modellen

Lasst sich die Annahme v(µ) = 1 nicht halten, so gilt unter Modell M1

V ar(εi) = V ar[yi|xi] = σ2v(µi), i = 1, . . . , n. (3.56)

Die Fehler weisen somit keine Varianzhomogenitat mehr auf. In diesem Fall erfolgt die Schatzung

der Regressionsparameter β = (β1, β2, . . . , βp)t durch gewichtete kleinste Quadrate. Dabei wer-

den die Residuen mit wi = 1v(µi)

, i = 1, . . . , n, gewichtet. Das bedeutet, die i-te Beobachtung

erhalt mehr Gewicht, wenn V ar(yi) = σ2v(µi) klein ist. Anders ausgedruckt wird Beobachtun-

gen mit großer Fehlervarianz ein geringeres Gewicht zugewiesen. Diese Gewichtung ergibt sich

als Spezialfall der Gewichtung in GLMs. Nach Unterabschnitt 2.3.2 bzw. [McCu 89] ist in GLMs

die Schatzung der Regressionsparameter aquivalent zu gewichteten kleinsten Quadraten mit den

Gewichten wi = 1v(µi)

g′(µi)−2, i = 1, . . . , n. Durch die Annahme g(µ) = µ ergibt sich daraus

die Gewichtung wi = 1v(µi)

, i = 1, . . . , n.


Die Uberlegungen aus Unterabschnitt 3.4.1 fuhren uns wieder zur a priori Verteilung (3.55);

allerdings muss nun gewichtet werden. Es gilt

y =∑n

i=1wiyi∑ni=1wi

,

xj =∑n

i=1wixij∑ni=1wi

,

s20 =∑n

i=1wi(yi − y)2∑ni=1 wi

und

s2j =∑n

i=1wi(xij − xj)2∑ni=1wi

, i = 1, . . . , n, j = 1, . . . , p. (3.57)

3.4.3 A Priori Verteilungen in GLMs

Weisen die Fehler des Modells heterogene Varianzen auf, und sind auch andere Linkfunktioen

ausser der Identiat zugelassen, so bewegen wir uns nicht mehr in der Klasse der linearen Mo-

delle, sondern betrachten GLMs. Wie bereits im vorhergehenden Unterabschnitt erwahnt, ist in

GLMs die Schatzung der Regressionsparameter β = (β1, β2, . . . , βp)t aquivalent zu gewichteten

kleinsten Quadraten mit der adjusted dependent variable zi = g(µi) + (yi − µi)g′(µi) und Ge-

wichten wi = 1v(µi)

g′(µi)−2, i = 1, . . . , n (siehe Kapitel 2, Unterabschnitt 2.3.2 bzw. [McCu 89]).

Wiederum fuhren uns die Uberlegungen aus Unterabschnitt 3.4.1 zur a priori Verteilung (3.55).

Es wird nur y durch z ersetzt, und alle ”summary-Kenngroßen“ werden mit w = 1v(µ)

g′(µ)−2

gewichtet. In der a priori Verteilung (3.55) gilt nun

v = (z, 0, . . . , 0)t, wobei z =∑n

i=1 wizi∑ni=1 wi

.

Fur die Stichprobenvarianzen in der Matrix Q (3.54) gilt

s20 =∑n

i=1wi(zi − z)2∑ni=1wi

=∑n

i=1wi(z2i − 2ziz + z2)∑ni=1 wi

=∑n

i=1wiz2i − 2z

∑ni=1 wizi + z2

∑ni=1 wi∑n

i=1wi

=∑n

i=1wiz2i∑n

i=1wi− 2

∑ni=1 wizi∑ni=1 wi

∑ni=1wizi∑n

i=1 wi+

[∑ni=1 wizi∑ni=1 wi

]2 ∑ni=1 wi∑n

i=1wi

=∑n

i=1wiz2i∑n

i=1wi− 2

[∑ni=1wizi∑ni=1 wi

]2

+[∑n

i=1wizi∑ni=1wi

]2

=∑n

i=1wiz2i∑n

i=1wi−[∑n

i=1 wizi∑ni=1wi

]2

, und mit (3.58)

xj =∑n

i=1wixij∑ni=1 wi

folgt analog

s2j =∑n

i=1wi(xij − xj)2∑ni=1 wi

=

∑ni=1wix

2ij∑n

i=1 wi−[∑n

i=1 wixij∑ni=1 wi

]2

. (3.59)


Werden mehrere Modelle betrachtet, ist es wunschenswert, dass die a priori Verteilungen unter-

einander konsistent sind. Wenn also das Modell M0 durch Setzen der Restriktionen ρ(β) = 0

auf Parameter des Modells M1 definiert ist, soll P (β|M1) = P (β|M0, ρ(β) = 0) gelten. Um

dies sicherzustellen, empfiehlt es sich, eine a priori Verteilung fur das großte Modell wie oben

gezeigt zu bestimmen, und dann die a priori Verteilungen fur die anderen Modelle herzuleiten,

indem man entsprechend der modelldefinierenden Annahmen bedingt. Angenommen die a priori

Verteilung (3.55) korrespondiert zum großten Modell MK und Mk ist definiert durch Nullsetzen

verschiedener βj in MK . Dann gilt (β|Mk) ∼ N(ν [k] + v[k], Q[k]U [k]Q[k]t). Dabei bezeichnet der

Exponent [k] die Elemente von Vektoren bzw. die Zeilen und die Spalten von Matrizen, die zu

den in Modell Mk nullgesetzten Parametern von MK korrespondieren und somit entfernt worden

sind.

3.4.4 Wahl der Hyper-Parameter in der a priori Verteilung der Regressions-

parameter in GLMs

Die a priori Verteilung (3.55) hat drei benutzerspezifische Parameter, ν1, ψ und φ. A.E. Raftery

zeigt in seinem Artikel [Raft 94] (Seite 9ff), welche Werte diese Parameter annehmen sollten,

um angemessen die Situation zu reprasentieren, wenn a priori wenig Information vorhanden

ist. Ausgehend von den Gleichungen (3.37) und (3.38) erhalten wir eine weitere approximative

Darstellung des Bayes-Faktors:

2 logB10 ≈ χ2 + (E∗1 − E∗

0)

= 2L1(θ1) −L0(θ0)− log |F1| + 2λ1(θ1) + p1 log(2π) + log |F0| − 2λ0(θ0) − p0 log(2π)

= 2

log(P (D|θ1,M1)) − log(P (D|θ0,M0))

− log |F1| + 2 log(P (θ1|M1)) + p1 log(2π) + log |F0| − 2 log(P (θ0|M0))

−p0 log(2π)

logB10 ≈ log P (D|θ1,M1) − log P (D|θ0,M0)

−12

log |F1| + log(P (θ1|M1)) +12p1 log(2π)

+12

log |F0| − log(P (θ0|M0)) − 12p0 log(2π)

= log

P (D|θ1,M1)P (D|θ0,M0)

+

12

log|F0||F1| + log

P (θ1|M1)P (θ0|M0)

+(p1 − p0

2

)log(2π)

= log

P (D|θ1,M1)P (D|θ0,M0)

[ |F0||F1|

]1/2 P (θ1|M1)P (θ0|M0)

(2π)(p1−p0)/2


Beidseitige Anwendung der Exponentialfunktion bringt uns zu:

B10 ≈ (2π)(p1−p0)/2[P (D|θ1,M1)P (D|θ0,M0)

][ |F0||F1|

] 12

[P (θ1 |M1)P (θ0 |M0)

](3.60)

Diese Darstellung hat den Vorteil, dass die a priori Verteilung des Parametervektors θ nur im

letzten Faktor auftritt. Diesen letzten Faktor bezeichnet man als ratio of prior ordinates (RPO)

an der Stelle der MLE’s.

Zur Vereinfachung diskutieren wir die Wahl der Hyper-Parameter wieder unter Annahme

der kanonischen Situation des Unterabschnitts 3.4.1, d.h. es sei wieder g(µ) = µ und v(µ) = 1.

Die Variablen seien auf empirischen Erwartungswert 0 und empirische Varianz 1 standardisiert.

Das Modell M1 enthalte nur eine unabhangige Variable x2 = (x12, x22, . . . , xn2)t. Die zwei zu

vergleichenden Modelle lauten damit

M1 : Yi = β1 + β2xi2 + εi und M0 : Yi = β1 + εi, i = 1, . . . , n

mit

y =1n

n∑i=1

yi = 0,

x2 =1n

n∑i=1

xi2 = 0,

s20 =1n

n∑i=1

(yi − y)2 =1n

n∑i=1

y2i = 1 und

s22 =1n

n∑i=1

(xi2 − xj)2 =1n

n∑i=1

x2i2 = 1. (3.61)

Mit θ1 = (β1, β2)t und θ0 = β1 und wegen (θ1|M1)(3.55)∼ N((ν1 + y, 0), QUQt) mit

QUQt =

(1 −x2

s2

0 1s2

)(ψ2 0

0 φ2

)(1 0

−x2s2

1s2

)x2=0,s2=1

=

(ψ2 0

0 φ2

)

gilt

RPO10(φ; β2) :=P (θ1|M1)

P (θ0|M0)=

12πψφ exp

−1

2(β1 − ν1, β2)

(1ψ2 0

0 1φ2

)(β1 − ν1

β2

)

1√2πψ

exp−1

2(β1−ν1)2

ψ2

=1√2πφ

exp−1

2

((β1−ν1)2

ψ2 + β22φ2

)exp

−1

2(β1−ν1)2

ψ2

=

1φ√

2πexp

(− β2

2

2φ2

). (3.62)


Dabei wurde benutzt, dass unter den getroffenen Annahmen β1 = 0 in beiden Modellen gilt. Da

es sich bei den von uns betrachteten Modellen M0 und M1 aufgrund der kanonischen Annahmen

um lineare Modelle handelt, wird der Kleinste-Quadrate-Ansatz (Least Squares) zur Schatzung

der Parameter herangezogen. Dabei minimiert man die Summe der quadratischen Abweichungen

(Residual Sum of Squares)

Q(β1, β2) =n∑i=1

(yi − β1 − β2xi2)2.

Damit ist (β1, β2) durch Q(β1, β2) = minβ1,β2Q(β1, β2) definiert. Partielles Ableiten der Residu-

enquadratsumme Q nach β1 und β2 und anschließendes Nullsetzen liefert unter Beachtung der

getroffenen Annahmen:

∂Q

∂β1(β1, β2) = −2

n∑i=1

(yi − β1 − β2xi2)!= 0

⇔ 0 =n∑i=1

yi − nβ1 − β2

n∑i=1

xi2

⇔ nβ1 =n∑i=1

yi − β2

n∑i=1

xi2

⇔ β1 =1n

n∑i=1

yi − β21n

n∑i=1

xi2

= y − β2x2x2=y=0

= 0 (3.63)

∂Q

∂β2(β1, β2) = −2

n∑i=1

(yi − β1 − β2xi2)xi2!= 0

⇔ 0 =n∑i=1

yixi2 − β1

n∑i=1

xi2 − β2

n∑i=1

x2i2

⇔ β2

n∑i=1

x2i2 =

n∑i=1

yixi2 − β1

n∑i=1

xi2

⇔ β2 =∑n

i=1 yixi2 − β1∑n

i=1 xi2∑ni=1 x

2i2

=1n

∑ni=1 yixi2 − β1

1n

∑ni=1 xi2

1n

∑ni=1 x

2i2

1n

∑x2

i2=1=

1n

n∑i=1

yixi2 − β1x2x2=0=

1n

n∑i=1

yixi2 (3.64)

Bis zu diesem Approximationslevel (Gleichung (3.60) mit Gleichung (3.62)) haben ν1 und ψ kei-

nen Einfluss auf den Bayes-Faktor B10. Ausserdem wissen wir aus [Raft 94], dass numerischen

Experimenten zufolge die exakten Bayes-Faktoren ziemlich - wenn auch nicht komplett - insen-

sitiv gegenuber der genauen Wahl der Werte von ν1 und ψ sind. A.E. Raftery wahlte ν1 = 0

und ψ = 1 (siehe [Raft 94], Seite 9). Scheinbar besteht nur wenig Bedarf, Wertebereiche fur ν1

und ψ zu betrachten.


Deshalb untersuchen wir nun die Wahl von φ. Der Cauchy-Schwarz-Ungleichung (CSU) zufolge

gilt |β2| ≤ 1 :

n|β2| (3.64)= |

n∑i=1

yixi2| = |ytx2|CSU≤ ||y||||x2|| =

√√√√ n∑i=1

y2i

√√√√ n∑i=1

x2i2

(3.61)=

√n√n = n

⇒ |β2| ≤ 1 (3.65)

Dabei bezeichnt ||.|| die euklidische Norm. In einer ersten Uberlegung suchen wir ein φ, so dass

RPO10(φ; β2) fur alle moglichen Werte von β2 so nahe wie moglich bei 1 liegt. Damit soll der

Effekt der a priori Verteilung auf den Bayes-Faktor minimiert werden. Fur alle (positiven) Werte

von φ ist RPO10(φ; β2) bezuglich |β2| ≤ 1 minimal, wenn | β2 |= 1, da das Maximum fur β2 = 0

erreicht wird, und RPO10(φ; β2) fur β2 gegen plus und minus unendlich von oben gegen Null

geht. Dies kann man wie folgt sehen:

∂

∂β2

RPO10(φ; β2) =1

φ√

2πexp

(− β2

2

2φ2

)(− β2

φ2

)!= 0

⇒ β2 = 0 stationarer Punkt,

∂2

∂β2∂β2

RPO10(φ; β2)|β2=0 =1

φ√

2πexp

(− β2

2

2φ2

)(− β2

φ2

)(− β2

φ2

)

+1

φ√

2πexp

(− β2

2

2φ2

)(− 1φ2

)|β2=0 < 0

⇒ Maximum bei β2 = 0.

limβ2→∞

RPO10(φ; β2) = limβ2→∞

1φ√

2πexp

(− β2

2

2φ2

)= 0

limβ2→−∞

RPO10(φ; β2) = limβ2→−∞

1φ√

2πexp

(− β2

2

2φ2

)= 0

Mit Abbildung 3.1 soll dies veranschaulicht werden. Fur |β2| = 1 ist RPO10(φ; β2) immer klei-

ner als 1 (siehe Abbildung 3.2). Wir wahlen also |β2| = 1 und versuchen sicherzustellen, dass

RPO10(φ; 1) nicht zu klein ist. RPO10(φ; 1) ist wiederum maximal bezuglich φ fur φ = 1. Dies

lasst sich wie folgt zeigen.

∂

∂φRPO10(φ; 1) =

1√2π

exp(− 1

2φ2

)(φ−4 − φ−3

) != 0

φ positiv⇒ φ = 1 stationarer Punkt,∂2

∂φ∂φRPO10(φ; 1) =

1√2π

exp(− 12φ2

)(−4φ−5 − 3φ−4) < 0

φ positiv⇒ Maximum beiφ = 1.


00.1

0.2

−0.1−0.08−0.06−0.04−0.0200.020.040.060.080.10

5

10

20

30

40

phibeta2dach

RP

O(p

hi;b

eta2

dach

)

Abbildung 3.1: RPO10(φ; β2)

limφ→∞

RPO10(φ; 1) = limφ→∞

1φ√

2π︸︷︷︸→0

exp(− 1

2φ2

)︸︷︷︸

→1

= 0 (3.66)

Ausserdem gilt mit der Regel von L’Hospital (L’H):

limφ→0

RPO10(φ; 1) = limφ→0

1φ√

2π

exp(

12φ2

) L′H= limφ→0

− 1√2πφ−2

exp(

12φ2

)12(−2)φ−3

= limφ→0

1√2π

1

exp(

12φ2

)︸︷︷︸

→∞

φ−1︸︷︷︸→∞

= 0

Im Punkt φ = 1 gilt RPO10(φ; 1) = (2πe)−12 = 0.24 (siehe Abbildung 3.2). Die Tatsache, dass

dieser Wert kleiner als 1 ist, reflektiert die ”Strafe“ fur den Mangel an Sparsamkeit, die die

Bayes-Prozedur automatisch dem großeren Modell M1 auferlegt. Fur steigendes φ (φ > 1) sinkt

die Evidenz fur das Modell M1 weiter (siehe (3.66) und Abbildung 3.2).

Um das Ausmaß, mit dem ein gegebenes φ schlechter ist als φ = 1, mit den Bezeichnungen aus

dieser ersten Uberlegung zu messen, definierte A. E. Raftery den Wert

R(φ) :=RPO10(1; 1)RPO10(φ; 1)

=1√2π

exp(−12)

1φ√

2πexp(− 1

2φ2 )= φ exp(φ−2 − 1)/2.


0 0.2 0.4 0.6 0.8 1 1.2 1.40

0.05

0.1

0.15

0.2

0.25

phi

RP

O(p

hi;1

)

Abbildung 3.2: RPO10(φ; 1)

R(φ) misst den maximalen Einfluss der a priori Verteilung auf den Bayes-Faktor B10. Wir wollen

alle Werte von φ ausschließen, fur die R(φ) zu groß wird, und beschranken deshalb R(φ) nach

oben durch R(φ) < c. Dabei reprasentiert c ”wenig Evidenz“. A.E. Raftery wahlte in Gedenken

an H. Jeffreys [Jeff 61] c =√

10 = 3.16. Dieser Wert korrespondiert zur Evidenz ”nicht mehr

wert als bloße Erwahnung“ (vgl. Tabelle 3.1). Es wird also vorgeschlagen, φ so zu wahlen, dass

R(φ) <√

10 gilt.

R(φ) nimmt sein Minimum bei φ = 1 an (R(1) = 1) und geht fur φ gegen Null und fur φ gegen

unendlich gegen unendlich. Dies lasst sich wie folgt sehen:

limφ→∞

R(φ) = limφ→∞

φ︸︷︷︸→∞

exp(φ−2 − 1)/2

︸︷︷︸→exp(− 1

2)

= ∞

Mit der Regel von L’Hospital (L’H) gilt:

limφ→0

R(φ) = limφ→0

φ exp(φ−2 − 1)/2

= lim

φ→0

exp(φ−2 − 1)/2

φ−1

(L′H)= lim

φ→0exp

(φ−2 − 1)/2

︸︷︷︸→∞

φ−1︸︷︷︸→∞

= ∞


∂

∂φR(φ) = exp

(φ−2 − 1)/2

+ φ exp

(φ−2 − 1)/2

12(−2)φ−3 != 0

⇔ exp(φ−2 − 1)/2

(1 − φ−2) = 0

⇔ φ = 1

⇒ φ = 1 stationarer Punkt

∂2

∂φ2R(φ)|φ=1 = exp

(φ−2 − 1)/2

12(−2)φ−3(1 − φ−2) + exp

(φ−2 − 1)/2

2φ−3|φ=1

= exp(φ−2 − 1)/2

(φ−5 + φ−3)|φ=1 = 2 > 0

⇒ Minimum beiφ = 1

Das Kriterium R(φ) <√

10 = 3.16 schließt danach sowohl φ-Werte, die viel großer als 1 sind,

als auch φ-Werte, die viel kleiner als 1 sind, aus.

Eine zweite Uberlegung beschaftigt sich mit dem Vergleich nicht-genesteter Modelle M1

und M2. Wir betrachten den Fall, bei welchem Modell M1 und Modell M2 jeweils eine un-

abhangige Variable x2 = (x12, . . . , xn2)t bzw. x3 = (x13, . . . , xn3)t besitzen. Die zu vergleichen-

den Modelle lauten damit

M1 : Yi = β1 + β2xi2 + εi und M2 : Yi = β1 + β3xi3 + εi, i = 1, . . . , n.

Es gilt mit θ1 = (β1, β2)t und θ2 = (β1, β3)t

RPO12(φ; β2, β3) =P (θ1 |M1)P (θ2 |M2)

=

P (θ1|M1)

P (θ0|M0)

P (θ2|M2)

P (θ0|M0)

(3.62)=

1φ√

2πexp(− β2

22φ2 )

1φ√

2πexp(− β2

32φ2 )

= exp

(− β2

2

2φ2

)exp

(β2

3

2φ2

)

= exp

[−(β2

2 − β23)

2φ2

]. (3.67)

Wie vorher sichert die Cauchy-Schwarz-Ungleichung, dass |β2| ≤ 1 und |β3| ≤ 1 (siehe Unglei-

chung (3.65)). Wir suchen nun wieder das φ, fur das RPO12(φ; β2, β3) uber dem Wertebereich

von β2 und β3 moglichst nahe bei 1 liegt, um den Einfluss der a priori Verteilung zu minimieren.

RPO12(φ; β2, β3) ist am weitesten von 1 entfernt, wenn β2 = 0 und β3 = 1. Deshalb definierte

A.E. Raftery den Wert:

S(φ) := RPO12(φ; 0, 1) = exp(φ−2/2) (3.68)

Wir fordern wieder S(φ) < c =√

10. Man bemerke, dass S(φ) ↓ 1 fur φ→ ∞. Somit wird dieses

Kriterium kleine Werte von φ ausschließen, große hingegen nicht.


Hat man einen einzigen Wert fur φ zu wahlen, bietet sich φ = e1/2 = 1.65 an, da fur diesen Wert

R(φ) = S(φ) gilt.

φ exp(φ−2 − 1)/2 = exp(φ−2/2)

⇔ log φ+ (φ−2 − 1)/2 = φ−2/2

⇔ 2 log φ+ φ−2 − 1 = φ−2

⇔ log φ = 1/2

⇔ φ = e1/2 ≈ 1.65

Zur Veranschaulichung betrachte man die Abbildung 3.3.

phi

R(ph

i) bzw

. S(ph

i)

1 2 3 4 5 6

01

23

45

phi

R(ph

i) bzw

. S(ph

i)

1 2 3 4 5 6

01

23

45

R(phi)S(phi)

Abbildung 3.3: R(φ) und S(φ) als Funktionen von φ

Im Allgemeinen ist es allerdings besser, die Ergebnisse fur einen vernunftigen Wertebereich von

φ zu betrachten. Es gilt R(φ) < c =√

10 und gleichzeitig S(φ) < c =√

10 fur 0.67 ≤ φ ≤ 5.10:

S(φ) = exp(

12φ2

)<

√10

⇔ 12φ2

< ln(√

10)

⇔ φ >

√1

2 ln√

10= 0.66

3.5. Unbekannter Dispersionsparameter und Uberdispersion 71

R(φ) = φ exp(φ−2 − 1)/2

<

√10

mathematica⇔ 0.46 ≤ φ ≤ 5.10

Wie in Abbildung 3.3 zu erkennen ist, steigen sowohl R(φ) als auch S(φ) bereits fur φ < 1

extrem an. Deshalb ist es nicht notig, Werte kleiner 1 fur φ zu untersuchen. Die Uberlegungen

dieses Unterabschnitts zur Wahl der Hyperparameter lassen sich auch auf andere Link- und

Varianzfunktionen ubertragen. Dementsprechend legte A.E. Raftery folgenden Wertebereich fur

φ als den relevanten fest:

1 ≤ φ ≤ 5 mit φ = 1.65 als Zentralwert (3.69)

Der Log-Bayes-Faktor als Funktion von φ andert sich rapide fur φ < 1 und sehr viel langsamer

fur φ im von A.E. Raftery hergeleiteten Wertebereich (3.69).

3.5 Unbekannter Dispersionsparameter und Uberdispersion

Bis jetzt haben wir angenommen, dass der Dispersionsparameter σ2 bekannt ist. Dies ist der

Fall in Poisson- und Binomialmodellen ohne Uberdispersion. Wenn σ2 unbekannt ist, ware ein

offensichtlicher Losungsweg, zu verfahren wie oben und nur σ2 durch σ2 zu ersetzen, wie es

in [McCu 89] vorgeschlagen wird. Ein vernunftiger Schatzer ware σ2 = P/(n − p), wobei P

die Pearson-Goodness-of-Fit-Statistik fur das komplexeste unter allen betrachteten Modellen ist

(siehe [McCu 89], Seite 91 und 127).

Ein genauerer und komplett bayesianischer Ansatz besteht allerdings darin, den Dispersionspa-

rameter σ2 in gleicher Weise wie die Regressionsparameter β = (β1, . . . , βp)t als einen Parameter

zu betrachten, indem man ihm eine a priori Verteilung gibt und ausintegriert.

Dazu konstruieren wir sogenannte Durchschnitts-Bayes-Faktoren (overall Bayes factors) zum

Vergleich von Modell M1 mit Modell M0. Ist in beiden Modellen der Dispersionsparameter

unbekannt, so sind die Durchschnitts-Bayes-Faktoren gegeben durch

B10 :=∫ ∫

B10(σ21 , σ

20)π(σ2

1)π(σ20)dσ

21dσ

20 . (3.70)

Dabei bezeichnet π(σ2i ), i = 1, 0, jeweils die a priori Dichte der zu den betrachteten Model-

len gehorigen Dispersionsparameter σ21 bzw. σ2

0 . Im Gegensatz zu den Bayes-Faktoren B10, die

sich nur fur individuell gewahlte Dispersionsparameter berechnen lassen, berucksichtigen die

Durchschnitts-Bayes-Faktoren, dass die Dispersionsparameter unbekannt sind und beziehen so-

mit die Unsicherheit in der Schatzung der Dispersionsparameter mit ein.

Ist der Dispersionsparameter im Modell M1 unbekannt wir zum Beispiel in einem NB-Modell,

im Modell M0 hingegen bekannt wie zum Beispiel in einem Poissonmodell, so vereinfacht sich


die Berechnung der Durchschnitts-Bayes-Faktoren zu:

B10 =∫B10(σ2

1)π(σ21)dσ

21 (3.71)

Dabei bezeichnet π(σ21) die a priori Dichte des Dispersionsparameters im Modell M1. Vergleicht

man ein NB-Modell M1 mit einem Poissonmodell M0, und wird der Dispersionsparameter im

NB-Modell M1 nicht im Voraus gewahlt, sondern stattdessen geschatzt, so tendieren die Bayes-

Faktoren dazu, die Evidenz fur des NB-Modell M1 zu uberschatzen. In diesem Falle sind den

Bayes-Faktoren B10 die Durchschnitts-Bayes-Faktoren B10 aus den oben genannten Grunden

vorzuziehen.

Kapitel 4

Anwendungsbeispiele

In diesem Kapitel wollen wir nun fur konkrete Datenbeispiele Bayes-Faktoren zum Vergleich

genesteter und vor allem auch nicht-genesteter Regressionsmodelle berechnen. Dazu benutzen

wir die Approximation (3.23) unter der Annahme, dass die a priori Verteilung der Regressions-

parameter die Normalverteilung ist. Danach handelt es sich in den Ausfuhrungen dieses Kapitels

immer um approximierte Bayes-Faktoren, wenn von Bayes-Faktoren die Rede ist.

4.1 Patent-Daten

Im Folgenden werden die im Internet auf der Seite

http://coe.ubc.ca/users/marty/patent.data

bereitgestellten Patent-Daten untersucht. Diese Daten stammen von 70 amerikanischen High-

Tech-Firmen und wurden im Jahre 1976 erhoben. Drei Spalten sind fur unsere Analyse von

Interesse. Zum einen wahlen wir die Spalte patente, die uns die Anzahl der angemeldeten

Patente der Unternehmen angibt, als Response. Als erste Regressorvariable benutzen wir die

Spalte rdsales. Sie gibt uns Aufschluss uber das Verhaltnis der Geldbetrage, die in Forschung

und Entwicklung (research and development) investiert wurden, zu den durch die Entwicklung

der Patente erzielten Einnahmen (R&D/Sales). Die Spalte mit dem Namen log(rd) dient uns als

Basis fur unsere zweite Regressorvariable. Sie gibt uns den Logarithmus der fur Forschung und

Entwicklung benutzten Gelder an. Uns interessieren aber die reinen Investitionsgelder. Deshalb

transformieren wir diese Spalte mittels der Exponentialfunktion zuruck und verwenden die somit

erhalten Spalte rd als zweite Regressorvariable.

4.1.1 Explorative Datenanalyse der Patent-Daten

Die Tatsache, dass alle Kovariablen metrisch sind, und eine Gewichtung des Response nicht

notwendig ist, erleichtert uns die Arbeit. Bei der Analyse der Zahlvariable patente konnen wir

73

74 Kapitel 4. Anwendungsbeispiele

somit auf eine Faktorisierung der Regressoren verzichten, und auch die Anwendung eines offsets

entfallt. Alle Responsewerte haben Gewicht 1. Nehmen wir den Fehler als Poisson-verteilt an,

und wahlen wir die Log-Linkfunktion als Linkspezifizierung, so folgt mit Abschnitt 2.2.3:

patente ∼ Poi(µi), i = 1, . . . , 70, mit

µi = ti exp(xtiβ) = exp(ln(ti) + xt

iβ) bzw. ln(µi) = ln(ti)︸︷︷︸offset

+xtiβ

Fur die von uns betrachteten Patent-Daten gilt ti = 1, i = 1, . . . , n.

Plottet man also die Kovariablen gegen den Logarithmus der Patente, sollte ein linearer Zu-

sammenhang zu erkennen sein. Wie aus der explorativen Analyse der Patent-Daten in [Siko 02]

(Abbildung 5.9, Seite 134) ersichtlich, verbessert sich die Linearitat zwischen der Kovariablen

rd und dem logarithmierten Response durch Transformation mittels der 5. Wurzel. Die transfor-

mierte Kovariable 5. Wurzel aus rd bezeichnen wir mit rdw5. Die Kovariable rdsales weist nur

einen undeutlichen linearen Zusammenhang mit dem logarithmierten Response auf. Trotzdem

wollen wir sie als Modellvariable aufnehmen.

Tragt man die Interaktion zwischen rdsales und rdw5 als dritte Kovariable gegen den Lo-

garithmus des Patente ab, so ist auch diese als relevant einzustufen, da ein leichter linearer

Zusammenhang mit dem Logarithmus der Patente erkennbar ist.

4.1.2 Modellvergleiche mit Hilfe der Bayes-Faktoren

Um einerseits Aussagen uber die Signifikanz der Interaktion zwischen den Kovariablen rdsales

und rdw5 zu treffen und andererseits beurteilen zu konnen, ob die von uns betrachteten Patent-

Daten Uberdispersion aufweisen, werden nun die vier Regressionsmodelle in Tabelle 4.1 mit

Hilfe der Bayes-Faktoren hinsichtlich ihrer Anpassungsgute verglichen. In dieser Tabelle bezeich-

Modell Fehlerverteilung Kovariablen

1 Poisson rdsales,rdw5

2 Poisson rdsales,rdw5,rdsales*rdw5

3 Negbin rdsales,rdw5

4 Negbin rdsales,rdw5,rdsales*rdw5

Tabelle 4.1: Zu vergleichenden Modelle zur Anpassung der Patent-Daten

net rdsales*rdw5 die Interaktion zwischen den Kovariablen rdsales und rdw5. Als Response

wahlen wir in jedem Modell die Anzahl der angemeldeten Patente (patente). Die Linkfunktion

sei der Logarithmus.

Im Folgenden bezeichnen wir mit 2logBij(ai, aj , φ) die Approximation (3.23) des skalierten

4.1. Patent-Daten 75

Bayes-Faktors unter der Annahme, dass die a priori Verteilung der Regressionsparameter die

Normalverteilung ist mit Hyperparameter φ. Dabei ist Modell M1 gleich Modell Mi mit Dis-

persionsparameter ai, und Modell M2 ist gleich Modell Mj mit Dispersionsparameter aj . In

Poissonmodellen hat der Dispersionsparameter a den Wert Null.

Um die Signifikanz der Interaktion zu untersuchen, betrachten wir zunachst den Vergleich der

Poissonmodelle 1 und 2. Fur die Approximation des skalierten Bayes-Faktors zum Vergleich

zweier Poissonmodelle wurde die Funktion BFpois2ln (siehe Anhang, Seite 90ff) geschrieben.

Dieser Funktion sind der gewahlte Wert fur den Hyperparameter φ der a priori Verteilung, die

zu vergleichenden Modelle, der Response und der Offset zu ubergeben. Dabei soll der Hyperpa-

rameter φ nach den Ausfuhrungen in Unterabschnitt 3.4.4 im Wertebereich (3.69) liegen. Um

auch die Auswirkung der Wahl von φ auf die Hohe des Bayes-Faktors untersuchen zu konnen,

berechnen wir die Bayes-Faktoren fur die Modellvergleiche immer fur φ = 1, φ = 1.65 und φ = 5.

Es ist zu beachten, dass der Funktion BFpois2ln auch dann ein Offset zu ubergeben ist, wenn

alle Responsewerte Gewicht 1 haben. In diesem Fall wahlt man als Offset den Logarithmus eines

Eins-Vektors der Lange des Datensatzes. Beispielsweise ist fur die Berechnung des skalierten

Bayes-Faktors zum Vergleich der Poissonmodelle 1 und 2 mit φ = 1 folgende Splus-Eingabe

notwendig:

patente21<-BFpois2ln(1,patente∼offset(log(rep(1,70)))+rdsales+rdw5+rdsales*rdw5,patente∼offset(log(rep(1,70)))+rdsales+rdw5,patente,

log(rep(1,70))

Mit φ = (1, 1.65, 5) erhalten wir fur den skalierten Bayes-Faktor die Werte

2logB21(0, 0, φ) = (7.43, 5.75, 5.69). (4.1)

Diese Werte sprechen nach Tabelle 3.1 stark fur das Poissonmodell 2. Die Hinzunahme der Inter-

aktion zwischen den Kovariablen rdsales und rdw5 verbessert somit die Modellanpassung. Mit

steigendem φ sinkt die Evidenz fur das Poissonmodell 2. Da die Bayes-Prozedur kleinere Modelle

bevorzugt, wird das Poissonmodell 2 dafur bestraft, dass es im Vergleich zum Poissonmodell 1

eine Kovariable mehr enthalt (vgl. Kapitel 3, Seite 67). Die Wahl von φ im Wertebereich (3.69)

hat aber wie gewunscht so geringfugig Einfluss auf die Hohe des skalierten Bayes-Faktors, dass

die Evidenzaussage davon nicht beeinflusst wird.

Da die Poissonmodelle 1 und 2 genestet sind, lassen sie sich auch mit herkommlichen, auf

p-Werten basierenden Signifikanztests vergleichen. Dafur berechnen wir die Poissonmodelle 1

und 2 in Splus:

patent1<-glm(patente∼rdsales+rdw5,family=poisson,link=log,x=T) (4.2)


patent2<-glm(patente∼rdsales+rdw5+rdsales*rdw5,family=poisson,link=log,x=T) (4.3)

Zu diesen Modellen konnen wir uns nun das jeweilge Splus-summary ausgeben lassen (siehe Ta-

belle 4.2 und Tabelle 4.3). In der zweiten Spalte Value lassen sich die geschatzten Werte der

Coefficients:

Value Std. Error t value Pr(>|t|)

(Intercept) -1.585 0.1272 -12.46 0

rdsales 0.421 0.0496 8.49 0

rdw5 2.564 0.0566 45.29 0

(Dispersion Parameter for Poisson family taken to be 1 )

Null Deviance: 3325 on 69 degrees of freedom

Residual Deviance: 387 on 67 degrees of freedom

Number of Fisher Scoring Iterations: 4

Tabelle 4.2: Splus-summary des Poissonmodells 1 fur die Patent-Daten

Coefficients:

Value Std. Error t value P(>|t|)

(Intercept) -1.17 0.185 -6.30 2.98e-010

rdsales -5.50 1.991 -2.76 0.006

rdw5 2.30 0.104 22.20 0

rdsales:rdw5 3.90 1.312 2.98 0.003

(Dispersion Parameter for Poisson family taken to be 1 )


Residual Deviance: 377 on 66 degrees of freedom


Tabelle 4.3: Splus-summary des Poissonmodells 2 fur die Patent-Daten

Regressionsparameter β ablesen. Betrachten wir das summary fur das Poissonmodell 2 (siehe

Tabelle 4.3), so lasst sich aus den zum Wald-Test (t-Test, siehe Kapitel 2, Seite 21) korrespon-

dierenden p-Werten in der letzten Spalte schließen, dass alle Kovariablen des Poissonmodells 2

hoch signifikant sind. Der p-Wert gibt uns das kleinste α-Niveau an, fur das die Hypothese H0

gerade noch verworfen werden kann. Die p-Werte in Tabelle 4.3 sind sehr viel kleiner als das

Standard-Signifikanzniveau 0.05. Das heisst, fur j = 1, . . . , 4 kann die Hypothese Hj : βj = 0 zu

einem sehr kleinen Signifikanzniveau verworfen werden. Die Kovariablen haben einen signifikan-


ten Einfluss auf die Zielvariable patente. Mit einem p-Wert von 0.003 ist insbesondere auch die

Interaktion hoch signifikant und sollte nicht aus dem Modell entfernt werden.

Auch der Partial-Deviance-Test (siehe Kapitel 2, Seite 25) kann zur Uberprufung der Signi-

fikanz der Interaktion benutzt werden. Dieser benotigt die Devianzen der Poissonmodelle 1 und

2, die den summaries (siehe Tabelle 4.2 und Tabelle 4.3) entnommen werden konnen. Basierend

auf der Partial-Deviance-Statistik mit einem Freiheitsgrad erhalt man einen p-Wert von 0.002.

Damit kann auch bei diesem Test die Nullhypothese H0 : β Interaktion = 0 zu einem sehr klei-

nen Signifikanzniveau verworfen werden. Die Interaktion ist hoch signifikant.

Das von den skalierten Bayes-Faktoren gelieferte Ergebnis bestatigt sich damit voll und ganz.

Das Poissonmodell 2 liefert eine bessere Anpassung als das Poissonmodell 1.

Betrachten wir noch einmal das Splus-summary fur das Poissonmodell 2 (siehe Tabelle 4.3),

stellen wir allerdings fest, dass trotz der hohen Signifikanz der Kovariablen die Anpassung der

Daten durch dieses Modell als weniger gelungen einzustufen ist. Die Devianz ist mit 377 viel

großer als die Anzahl der Freiheitsgrade (66). Die mangelnde Anpassungsgute kann nun meh-

rere Grunde haben. Sie kann z.B duch Ausreisser, Linkmissspezifizierung oder Uberdispersion

in den Daten entstehen. Eine bereits in [Siko 02] durchgefuhrte Residuenanalyse bestatigt, dass

das Poissonmodell 2 eine eher schlechte Anpassung an die Daten liefert. Eine Linkmissspezi-

fizierung kann allerdings ausgeschlossen werden. Es liegt somit die Vermutung nahe, dass die

Patent-Daten Uberdispersion aufweisen. In Abbildung 4.1 wurde deshalb versucht, ein mogliches

Vorhandensein von Uberdispersion in den Patent-Daten grafisch dazustellen. Dazu haben wir die

geschatzten Erwartungswerte exp(xtiβ) gegen die geschatzten Varianzen abgetragen. Dabei wur-

den als Schatzer fur die Varianzen die quadrierten Rohresiduen (yi− µi)2, i=1,. . . ,n, verwendet.

Entsprechend der Aquidispersionsannahme im Poissonmodell (E(Yi) = V ar(Yi), i=1,. . . ,n) soll-

ten alle Werte in der Grafik auf der als durchgezogene Linie eingezeichneten Winkelhalbierenden

liegen. Die mit dem Splus-Befehl lowess erzeugte Glattungsfunktion (gestrichelte Linie) verlauft

allerdings weit oberhalb der Winkelhalbierenden. Die Varianz der einzelnen Beobachtungen ist

deutlich großer als ihr Erwartungswert. Eine entscheidende Modellannahme der Poissonvertei-

lung ist damit verletzt.

Deshalb vergleichen wir nun die NB-Modelle 3 und 4 mit den Poissonmodellen 1 und 2 (sie-

he Tabelle 4.1). Wir wollen herauszufinden, ob ein Ubergang zur NB-Verteilung die Anpassung

verbessert. Bezieht man NB-Modelle mit in den Modellvergleich ein, hangt der Bayes-Faktor

neben der Wahl von φ insbesondere auch vom im NB-Modell gewahlten Dispersionsparame-

ter ab. Die Funktion BFnb2ln (siehe Anhang, Seite 94ff) berucksichtigt dies. Im Vergleich zur

Funktion BFpois2ln, die eine Gegenuberstellung zweier Poissonmodelle ermoglicht, sind dieser

Splus-Funktion zusatzlich die Dispersionsparameter der beiden Modelle und die Lange des Da-

tensatzes zu ubergeben.


geschaetzter Erwartungswert der Patente

ge

sch

ae

tzte

Va

ria

nz

de

r P

ate

nte

0 50 100 150 200

01

00

02

00

03

00

0

Abbildung 4.1: Uberdispersion im Poissonmodell 2 der Patent-Daten (— Winkelhalbierende,

– – – geglatteter Zusammenhang)

Mit folgender Splus-Eingabe lasst sich beispielsweise der skalierte Bayes-Faktor fur den Ver-

gleich des NB-Modells 3 mit Dispersionsparameter a3 = 0.3 mit dem Poissonmodell 1 (a1 = 0)

fur φ = 1 berechnen.

patente31<-BFnb2ln(0.3,0,1,patente∼offset(log(rep(1,70)))+rdsales+rdw5,patente∼offset(log(rep(1,70)))+rdsales+rdw5,patente,log(rep(1,70)),70)

In den Grafiken in Abbildung 4.2 sind die logarithmierten Bayes-Faktoren (Log-Bayes-Faktoren)

fur die vier Vergleiche zwischen den NB-Modellen 3 und 4 und den Poissonmodellen 1 und 2 in

Abhangigkeit vom im NB-Modell gewahlten Dispersionparameter und jeweils fur φ = (1, 1.65, 5)

aufgetragen. Fur die Beschriftung der Grafiken wurde folgende Bezeichnung gewahlt:

logBij(ai, aj, phi) := logBij(ai, aj , φ), i = 3, 4; j = 1, 2

Ein logarithmierter Bayes-Faktor großer als Null spricht somit fur das jeweils gewahlte NB-

Modell. Betrachten wir zunachst nur die Grafik links oben zum Vergleich des NB-Modells 3


a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB31(a3,0,phi)

a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB31(a3,0,phi)

a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

phi=1phi=1.65phi=5

a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB32(a3,0,phi)

a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB32(a3,0,phi)

a3

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

phi=1phi=1.65phi=5

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB41(a4,0,phi)

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

logB41(a4,0,phi)

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

phi=1phi=1.65phi=5

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0logB42(a4,0,phi)

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0logB42(a4,0,phi)

a4

log

Bay

es fa

ctor

0 5 10 15 20 25 30

020

4060

8010

0

phi=1phi=1.65phi=5

Abbildung 4.2: Log-Bayes-Faktoren als Funktion des zum jeweiligen NB-Modell gehorigen Uber-

dispersionsparameters fur die NB-Poisson-Modellvergleiche der Patent-Daten in Abhangigkeit

vom gewahlten Hyperparameter φ

mit dem Poissonmodell 1. Fur einen sehr kleinen Dispersionsparameter (a3 → 0) geht das NB-

Modell 3 in das Poissonmodell 1 uber, da beide Modelle die gleichen Kovariablen enthalten

und die NB-Verteilung die Poissonverteilung als Spezialfall enthalt (siehe Kapitel 2, Seite 32).

Deshalb startet die Kurve mit einem Log-Bayes-Faktor von nahezu Null. Bei steigendem Dis-

persionsparameter a3 erhoht sich die Evidenz fur das NB-Modell 3 sehr schnell. Ubersteigt

der Dispersionsparameter allerdings einen gewissen Optimalwert, sinkt der Log-Bayes-Faktor


wieder. Die Evidenz fur das NB-Modell 3 wird schwacher. Die Kurven zeigen deutlich, dass ein

Ubergang zur NB-Verteilung die Anpassung verbessert. Es bestatigt sich unsere Vermutung,

dass die Patent-Daten Uberdispersion aufweisen. Die Kurven deuten gleichzeitig an, in welchem

Bereich der Dispersionsparameter zu wahlen ist, um eine moglichst gute Anpassung des NB-

Modells 3 an die Daten zu erreichen. Diese Interpretation trifft auch auf die Grafik rechts unten

fur den Vergleich des NB-Modells 4 mit dem Poissonmodell 2 zu.

Betrachtet man die Grafiken rechts oben und links unten genauer, so stellt man fest, dass

diese Kurven nicht bei Null beginnen. Das liegt daran, dass mit diesen Kurven Modelle ver-

glichen werden, die nicht dieselben Kovariablen aufweisen. Die Grafik links unten dient dem

Vergleich von NB-Modell 4 mit dem Poissonmodell 1. Weil das NB-Modell 4 im Gegensatz zum

Poissonmodell 1 zusatzlich zu den Kovariablen rdsales und rdw5 auch noch die Interaktion

zwischen diesen Kovariablen enthalt, ist es schon unabhangig von der NB-Verteilungsannahme

besser als das Poissonmodell 1. Lassen wir den Dispersionsparameter a4 nahezu Null werden,

so geht der Vergleich zwischen dem NB-Modell 4 und dem Poissonmodell 1 namlich in den Ver-

gleich zwischen Poissonmodell 2 und Poissonmodell 1 uber. Wie bereits gezeigt wurde, liefert das

Poissonmodell 2 eine wesentlich bessere Anpassung an die Daten. Fur den Log-Bayes-Faktor zu

diesem Vergleich erhalt man entsprechend den Werten in (4.1) logB21(0, 0, φ) = (3.72, 2.88, 2.85).

Deshalb beginnt die Kurve fur sehr kleines a4 nicht bei Null, sondern bei einem Wert von ca. 3.

Analog dazu beginnt auch die Kurve fur den Vergleich zwischen dem NB-Modells 3 und

dem Poissonmodell 2 nicht bei Null. Bei diesem Modellvergleich enthalt allerdings das Poisson-

modell 2 die Interaktion zwischen rdsales und rdw5 als zusatzliche Kovariable. Wahlen wir

den Dispersionsparameter a3 des NB-Modells 3 sehr klein (a3 → 0) so geht der Vergleich zwi-

schen dem NB-Modell 3 und dem Poissonmodell 2 in den Vergleich zwischen dem Poissonmodell

1 und dem Poissonmodell 2 uber. Der Log-Bayes-Faktor fur diesen Modellvergleich liegt bei

logB12(0, 0, φ) = (−3.72,−2.88,−2.85). Dementsprechend beginnt die Kurve fur sehr kleines a3

nicht bei Null, sondern bei einem Wert von ca. -3.

Wie sich insbesondere in den oberen beiden Grafiken der Abbildung 4.2 zeigt, hat die Wahl

des Hyperparameters φ nur sehr geringfugig Einfluss auf die Hohe der Log-Bayes-Faktoren. Erst

wenn man den Dispersionsparameter sehr groß wahlt, weichen die zu den drei verschiedenen

φ-Werten korrespondierenden Kurven etwas voneinander ab.

Fur die Berechnung der Log-Bayes-Faktoren benotigen wir die Log-Likelihoods der zu verglei-

chenden Modelle. Diese und die korrespondierenden ML-Schatzer fur die Dispersionsparameter

(aML) findet man in Tabelle 4.4. Die Log-Likelihoods lassen sich mit der Funktion BFnb2ln

berechnen, die ML-Schatzer fur die Dispersionsparameter der NB-Modelle 3 und 4 wurden mit

der Splus-Funktion glm.nb bestimmt. Dabei ist zu beachten, dass diese Funktion den geschatz-

ten Dispersionsparameter unter dem Namen Theta ausgibt, wobei Theta= 1aML gilt. Der nur fur

genestete Modelle definierte Likelihood-Quotienten-Test (siehe Kapitel 2, Seite 21) liefert fur die


Modell- M1 M0

vergleich

M1 −M0 Log-Likelihood aML1 Log-Likelihood aML

0

3 - 1 5115.45 0.31 5007.65 0

4 - 1 5116.81 0.29 5007.65 0

3 - 2 5115.45 0.31 5012.59 0

4 - 2 5116.81 0.29 5012.59 0

Tabelle 4.4: Log-Likelihoods und zugehorige ML-Schatzer der Dispersionsparameter fur die NB-

Poisson-Modellvergleiche der Patent-Daten

Vergleiche zwischen NB-Modell 3 und Poissonmodell 1, zwischen NB-Modell 4 und

Poissonmodell 1 und zwischen NB-Modell 4 und Poissonmodell 2 jeweils einen p-Wert von Null.

Die Hypothese βr = 0 lasst sich somit bei jedem dieser drei Modellvergleiche ablehnen. Dabei

bezeichnet βr die Parameter, die im NB-Modell zusatzlich zu denen im Poissonmodell enthalten

sind. Das jeweilige NB-Modell passt die Daten wesentlich besser an.

Beim Vergleich von NB-Modell 3 mit Poisson-Modell 2 handelt es sich um nicht-genestete

Modelle. Dementsprechend kann der Likelihood-Quotienten-Test nicht benutzt werden, um die

Anpassungsgute von NB-Modell 3 und Poissonmodell 2 zu vergleichen. An dieser Stelle stehen

uns die Bayes-Faktoren zur Verfugung. Um einerseits die Aussagen des Likelihood-Quotienten-

Tests fur den Vergleich zwischen NB-Modell 3 und Poissonmodell 1, zwischen NB-Modell 4 und

Poissonmodell 1 und zwischen NB-Modell 4 und Poissonmodell 2 zu bestatigen und andererseits

auch eine Evidenzaussage fur den Vergleich zwischen NB-Modell 3 und Poissonmodell 2 treffen

zu konnen, wurden in Tabelle 4.5 die zu den Grafiken in Abbildung 4.2 korrespondierenden

maximalen Log-Bayes-Faktoren logBmax10 mit den zugehorigen Dispersionsparametern aBF zu-

sammengefasst.

Modell- φ = 1 φ = 1.65 φ = 5

vergleich logBmax10 aBF1 aBF0 logBmax

10 aBF1 aBF0 logBmax10 aBF1 aBF0

3 - 1 110.54 0.33 0 110.52 0.33 0 110.51 0.33 0

4 - 1 114.67 0.31 0 111.59 0.32 0 110.47 0.32 0

3 - 2 106.83 0.33 0 107.65 0.33 0 107.66 0.33 0

4 - 2 110.95 0.31 0 108.72 0.32 0 107.62 0.32 0

Tabelle 4.5: Maximale Log-Bayes-Faktoren und deren Uberdispersionsparameter fur die NB-

Poisson-Modellvergleiche der Patent-Daten in Abhangigkeit vom gewahlten Hyperparameter φ

Bei allen vier Modellvergleichen spricht sich der maximale Log-Bayes-Faktor entscheidend fur das


jeweils gewahlte NB-Modell aus. Wie zuvor beim Vergleich der Poissonmodelle 1 und 2 zeigt sich,

dass die Hinzunahme der Interaktion die Evidenz fur das jeweilige NB-Modell nochmals erhoht,

allerdings nicht in signifikantem Ausmaße. Betrachten wir beispielsweise den Vergleich zwischen

NB-Modell 3 und Poissonmodell 1, so liegt der maximale Log-Bayes-Faktor fur φ = 1 bei 110.54.

Dieser Wert spricht entscheidend fur das NB-Modell 3. Fur den Vergleich des NB-Modells 4 mit

dem Poissonmodell 1 erhalten wir einen noch etwas hoheren maximalen Log-Bayes-Faktor. Die

Hinzunahme der Interaktion lasst den maximalen Log-Bayes-Faktor auf 114.67 ansteigen. Ana-

loges beobachten wir beim Vergleich der NB-Modelle 3 und 4 mit dem Poissonmodell 2. Fur

φ = 1.65 und φ = 5 bewirkt die Hinzunahme der Interaktion keine wesentlich Erhohung bzw.

sogar eine leichte Senkung des maximalen Log-Bayes-Faktors. Wie bereits in Kapitel 3 (Seite 67)

erwahnt, verlagert sich die Evidenz mit steigendem φ immer mehr in Richtung des kleineren Mo-

dells. Aus diesem Grund sinkt der maximale Log-Bayes-Faktor fur den Vergleich von NB-Modell

4 mit dem Poissonmodell 1 mit steigendem φ. Fur den Vergleich zwischen dem NB-Modell 3 und

dem Poissonmodell 2 geht der maximalen Log-Bayes-Faktor fur steigendes φ nach oben, weil in

diesem Fall das Poissonmodell 2 eine Kovariable mehr besitzt. Die Signifikanzaussage wird aber

bei keinem der in Tabelle 4.5 aufgefuhrten Modellvergleiche durch die Wahl von φ beeinflusst.

Die NB-Modelle 3 und 4 liefern eine wesentlich bessere Anpassung als die Poissonmodelle 1 und

2. Da die Daten demnach starke Uberdispersion aufweisen, gilt es nun, das NB-Modell 4 dem

NB-Modell 3 gegenuberzustellen.

Die Log-Bayes-Faktoren fur den direkten Vergleich zwischen NB-Modell 4 und NB-Modell 3

lassen sich aus den bereits bekannten maximalen Log-Bayes-Faktoren in Tabelle 4.5 berechnen.

Aufgrund der Definition des Bayes-Faktors (siehe Gleichung (3.5)) gilt:

logBij = log(P (D|Mi)P (D|Mj)

)= log

⎛⎝ P (D|Mi)

P (D|Mv)

P (D|Mj)P (D|Mv)

⎞⎠

= log(P (D|Mi)P (D|Mv)

)− log

(P (D|Mj)P (D|Mv)

)

= log(Biv) − log(Bjv) (4.4)

Dabei bezeichne Mv ein beliebiges Vergleichsmodell. Wahlen wir das Poissonmodell 1 als Ver-

gleichsmodell, so konnen wir nun mit Hilfe der maximalen Log-Bayes-Faktoren in Tabelle 4.5

die Log-Bayes-Faktoren fur den Vergleich zwischen NB-Modell 4 mit Dispersionsparameter aBF4

und NB-Modell 3 mit Dispersionsparameter aBF3 berechnen:

logB43(aBF4 , aBF3 , φ) = logB41(aBF4 , 0, φ) − logB31(aBF3 , 0, φ)

= (logB41(0.31, 0, 1) − logB31(0.33, 0, 1), log B41(0.32, 0, 1.65)

− logB31(0.33, 0, 1.65), log B41(0.32, 0, 5) − logB31(0.33, 0, 5))


= (114.67 − 110.54, 111.59 − 110.52, 110.47 − 110.51)

= (4.13, 1.07,−0.04) (4.5)

Die entsprechenden skalierten Bayes-Faktoren fur den Vergleich zwischen NB-Modell 4 und NB-

Modell 3 lauten

2 logB43(aBF4 , aBF3 , φ) = (8.26, 2.14,−0.08). (4.6)

Mit steigendem φ sinkt die Evidenz fur das großere NB-Modell 4. Fur φ = 1 spricht der skalierte

Bayes-Faktor nach Tabelle 3.1 stark fur das NB-Modell 4, fur φ = 1.65 ist die Evidenz fur das

NB-Modell 4 nicht mehr wert als bloße Erwahnung, und fur φ = 5 spricht der skalierte Bayes-

Faktor sogar fur das NB-Modell 3. Die Wahl von φ beeinflusst die Evidenzaussage entscheidend.

Wenn wir die maximalen Log-Bayes-Faktoren zur Beurteilung der Evidenz benutzen, ist al-

lerdings zu beachten, dass wir die Unsicherheit in der Parameterschatzung missachten. Wird

der Dispersionsparameter nicht im voraus gewahlt, sondern stattdessen geschatzt, tendieren die

maximalen Log-Bayes-Faktoren in Tabelle 4.5 dazu, die Evidenz fur das jeweilige NB-Modell

zu uberschatzen. Der Dispersionsparameter ist nicht - wie bei der Berechnung der maximalen

Log-Bayes-Faktoren vorausgesetzt - bekannt, sondern unterliegt einer Verteilung. Deshalb ist es

besser, die in Abschnitt 3.5 eingefuhrten Durchschnitts-Bayes-Faktoren fur die Modellvergleiche

zu benutzen. Sie berucksichtigen, dass der Dispersionsparameter unbekannt ist und beziehen

somit die Unsicherheit in der Schatzung des Dispersionsparameters mit ein. Der Dispersionspa-

rameter wird in gleicher Weise wie die Regressionsparameter βj (j=1,. . . ,p) als ein Parameter

betrachtet. Dazu gibt man ihm eine a priori Verteilung und integriert aus. Die durchschnittlichen

Log-Bayes-Faktoren fur den Vergleich der NB-Modelle 3 und 4 mit den Poissonmodellen 1 und

2 berechnen wir demnach folgendermaßen:

log Bij(φ) :=∫

logBij(ai, 0, φ)π(ai)dai, i = 3, 4; j = 1, 2 (4.7)

Dabei bezeichnet π(ai) die a priori Dichte des zum NB-Modell Mi gehorigen Dispersionsparame-

ters ai. Entsprechend der Kurvenverlaufe in der Abbildung 4.2 integrieren wir uber ai ∈ (0, 30].

Da die Gammaverteilung analog zum Dispersionsparameter in einem NB-Modell nur fur positive

Werte definiert ist, nehmen wir a priori an, dass der Dispersionsparameter ai Γ(n, λ)-verteilt ist

mit Erwartungswert nλ = 0.25 (alternativ 0.1) und Varianz n

λ2 = 0.5 (n, λ > 0). Damit gilt fur

die a priori Dichte π(ai), i=3,4:

π(ai) =λn

Γ(n)an−1i e−λai1(0,∞)(ai) mit n = 0.125, λ = 0.5 bzw. n = 0.02, λ = 0.2 (4.8)

Die errechneten durchschnittlichen Log-Bayes-Faktoren wurden in Tabelle 4.6 zusammenge-

fasst. Vergleicht man die durchschnittlichen Log-Bayes-Faktoren mit den maximalen Log-Bayes-


Modell- log B10(φ) mit log B10(φ) mit

vergleich a1 ∼ Γ(0.125, 0.5) a1 ∼ Γ(0.02, 0.2)

(a priori mehr Uberdispersion) (a priori weniger Uberdispersion)

M1 - M0 φ = 1 φ = 1.65 φ = 5 φ = 1 φ = 1.65 φ = 5

3 - 1 44.07 44.06 44.05 10.18 10.18 10.17

4 - 1 47.73 45.76 45.26 11.29 10.76 10.63

3 - 2 41.00 41.68 41.70 9.22 9.44 9.44

4 - 2 44.67 43.38 42.90 10.33 10.02 9.90

Tabelle 4.6: Durchschnittliche Log-Bayes-Faktoren fur die NB-Poisson-Modellvergleiche der

Patent-Daten in Abhangigkeit von φ und der gewahlten a priori Verteilung fur den Disper-

sionsparameter im jeweiligen NB-Modell

Faktoren in Tabelle 4.5, so sieht man sofort, dass die maximalen Log-Bayes-Faktoren die Evidenz

fur das jeweilige NB-Modell wesentlich uberschatzen. Doch auch die um einiges kleineren durch-

schnittlichen Log-Bayes-Faktoren sprechen eindeutig fur die NB-Modelle. Die Wahl von φ hat

keinen Einfluss auf die Evidenzaussagen. Unter der Annahme, dass der Dispersionsparameter ai(i=3,4) priori Gamma-verteilt ist mit Erwartungswert 0.25 und Varianz 0.5 liegen die Werte der

Durchschnitts-Bayes-Faktoren zwischen 40 und 50. Senkt man den Erwartungswert des

Dispersionsparameters auf 0.1, so fallen die durchschnittlichen Log-Bayes-Faktoren auf einen

Wert von ungefahr 10. Nimmt man a priori an, dass der Erwartungswert des Dispersionspara-

meters bei 0.25 liegt, so erhalt eventuell in den Daten vorhandene Uberdispersion relativ viel

Gewicht. Senkt man den a priori Erwartungswert auf 0.1, so gewichtet man mehr, dass kei-

ne bzw. vergleichsmaßig weniger Uberdispersion in den Daten vorliegt. Zur Veranschaulichung

betrachte man die Abbildung 4.3, in der die Gammaverteilungsdichten mit Erwartungswert

0.25 bzw. Erwartungswert 0.1 und Varianz 0.5 zu sehen sind. Wie sich den bisher berechneten

Bayes-Faktoren entnehmen lasst, weisen die von uns betrachteten Patent-Daten definitiv Uber-

dispersion auf. Deshalb sprechen die durchschnittlichen Log-Bayes-Faktoren auch dann noch

eindeutig fur das jeweilige NB-Modell, wenn man den Erwartungswert in der a priori Dichte

relativ klein wahlt und somit mehr Gewicht auf wenig Uberdispersion legt.

Aus dem Vergleich der Poissonmodelle 1 und 2 wissen wir bereits, dass das Poissonmodell 2 - das

mit der Interaktion - dem Poissonmodell 1 vorzuziehen ist. Der Tabelle 4.6 ist zu entnehmen,

dass die NB-Modelle 3 und 4 im Vergleich zu den Poissonmodellen 1 und 2 immer eine viel

bessere Anpassung liefern.

Die durchschnittlichen Log-Bayes-Faktoren fur den direkten Vergleich der NB-Modelle 4 und

3 lassen sich entsprechend Gleichung (4.4) aus den durchschnittlichen Log-Bayes-Faktoren in

Tabelle 4.6 berechnen. Da man beim Vergleich zweier NB-Modelle a priori von Uberdispersion


Gamma(0.125,0.5), EW=0.25

x

Dic

hte

0 1 2 3 4 5

0.0

0.0

20.0

40.0

60.0

80.1

0

Gamma(0.02,0.2), EW=0.1

x

Dic

hte

0 1 2 3 4 5

0.0

0.0

20.0

40.0

60.0

80.1

0

Abbildung 4.3: Gammaverteilungsdichte mit Erwartungswert 0.25 und 0.1 (Varianz = 0.5)

in den Daten ausgeht, seien die Dispersionsparameter a3 und a4 a priori Gamma-verteilt mit

Erwartungswert 0.25 und Varianz 0.5. Wir wahlen wieder das Poissonmodell 1 als Vergleichs-

modell. Dann gilt:

log B43(φ) =∫ ∫

logB43(a4, a3, φ)π(a4)π(a3)da4da3

(4.4)=

∫ ∫logB41(a4, 0, φ)π(a4)π(a3)da4da3

−∫ ∫

logB31(a3, 0, φ)π(a4)π(a3)da4da3

=∫

logB41(a4, 0, φ)π(a4)[∫

π(a3)da3

]︸︷︷︸

=1, da π(a3)Dichte

da4

−∫

logB31(a3, 0, φ)π(a4)[∫

π(a4)da4

]︸︷︷︸

=1, da π(a4)Dichte

da3


=∫

logB41(a4, 0, φ)π(a4)da4 −∫

logB31(a3, 0, φ)π(a3)da3 =

(4.7)= log B41(φ) − log B31(φ)

Tabelle4.6= (47.73 − 44.07, 45.76 − 44.06, 45.26 − 44.05)

= (3.66, 1.70, 1.21) (4.9)

Wie erwartet sinkt der durchschnittliche Log-Bayes-Faktor fur steigendes φ, da das NB-Modell

4 eine Kovariable mehr enthalt. Als Werte fur den skalierten durchschnittlichen Bayes-Faktor

erhalten wir

2 log B43(φ) = (7.32, 3.40, 2.42), φ = (1, 1.65, 5). (4.10)

Fur φ = 1 passt das NB-Modell 4 die Daten besser an als das NB-Modell 3. Der skalierte Bayes-

Faktor spricht stark fur das NB-Modell 4. Fur φ = 1.65 und φ = 5 ist die Evidenz fur das

NB-Modell 4 immer noch positiv. Insgesamt liefert also das NB-Modell 4 die beste Anpassung.

Um das Verhaltnis von Devianz und Freiheitsgraden bei diesem Modell zu uberprufen, lassen

wir in Splus folgendes Modell berechnen:

patente4<-glm.nb(patente∼rdsales+rdw5+rdsales*rdw5,link=log,x=T) (4.11)

Das zugehorige Splus-summary findet man in Tabelle 4.7.

Coefficients:

Value Std. Error t value P(>|t|)

(Intercept) -1.86 0.398 -4.66 3.15e-006

rdsales -8.23 5.998 -1.37 0.17

rdw5 2.64 0.251 10.48 0

rdsales:rdw5 5.74 3.963 1.45 0.15

(Dispersion Parameter for Negative Binomial family taken to be 1 )


Residual Deviance: 78.5 on 66 degrees of freedom


Theta: 3.40016

Std. Err.: 0.9251

2 x log-likelihood: 10233.61782

Tabelle 4.7: Splus-summary des NB-Modells 4 fur die Patent-Daten


Das Verhaltnis der Devianz zu den Freiheitsgraden hat sich im Vergleich zum Poissonmodell

2 (siehe Tabelle 4.3) mit 78.5 zu 66 immens verbessert. Die starke Senkung der Devianz von

377 im Poissonmodell 2 auf jetzt 78.5 weist darauf hin, dass das NB-Modell eine wesentlich

bessere Anpassung als das Poissonmodell 2 liefert. Die dazu bereits in [Siko 02] (Seite 140f)

durchgefuhrte Residuenanalyse spiegelt die verbesserte Anpassung wider.

Es bleibt zu erwahnen, dass der in Tabelle 4.7 angegebene, zum t-Test korrespondierende

p-Wert fur die Interaktion trotz seines mit 0.15 relativ hohen Wertes nicht dahingehend gedeu-

tet werden kann, dass die Interaktion nicht signifikant ist und dementsprechend aus dem Modell

entfernt werden sollte. Da das NB-Modell 3 und das NB-Modell 4 nicht genestet sind, darf der

t-Test nicht zum Vergleich dieser Modelle herangezogen werden. Die in Tabelle 4.7 aufgefuhrten

Teststatistiken (t-value) vergleichen das NB-Modell 4 unter Festhaltung des Dispersionspara-

meters mit dem NB-Modell ohne die jeweils zu testende Kovariable. Das zur Uberprufung der

Signifikanz der Interaktion herangezogene Vergleichsmodell ist somit nicht mit dem von uns be-

trachteten NB-Modell 3 identisch.

Der zweite Anwendungsteil behandelt die Anpassung verschiedener Regressions-

modelle an einen Kfz-Haftpflicht-Datensatz und den anschließenden Vergleich der

angepassten Modelle mit Hilfe der Bayes-Faktoren. Da der betrachtete Datensatz

von der Versicherungskammer Bayern zur Verfugung gestellt wurde, kann diese

Anwendung aus datenschutzrechtlichen Grunden nicht veroffentlicht werden.

Kapitel 5

Zusammenfassung und Ausblick

Zum Abschluss dieser Arbeit soll noch einmal hervorgehoben werden, dass die Bayes-Faktoren

eine entscheidende Erweiterung der bisher im Rahmen der herkommlichen Signifikanztests mogli-

chen Modellregression darstellen. Wie bereits in der Einleitung erwahnt, sind speziell in der

Kfz-Tarifierung zur Modellierung von Schadenanzahlen oft Vergleiche nicht-genesteter Regressi-

onsmodelle notwendig. Dies beruht zum einen darauf, dass neu entwickelte Merkmale bestimm-

ten Merkmalen, die bereits in die Tarifentwicklung einfließen, gegenubergestellt werden mussen.

Die zu vergleichenden Modelle sind dann nicht genestet, da sie unterschiedliche Kovariablen

enthalten. Zum anderen weisen Zahldaten oft Uberdispersion auf. Wird versucht, die auf Uber-

dispersion beruhende mangelnde Anpassungsgute von Poissonmodellen durch den Ubergang

auf NB-Modelle zu beheben, hat man in einem nachsten Schritt NB-Modelle hinsichtlich ihrer

Anpassungsgute zu vergleichen. Diese sind aufgrund des zusatzlich zu schatzenden Dispersions-

parameters unabhangig von den enthaltenen Kovariablen nicht genestet. Die herkommlichen, auf

p-Werten basierenden Signifikanztest, die standardmaßig zum Vergleich verschiedener Regressi-

onsmodelle benutzt werden, sind nur fur genestete Modelle definiert. Die in dieser Diplomarbeit

eingefuhrten Bayes-Faktoren hingegen stellen diese Anforderung nicht. Mit ihrer Hilfe lassen sich

auch nicht-genestete Regressionsmodelle bzgl. ihrer Anpassungsgute vergleichen.

In der statistischen Praxis verwendet man auf p-Werten basierende Ruckwarts- oder Vorwarts-

analysen, um die Signifikanz einzelner Kovariablen in einem Poissonmodell zu testen. Diese Ver-

fahren zur Modellanpassung sind aus den bereits genannten Grunden nicht geeignet, wenn man

ein NB-Modell untersucht. Eine Software-Routine zur Anpassung von NB-Modellen, die die not-

wendigen Modellvergleiche mit Hilfe der Bayes-Faktoren durchfuhrt, wurde es ermoglichen, auch

die Signifikanz einzelner Merkmale in einem NB-Modell zu beurteilen.

Im Anwendungsteil dieser Diplomarbeit ist der Fall aufgetreten, dass die Bayes-Faktoren

die im Rahmen einer Ruckwartsanalyse bestatigte Signifikanz einzelner Kovariablen in einem

Poissonmodell widerlegen. Zu einem bereits angepassten Modell wurden Kovariablen hinzuge-

nommen, deren p-Werte das geforderte Signifikanzniveau unterschreiten. Die Bayes-Faktoren

88

89

sagen allerdings aus, dass die Hinzunahme der Kovariablen die Modellanpassung verschlechtert.

Dies legt die Vermutung nahe, dass die zusatzlich in das Modell aufgenommenen Kovariablen

zu einer Uberanpassung des Modells fuhren (Overfitting). Die Bayes-Faktoren scheinen darauf

sensibler zu reagieren als die herkommlichen auf p-Werten basierenden Signifikanztests. Eine

umfangreiche Simulationsstudie konnte Aufschluss daruber geben, inwieweit diese Vermutung

gerechtfertigt ist.

Anhang

1. Funktionen zur Berechnung des skalierten Bayes-Faktors fur den Vergleich von

Modell M1 mit Modell M0 mittels Approximation (3.23) unter der Annahme, dass

die Regressionsparameter a priori Normal-verteilt sind.

Mit folgender Splus-Funktion lasst sich der approximierte skalierte Bayes-Faktor fur den Ver-

gleich zweier Poissonmodelle berechnen. Dieser Funktion sind der Hyperparameter φ, die

Modelle, der Response und der Offset zu ubergeben.

[1] BFpois2ln<-function(phi,model1,model0,response,offs)

[2]

[3] psi=1 #fest gewaehlte a priori Parameter

[4] nu1=0

[5] Modell<-glm(model1,family=poisson,link=log,x=T) #PoisGLM M1

[6] coef<-matrix(Modell$coef,ncol=1) #Parameterschaetzer

[7] eta<-Modell$x%*%coef #linearer Praediktor

[8] mu<- care.exp(offs + eta) #Link

[9] deta<- 1/mu #deta=d(eta)/d(mu)

[10] v <- mu #Varianzfunktion

[11] w<-1/(deta^2 * v) #Gewichtung

[12] w<-w/sum(w)

[13] z<-eta + (response - mu) *deta #adjusted dependent variable z

[14] xbar<-t((Modell$x)[1:length(response),-1]) %*%w

[15] #gewichtete Stichprobenerwartungswerte der Kovariablen

[16] xbar2<-t(((Modell$x)[1:length(response),-1])^2)%*%w

[17] s2x<-xbar2-xbar^2 #gewichtete Stichprobenvarianzen der Kovariablen

[18] zbar<-as.numeric(weighted.mean(z,w)) #gewichteter

[19] #Stichprobenerwartungswert der adjusted dependent variable z

[20] zbar2<-t(z^2)%*%w

90

Anhang 91

[21] s2z <- zbar2 - zbar^2 #gewichtete Stichprobenvarianz der adjusted

[22] #dependend variable z

[23] varz<- as.numeric(s2z)

[24] m<-ncol(as.matrix((Modell$x)[1:length(response),-1]))

[25] pmean<-c(nu1 + zbar,rep(0,times=m)) #a priori Erwartungswert

[26] #Erzeugung der Matrix U

[27] U<-diag(c(psi^2,rep(phi^2,times=m)),nrow=m+1)

[28] #Erzeugung der Matrix Q

[29] Q<-diag(1/sqrt(s2x),nrow=length(s2x))

[30] Q<-rbind(t(-xbar/sqrt(s2x)),Q)

[31] Qcol1<-c(1,rep(0,times=m))

[32] Q<-sqrt(varz) * cbind(Qcol1,Q)

[33] pvar<- Q %*% U %*% t(Q) #a priori Kovarianzmatrix

[34] CovMatrix<-vcov.glm(Modell) #Kovarianzmatrix der Parameterschaetzer

[35] V<-CovMatrix

[36] A<-solve(V) #asymptotische Fisher-Informationsmatrix

[37] thetahat<-matrix(Modell$coef,ncol=1) #Vektor der Parameterschaetzer

[38] G <- solve(pvar) #Inverse der a priori Kovarianzmatrix

[39] H<-solve(A + G)

[40] I<-diag(rep(1,times=m+1),nrow=m+1) #Einheitsmatrix (Rang=m+1)

[41] logdetW<-sum(log(eigen(pvar)$values))

[42] C<- G %*% (I - H %*% (2*I - A %*% H) %*% G)

[43] E1<- -(t(thetahat - pmean))%*% C %*% (thetahat - pmean) - logdetW

[44] - sum(log(eigen(A+G)$values))

[45]

[46] #Vergleichsmodell: Poissonmodell M0

[47]

[48] Modell0<-glm(model0,family=poisson,link=log,x=T) #PoisGLM M0

[49] coef0<-matrix(Modell0$coef,ncol=1) #Parameterschaetzer

[50] eta0<-Modell0$x%*%coef0 #linearer Praediktor

[51] mu0<- care.exp(offs + eta0) #Link

[52] deta0<- 1/mu0 #deta=d(eta)/d(mu)

[53] v0 <- mu0 #Varianzfunktion

[54] w0<-1/(deta0^2 * v0) #Gewichtung

[55] w0<-w0/sum(w0)

[56] z0<-eta0 + (response - mu0) *deta #adjusted dependent variable z

[57] xbar0<-t((Modell0$x)[1:length(response),-1]) %*%w0

92 Anhang


[59] xbar20<-t(((Modell0$x)[1:length(response),-1])^2)%*%w0

[60] s2x0<-xbar20-xbar0^2 #gewichtete Stichprobenvarianzen der Kovariablen

[61] zbar0<-as.numeric(weighted.mean(z0,w0)) #gewichteter


[63] zbar20<-t(z0^2)%*%w0

[64] s2z0 <- zbar20 - zbar0^2 #gewichtete Stichprobenvarianz der adjusted


[66] varz0<- as.numeric(s2z0)

[67] m0<-ncol(as.matrix((Modell0$x)[1:length(response),-1]))

[68] pmean0<-c(nu1 + zbar0,rep(0,times=m0)) #a priori Erwartungswert

[69] #Erzeugung der Matrix U0

[70] U0<-diag(c(psi^2,rep(phi^2,times=m0)),nrow=m0+1)

[71] #Erzeugung der Matrix Q0

[72] Q0<-diag(1/sqrt(s2x0),nrow=length(s2x0))

[73] Q0<-rbind(t(-xbar0/sqrt(s2x0)),Q0)

[74] Qcol10<-c(1,rep(0,times=m0))

[75] Q0<-sqrt(varz0) * cbind(Qcol10,Q0)

[76] pvar0<- Q0 %*% U0 %*% t(Q0) #a priori Kovarianzmatrix

[77] CovMatrix0<-vcov.glm(Modell0) #Kovarianzmatrix der Parameterschaetzer

[78] V0<-CovMatrix0

[79] A0<-solve(V0) #asymptotische Fisher-Informationsmatrix

[80] thetahat0<-matrix(Modell0$coef,ncol=1) #Vektor der Parameterschaetzer

[81] G0 <- solve(pvar0) #Inverse der a priori Kovarianzmatrix

[82] H0<-solve(A0 + G0)

[83] I0<-diag(rep(1,times=m0+1),nrow=m0+1) #Einheitsmatrix (Rang=m0+1)

[84] logdetW0<-sum(log(eigen(pvar0)$values))

[85] C0<- G0 %*% (I0 - H0 %*% (2*I0 - A0 %*% H0) %*% G0)

[86] E0<- -(t(thetahat0 - pmean0))%*% C0 %*% (thetahat0 - pmean0) - logdetW0

[87] - sum(log(eigen(A0+G0)$values))

[88] LogLikeM1<- sum( response * (offs+eta) - mu) #Log-Likelihood Modell 1

[89] LogLikeM0<- sum( response * (offs+eta0) - mu0) #Log-Likelihood Modell 0

[90] E<-E1-E0

[91] BF2ln<- 2*(LogLikeM1 - LogLikeM0) + (E1 - E0) #skalierter Bayes-Faktor

[92] list(BF2ln=BF2ln,model1=model1,model0=model0,

[93] LogLikeM1=LogLikeM1,LogLikeM0=LogLikeM0,E=E) #Ausgabe

[94]

Anhang 93

Will man den Vergleich zum Nullmodell, d.h. zum Modell nur mit Intercept ziehen, muss man

in der Funktion BFpois2ln die Zeilen 46 bis einschließlich 87 durch folgende Programmzeilen

ersetzen.

[46] #Poisson-NULLMODELL (nur mit Intercept):

[47]




[51] mu0<-care.exp(offs + eta0) #Linkfunktion

[52] deta0<-1/mu0 #deta=d(eta)/d(mu)

[53] v0<-mu0 #Varianzfunktion

[54] w0<-1/(deta0*v0) #Gewichtung

[55] w0<-w0/sum(w0)

[56] z0<-eta0 + (response - mu0) *deta0 #adjusted dependent variable z


[58] #Stichprobenerwartungswert der adjusted dependent variable

[59] zbar20<-t(z0^2)%*%w0




[63] pmean0<- nu1 + zbar0 #a priori Erwartungswert

[64] CovMatrix0<-vcov.glm(Modell0) #Kovarianzmatrix der Parameterschaetzer

[65] V0<-CovMatrix0

[66] A0<-1/V0 #asymptotische Fisher-Informations"matrix"



[69] Qcol10<-1

[70] Q0<-sqrt(varz0) *Qcol10


[72] U0<-psi^2

[73] pvar0<- Q0 * U0 * Q0 #a priori Kovarianzmatrix

[74] G0<- 1/pvar0 #Inverse der a priori Kovarianzmatrix

[75] H0<-1/(A0 + G0)

[76] C0<- G0 * (1 - H0 * (2 - A0 * H0) * G0)

[77] E0<- -(thetahat0 - pmean0)* C0 * (thetahat0 - pmean0) - log(pvar0)

[78] - log(A0 + G0)

94 Anhang

Mit folgender Splus-Funktion lasst sich der approximierte skalierte Bayes-Faktor fur den Ver-

gleich zweier NB-Modelle berechnen. Dieser Funktion sind die Dispersionsparameter a1 und

a2 der NB-Modelle, der Hyperparameter φ, die NB-Modelle, der Response, der Offset und die

Lange des Datensatzes zu ubergeben.

[1] BFnb2ln<-function(a1,a0,phi,model1,model0,response,offs,n)

[2]

[3] n<<-n

[4] a1<<-a1

[5] a0<<-a0

[6] psi=1 #fest gewaehlte a priori Parameter

[7] nu1=0

[8] response<<-response

[9] Modell<-glm(model1,family=neg.bin(theta=1/a1),x=T) #NegBinGLM M1

[10] coef<-matrix(Modell$coef,ncol=1) #Vektor der Parameterschaetzer

[11] eta<-Modell$x%*%coef #linearer Praediktor

[12] mu<<-care.exp(offs + eta) #Link

[13] deta<- 1/mu #deta=d(eta)/d(mu)

[14] v <- mu * (1 + mu*a1) #Varianzfunktion

[15] w<-1/(deta^2 * v) #Gewichtung

[16] w<-w/sum(w)

[17] z<-eta + (response - mu) *deta #adjusted dependent variable z

[18] xbar<-t((Modell$x)[1:length(response),-1]) %*%w


[20] xbar2<-t(((Modell$x)[1:length(response),-1])^2)%*%w

[21] s2x<-xbar2-xbar^2 #gewichtete Stichprobenvarianzen der Kovariablen

[22] zbar<-as.numeric(weighted.mean(z,w)) #gewichteter


[24] zbar2<-t(z^2)%*%w

[25] s2z <- zbar2 - zbar^2 #gewichtete Stichprobenvarianz der adjusted


[27] varz<- as.numeric(s2z)

[28] m<-ncol(as.matrix((Modell$x)[1:length(response),-1]))

[29] pmean<-c(nu1 + zbar,rep(0,times=m)) #a priori Erwartungswert

[30] #Erzeugung der Matrix U

[31] U<-diag(c(psi^2,rep(phi^2,times=m)),nrow=m+1)

[32] #Erzeugung der Matrix Q

[33] Q<-diag(1/sqrt(s2x) ,nrow=length(s2x))

Anhang 95

[34] Q<-rbind(t(-xbar/sqrt(s2x)),Q)

[35] Qcol1<-c(1,rep(0,times=m))

[36] Q<-sqrt(varz) * cbind(Qcol1,Q)

[37] pvar<- Q %*% U %*% t(Q) # a priori Kovarianzmatrix

[38] Summa<-summary(Modell,dispersion=1)

[39] CovMatrix<-Summa$cov.unscaled #Kovarianzmatrix der Parameterschaetzer

[40] V<-CovMatrix

[41] A<-solve(V) #asymptotische Fisher-Informationsmatrix

[42] thetahat<-matrix(Modell$coef,ncol=1) #Vektor der Parameterschaetzer

[43] G <- solve(pvar) #Inverse der a priori Kovarianzmatrix

[44] H<-solve(A + G)

[45] I<-diag(rep(1,times=m+1),nrow=m+1) #Einheitsmatrix (Rang=m+1)

[46] logdetW<-sum(log(eigen(pvar)$values))

[47] C<- G %*% (I - H %*% (2*I - A %*% H) %*% G)

[48] E1<- -(t(thetahat - pmean))%*% C %*% (thetahat - pmean) - logdetW

[49] - sum(log(eigen(A+G)$values))

[50]

[51] #Vergleichsmodell: NB-Modell M0

[52]

[53] Modell0<-glm(model0,family=neg.bin(theta=1/a0),x=T) #NegBinGLM M0



[56] mu0<<-care.exp(offs + eta0) #Link

[57] deta0<- 1/mu0 #deta=d(eta)/d(mu)

[58] v0 <- mu0 * (1 + mu0*a0) #Varianzfunktion

[59] w0<-1/(deta0^2 * v0) #Gewichtung

[60] w0<-w0/sum(w0)

[61] z0<-eta0 + (response - mu0) *deta #adjusted dependent variable z

[62] xbar0<-t((Modell0$x)[1:length(response),-1]) %*%w0


[64] xbar20<-t(((Modell0$x)[1:length(response),-1])^2)%*%w0

[65] s2x0<-xbar20-xbar0^2 #gewichtete Stichprobenvarianzen der Kovariablen



[68] zbar20<-t(z0^2)%*%w0



96 Anhang


[72] m0<-ncol(as.matrix((Modell0$x)[1:length(response),-1]))

[73] pmean0<-c(nu1 + zbar0,rep(0,times=m0)) #a priori Erwartungswert


[75] U0<-diag(c(psi^2,rep(phi^2,times=m0)),nrow=m0+1)


[77] Q0<-diag(1/sqrt(s2x0) ,nrow=length(s2x0))

[78] Q0<-rbind(t(-xbar0/sqrt(s2x0)),Q0)

[79] Qcol10<-c(1,rep(0,times=m0))

[80] Q0<-sqrt(varz0) * cbind(Qcol10,Q0)

[81] pvar0<- Q0 %*% U0 %*% t(Q0) #a priori Kovarianzmatrix

[82] Summa0<-summary(Modell0,dispersion=1)

[83] CovMatrix0<-Summa0$cov.unscaled #Kovarianzmatrix der Parameterschaetzer

[84] V0<-CovMatrix0

[85] A0<-solve(V0) #asymptotische Fisher-Informationsmatrix


[87] G0 <- solve(pvar0) #Inverse der a priori Kovarianzmatrix

[88] H0<-solve(A0 + G0)

[89] I0<-diag(rep(1,times=m0+1),nrow=m0+1) #Einheitsmatrix (Rang=m0+1)

[90] logdetW0<-sum(log(eigen(pvar0)$values))

[91] C0<- G0 %*% (I0 - H0 %*% (2*I0 - A0 %*% H0) %*% G0)

[92] E0<- -(t(thetahat0 - pmean0))%*% C0 %*% (thetahat0 - pmean0) - logdetW0

[93] - sum(log(eigen(A0+G0)$values))

[94] LogLikeM1<-LogLikeNB(response,a1,n,mu) #Log-Likelihood Modell 1

[95] LogLikeM0<-LogLikeNB(response,a0,n,mu0) #Log-Likelihood Modell 0

[96] BF2ln<-2*(LogLikeM1 - LogLikeM0) + (E1 - E0) #skalierter Bayes-Faktor

[97] E<-E1 - E0

[98] list(BF2ln=BF2ln,model1=model1,model0=model0,

[99] LogLikeM1=LogLikeM1,LogLikeM0=LogLikeM0,E=E) #Ausgabe

[100]

Anhang 97

In der Funktion BFnb2ln benotigt man zur Berechnung der Log-Likelihoods (Zeile 94 und 95)

folgendes Programm.

[1] LogLikeNB<-function(y,a,n,mu)

[2]

[3] a<<-a

[4] ergebnis<-rep(0,n)

[5] for(i in 1:n)

[6]

[7] ergebnis[i]<-Summand1(y[i],a) + y[i]*log(mu[i])

[8] - (y[i] + (1/a))*log(1 + a*mu[i])

[9]

[10] ergebnissum<-sum(ergebnis)

[11] return(ergebnissum)

[12]

In der Funktion LogLikeNB wird zudem in Zeile 7 die Funktion Summandl benutzt.

[1] Summandl<-function(y,a)

[2]

[3] if (y == 0) erg<-0

[4] if(y == 1) erg<-0

[5] if (y > 1)

[6] erg<-rep(0,y)

[7] for ( j in 2:y)

[8]

[9] erg[j] <- log(a * (j-1) + 1)

[10]

[11]

[12] ergsum<-sum(erg)

[13] return(ergsum)

[14]

Fur den Vergleich eines NB-Modells mit einem Poissonmodell sind in der Funktion

BFnb2ln die Zeilen 53 und 95 durch


[95] LogLikeM0<- sum( response * (offs+eta0) - mu0)

zu ersetzen.

98 Anhang

Will man den Vergleich zu einem NB-Nullmodell, d.h. zu einem NB-Modell nur mit Intercept

ziehen, muss man in der Funktion BFnb2ln die Zeilen 51 bis einschließlich 93 durch folgende

Programmzeilen ersetzen.

[51] #NB-NULLMODELL (nur mit Intercept):

[52]

[53] Modell0<-glm(model0,family=neg.bin(1/a0),x=T) #NegBinGLM M0



[56] mu0<-care.exp(offs + eta0) #Linkfunktion

[57] deta0<-1/mu0 #deta=d(eta)/d(mu)

[58] v0<-mu0 #Varianzfunktion

[59] w0<-1/(deta0*v0) #Gewichtung

[60] w0<-w0/sum(w0)

[61] z0<-eta0 + (response - mu0) *deta0 #adjusted dependent variable z


[63] #Stichprobenerwartungswert der adjusted dependent variable

[64] zbar20<-t(z0^2)%*%w0




[68] pmean0<- nu1 + zbar0 #a priori Erwartungswert

[69] Summa0<-summary(Modell0,dispersion=1)

[70] CovMatrix0<-Summa0$cov.unscaled #Kovarianzmatrix der Parameterschaetzer

[71] V0<-CovMatrix0

[72] A0<-1/V0 #asymptotische Fisher-Informations"matrix"



[75] Qcol10<-1

[76] Q0<-sqrt(varz0) *Qcol10


[78] U0<-psi^2

[79] pvar0<- Q0 * U0 * Q0 #a priori Kovarianzmatrix

[80] G0<- 1/pvar0 #Inverse der a priori Kovarianzmatrix

[81] H0<-1/(A0 + G0)

[82] C0<- G0 * (1 - H0 * (2 - A0 * H0) * G0)

[83] E0<- -(thetahat0 - pmean0)* C0 * (thetahat0 - pmean0) - log(pvar0)

[84] - log(A0 + G0)

Anhang 99

2. Die SAS-Prozedur sql zur Verdichtung des KH-Datensatzes

[1] proc sql;

[2] create table Datensatzname as

[3] select

[4] SFKl, Geschlecht, Brandrabatt, Haltedauer, Erwerbsalter, AltVN, dreij,

[5] Fahrerkreis, kmKlasse, KFZ, Typklasse, KfzKwKlasse,

[6] sum(Jahreseinheiten) as SumJE, sum(AnzSchaeden) as SumAnzSchaeden

[7] from Quelldatensatz

[8] group by

[9] SFKl, Geschlecht, Brandrabatt, Haltedauer, Erwerbsalter, AltVN, dreij,

[10] Fahrerkreis, kmKlasse, KFZ, Typklasse, KfzKwKlasse;

[11] quit;

3. Splus-Programm zur Erstellung der Haupteffekte-Grafiken

[1] main<-

[2] function(y = SumAnzSchaeden, n = SumJE, x = Kovariable,

[3] label = "Ueberschrift")

[4]

[5] sy <- y/n #Berechnung der Schadenhaeufigkeit pro Zelle

[6] temp <- tapply(sy, x, mean) #Mittelwert der Schadenhaeufigkeit fuer

[7] #jede Auspraegung der betrachteten Kovariable

[8] #Erzeugung der Grafik

[9] plot(log(temp), axes = F, pch = "x", xlab = "", ylab =

[10] "empirical log mean", type = "b", lty = 3, col = 1, lwd

[11] = 1, sub = label,main = label)

[12] axis(2)

[13] axis(1, at = 1:length(temp), labels = levels(x))

[14]

In dieser Funktion dient die Zeile 13 der Beschriftung der x-Achse in der Grafik. Dabei ist zu

beachten, dass die Option labels=levels(x) nicht fur Kovariablen geeignet ist, deren Aus-

pragungen numerische Bezeichnungen haben. Nehmen wir beispielsweise an, der betrachtete

Datensatz enthalt die Kovariable Geschlecht, deren Auspragungen mit 1 fur mannlich und 2

fur weiblich kodiert sind. In diesem Fall ist es notwendig, die Zeile 13 durch eine der folgenden

Zeilen zu ersetzen.

[13] axis(1, at = 1:length(temp), labels = c("1", "2"))

[13] axis(1, at = 1:length(temp), labels = c("maennlich", "weiblich"))

Literaturverzeichnis

[Bron 99] Bronstein, I.N., Semendjajew, K.A., Musiol, G., Muhlig, H. (1999): Taschen-

buch der Mathematik, 4. Auflage, Harri Deutsch-Verlag

[Berg 85] Berger, J.O. (1985): Statistical Decision Theory and Bayesian Analysis, New

York: Springer-Verlag

[Berg 87] Berger, J.O. and Sellke, T. (1987): Testing a point null hypothesis: the ir-

reconcilability of P values and evidence, Journal of the American Statistical

Association, 82, 112-122

[Bick 77] Bickel, P.J., Doksum, K.A. (1977): Mathematical Statistics: basic ideas and

selected topics, Holden-Day

[Brui 70] de Bruijn, N.G. (1970): Asymptotic Methods in Analysis, Amsterdam: North-

Holland

[Came 98] Cameron, A.C., Trivedi, P.K. (1998): Regression analysis of count data, Cam-

bridge: University Press

[Case 90] Casella, G., Berger, R.L. (1990): Statistical Inference, Wadsworth &

Brooks/Cole

[Fahr 94] Fahrmeir, L., Tutz, G. (1994): Multivariate statistical modelling based on ge-

neralized linear models, New York, Springer-Verlag

[Fahr 99] Fahrmeir, L., Kunstler, R., Pigeot, I., Tutz, G. (1999): Statistik, der Weg zur

Datenanalyse, zweite verbesserte Auflage, Springer-Verlag

[Fran 01] Franzmann, A. (2001): Multivariate Klassifikation in der Kraftfahrzeughaft-

pflichtversicherung, Diplomarbeit an der Johannes Gutenberg-Universitat

Mainz

[Jeff 61] Jeffreys, H. (1961): Theory of Probability (3rd. ed.), Oxford: University Press

100

Literaturverzeichnis 101

[John 93] Johnson, N.L., Kotz, S., Kemp, A.W. (1993): Univariate discrete distributi-

ons, 2nd edition, New York, John Wiley & Sons

[Krab 83] Krabs, W. (1983): Einfuhrung in die lineare und nichtlineare Optimierung,

Stuttgart, Teubner-Verlag

[Kass 92a] Kass, R.E. and Wassermann, L. (1992): A reference Bayesian test for ne-

sted hypotheses with large samples, Technical Report no. 576, Department of

Statistics, Carnegie-Mellon University.

[Kass 92b] Kass, R.E. and Vaidyanathan, S. (1992): Approximate Bayes factors and or-

thogonal parameters, with application to testing equality of two Binomial pro-

portions, J. Royal Statist. Soc. B., 54, 129-144

[Lawl 87] Lawless, J.F (1987): Negative binomial and mixed Poisson regression, The

Canadian Journal of Statistics, Vol.15, No. 3, 209-225

[McCu 89] McCullagh, P., and Nelder, J.A. (1989): Generalised Linear Models, 2nd edi-

tion, London: Chapman and Hall

[Rade 97] Rade, L., Westergreen, B. (1997): Springers Mathematische Formeln, 2. Auf-

lage, Springer-Verlag

[Raft 93] Raftery, Adrian E. (1993): Bayesian model selection in structural equation

models. In Testing Structural Equation Models (K.A. Bollen und J.S. Long),

Beverly Hilly: Sage, Seiten 163-180

[Raft 94] Raftery, Adrian E.(August 1993; Revised March 1994): Approximate Bayes

Factors and Accounting for Model Uncertainty in Generalized Linear Models,

Technical Report no. 255, Department of Statistics, University of Washington

[Raft 96] Raftery, Adrian E. (1996): Approximate Bayes factors and accounting for

model uncertainty in generalised linear models, Biometrika (1996), 83, 2, pp.

251-266, Printed in Great Britain

[Siko 02] Sikora, I.-D. (2002): Quantifizierung von Uberdispersion, Diplomarbeit an der

Technischen Universitat Munchen

[Rao 73] Rao, C.R. (1973): Linear statistical inference and its applications, 2nd edition,

New York, Wiley

[Tier 86] Tierney, L., and Kadane, J.B. (1986): Accurate approximations for posterior

moments and marginal densities, J. Amer. Statist. Ass., 81, 82-86

Abbildungsverzeichnis

2.1 Darstellung der Poisson-Wahrscheinlichkeitsfunktion fur verschiedene λ-Werte . . 4

2.2 NB-Verteilungsdichten fur verschiedene Erfolgswahrscheinlichkeiten p . . . . . . . 29

2.3 NB-Verteilungsdichten fur verschiedene r-Werte . . . . . . . . . . . . . . . . . . . 30

3.1 RPO10(φ; β2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.2 RPO10(φ; 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3.3 R(φ) und S(φ) als Funktionen von φ . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.1 Uberdispersion im Poissonmodell 2 der Patent-Daten (— Winkelhalbierende, – – –

geglatteter Zusammenhang) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

4.2 Log-Bayes-Faktoren als Funktion des zum jeweiligen NB-Modell gehorigen Uber-

dispersionsparameters fur die NB-Poisson-Modellvergleiche der Patent-Daten in

Abhangigkeit vom gewahlten Hyperparameter φ . . . . . . . . . . . . . . . . . . 79

4.3 Gammaverteilungsdichte mit Erwartungswert 0.25 und 0.1 (Varianz = 0.5) . . . 85

102

Tabellenverzeichnis

3.1 Tabelle fur die Interpretation der Bayes-Faktoren . . . . . . . . . . . . . . . . . . 45

3.2 Fur die Approximation der skalierten Bayes-Faktoren notwendige Terme und de-

ren Erlauterung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.3 Approximationen zur Berechnung des skalierten Bayes-Faktors . . . . . . . . . . 56

4.1 Zu vergleichenden Modelle zur Anpassung der Patent-Daten . . . . . . . . . . . . 74

4.2 Splus-summary des Poissonmodells 1 fur die Patent-Daten . . . . . . . . . . . . . 76

4.3 Splus-summary des Poissonmodells 2 fur die Patent-Daten . . . . . . . . . . . . . 76

4.4 Log-Likelihoods und zugehorige ML-Schatzer der Dispersionsparameter fur die

NB-Poisson-Modellvergleiche der Patent-Daten . . . . . . . . . . . . . . . . . . . 81

4.5 Maximale Log-Bayes-Faktoren und deren Uberdispersionsparameter fur die NB-

Poisson-Modellvergleiche der Patent-Daten in Abhangigkeit vom gewahlten Hy-

perparameter φ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

4.6 Durchschnittliche Log-Bayes-Faktoren fur die NB-Poisson-Modellvergleiche der

Patent-Daten in Abhangigkeit von φ und der gewahlten a priori Verteilung fur

den Dispersionsparameter im jeweiligen NB-Modell . . . . . . . . . . . . . . . . . 84

4.7 Splus-summary des NB-Modells 4 fur die Patent-Daten . . . . . . . . . . . . . . . 86

103

Technische Universität M ünchen - mediaTUM · Technische Universität M ünchen Zentrum...

Documents

Transcript of Technische Universität M ünchen - mediaTUM · Technische Universität M ünchen Zentrum...