Generalisierte Lineare Modelle

(1)

Herwig FRIEDL

Institut f¨ ur Statistik Technische Universit¨ at Graz

September 2013

(2)

Inhaltsverzeichnis

1 Transformation auf Normalverteilung 1

1.1 Box-Cox Transformationsfamilie . . . 1

1.2 Maximum-Likelihood Sch¨atzung . . . 2

1.3 Beispiel . . . 4

2 Die Exponentialfamilie 9 2.1 Maximum Likelihood Sch¨atzung . . . 12

2.2 Mitglieder der Exponentialfamilie . . . 12

2.3 Die Quasi-Likelihood Funktion . . . 15

2.3.1 Quasi-Dichten . . . 18

3 Das Generalisierte Lineare Modell 19 3.1 Maximum-Likelihood-Sch¨atzungen . . . 20

3.2 Asymptotische Eigenschaften des Sch¨atzers . . . 23

3.3 Pearson Statistik . . . 23

3.4 Eigenschaften von Score und Quasi-Score . . . 24

3.5 Deviance und Quasi-Deviance . . . 25

3.6 Parametertests . . . 27

3.7 Beispiel: Konstante Varianz . . . 27

3.8 Beispiel: Gamma-Varianz . . . 28

4 Logistische Regression 37 4.1 Toleranzverteilungen – Linkfunktionen . . . 38

4.1.1 Beispiel (Venables & Ripley, Seite 218) . . . 39

4.2 Interpretation der Parameter . . . 43

4.2.1 Beispiel (Agresti, Seite 87) . . . 44

4.3 Logit-Modelle . . . 46

4.3.1 Beispiel . . . 46 i

(3)

4.4 Komplexe Logistische Modelle . . . 49

4.4.1 Beispiel . . . 49

5 Loglineare Modelle 51 5.1 Poisson Loglineare Modelle f¨ur Anzahlen . . . 51

5.1.1 Beispiel . . . 51

5.2 Loglineare und Multinomiale Response Modelle . . . 57

5.2.1 Die Multinomialverteilung . . . 57

5.2.2 Vergleich von Poisson-Erwartungen . . . 59

5.2.3 Multinomiale Responsemodelle . . . 60

5.2.4 Beispiel . . . 61

5.3 Zweidimensionale Kontingenztafeln . . . 64

5.3.1 Unabh¨angigkeitsmodell . . . 64

5.3.2 Saturiertes (volles) Modell . . . 65

5.3.3 Interpretation der Parameter . . . 67

5.3.4 Beispiel . . . 69

6 Modelle mit zufälligen Effekten 71 6.1 Zufällige Prädiktoren . . . 71

6.2 EM Sch¨atzer . . . 72

6.2.1 Beispiel: Endliche diskrete Mischungen . . . 74

6.3 Uberdispersionsmodelle . . . .¨ 77

6.3.1 Normalverteilte zuf¨allige Eﬀekte . . . 77

6.3.2 Zuf¨allige Eﬀekte aus unbekannter Verteilung . . . 79

6.3.3 Pr¨adiktionen bei NPML Sch¨atzung . . . 80

6.3.4 Beispiel: Matched Pairs . . . 81

(4)

Kapitel 1

Transformation auf Normalverteilung

Die statistische Analyse von Daten basiert häufig auf der Annahme, dass diese normal- verteilt sind und konstante Varianz widerspiegeln. Falls die Daten diese Annahme nicht unterstützen, besteht die Möglichkeit der Verwendung einer Transformation, um dadurch eine bessere Approximation zu einer konstanten Varianz zu erzielen. Dann könnten auch klassische Methoden wie die Varianzanalyse oder die Lineare Regression auf solche Daten angewendet werden.

1.1 Box-Cox Transformationsfamilie

Die Verwendbarkeit der Normalverteilung wird erweitert, indem diese in eine größere Familie von Verteilungsfunktionen eingebettet wird, der Box-Cox Transformationsfamilie (Box und Cox, 1964). Deren allgemeine Form kann für eine positive Response y > 0 repräsentiert werden durch

y(λ) =



 y^λ−1

λ , fallsλ̸= 0, logy, fallsλ= 0,

(1.1)

wobeiλ den Parameter der Transformation bezeichnet. Spezialfälle in dieser Familie sind y(−1) = 1−1/y undy(1) =y−1. Darüberhinaus strebt fürλ→0,y(λ)→logy, so dass y(λ) eine stetige Funktion in λ ist.

Für Daten (y_i, x_i), i = 1, . . . , n, nehmen wir nun an, dass es einen Wert von λ gibt, für den y_i(λ) einer Normalverteilung mit identischer Varianz genügt, d.h.

yi(λ)^ind∼ N(µi(λ), σ²(λ)).

Unter dieser Annahme kann man auf die Dichtefunktion der originalen Dateny schließen.

1

(5)

Diese ist gerade

f(y, λ, µ(λ), σ(λ)) = 1

√2πσ²(λ)exp (

−(y(λ)−µ(λ))² 2σ²(λ)

) d dyy(λ)

. (1.2)

Die Verwendung der Transformation (1.1) ergibt

f(y, λ, µ(λ), σ(λ)) =











√ 1

2πσ²(λ)exp (

−

((y^λ−1)/λ−µ(λ))2

2σ²(λ)

)

y^λ⁻¹, fallsλ̸= 0,

√ 1

2πσ²(λ)exp (

−(logy−µ(λ))² 2σ²(λ)

)

y⁻¹, fallsλ= 0.

(1.3)

In Regressions Analysen verwendet man häufig ein simples lineares Modell y_i(λ)înd∼ N(x^⊤_i β(λ), σ²(λ)),

wobeiβ(λ) = (β₁(λ), . . . , β_p(λ))^⊤denp×1 Vektor der unbekannten Parameter bezeichnet.

Für λ ̸= 0 definiert man den reparameterisierten Vektor β = (β₁, . . . , β_p)^⊤ bezüglich y^λ, mit Intercept β1 = 1 +λβ1(λ) und Effekten βj =λβj(λ), j = 2, . . . , p, sowie σ =|λ|σ(λ).

Damit kann die Dichte (1.3) umgeschrieben werden zu f(y, λ, β, σ) = 1

√2πσ² exp (

−(y^λ−x^⊤β)² 2σ²

)

|λ|y^λ⁻¹.

F¨urλ = 0 verwendeβ_j =β_j(λ),j = 1, . . . , p, und σ=σ(λ) als Parameter bez¨uglich eines linearen Modells f¨ur logy. Dann kann man (1.2) schreiben als

f(y,0, β, σ) = 1

√2πσ² exp (

−(logy−x^⊤β)² 2σ²

) y⁻¹.

1.2 Maximum-Likelihood Sch¨ atzung

Die Log-Likelihood Funktion der n Beobachtungen (y_i, x_i) basierend auf die untransfor- mierte Dichte (1.2) ist (auf der originalen Datenskala)

l(λ, µ(λ), σ(λ)) =

− n 2log

(

2πσ²(λ)

)− 1 2σ²(λ)

∑n i=1

(

yi(λ)−µi(λ) )2

+

∑n i=1

log d

dy_iyi(λ)

. (1.4) Mit der Box-Cox Transformation (1.1) ist die Summe der Logarithmen der Jacobi Terme

∑n i=1

log d

dy_iy_i(λ) =

{

(λ−1)∑

ilogyi, fallsλ ̸= 0,

−∑

ilogyi, fallsλ = 0.

(6)

1.2. MAXIMUM-LIKELIHOOD SCH ÄTZUNG 3 Bezüglich der Parameter (β, σ) erhält man

l(λ, β, σ) =

{−ⁿ₂ log (2πσ²)− _2σ¹2

∑

i

(y^λ_i −x^⊤_i β)2

+nlog|λ|+ (λ−1)∑

ilogyi, falls λ̸= 0,

−ⁿ₂ log (2πσ²)− _2σ¹2

∑

i

(logyi −x^⊤_i β)2

−∑

ilogyi falls λ= 0. (1.5) Für einen festen Wert vonλlösen die Maximum-Likelihood Schätzer ˆβ_λ und ˆσ_λ basierend auf (1.5) die Schätzgleichungen

∂l(λ, β, σ)

∂β =

{ 1 σ²

∑

ixi(y^λ_i −x^⊤_i β) = 0, falls λ̸= 0,

1 σ²

∑

ixi(logyi−x^⊤_i β) = 0, falls λ= 0,

∂l(λ, β, σ)

∂σ =

{−ⁿ_σ +_σ¹3

∑

i(y_i^λ−x^⊤_i β)² = 0, falls λ̸= 0,

−ⁿ_σ +_σ¹3

∑

i(logy_i−x^⊤_i β)² = 0, falls λ= 0, und sind deshalb gleich

βˆ_λ = {

(X^⊤X)⁻¹X^⊤y^λ, falls λ̸= 0, (X^⊤X)⁻¹X^⊤logy, falls λ= 0, ˆ

σ_λ² = 1

nSSE_λ( ˆβ_λ) = {1

n

∑

i(y^λ_i −x^⊤_i βˆλ)², falls λ̸= 0,

1 n

∑

i(logyi−x^⊤_i βˆλ)², falls λ= 0,

wobeiy^λ(resp. logy) elementeweise gerechnet sind und SSE_λ( ˆβ_λ) die Fehlerquadratsumme von y^λ (resp. logy) an der Stelle ˆβλ f¨ur ein festes λ bezeichnet. Bemerke, dass wegen der obigen Reparameterisierung die Fehlerquadratsumme SSE_λ( ˆβ_λ) in λ = 0 unstetig ist.

Substituiert man beide Schätzer in der Log-Likelihood Funktion (1.5) und lässt darin alle konstanten Terme weg, so erhält man die Profile (Log-) Likelihood Funktion

pl(λ) =l(λ,βˆ_λ,σˆ_λ) =

{−ⁿ₂ log SSE_λ( ˆβ_λ) +nlog|λ|+ (λ−1)∑

ilogy_i, falls λ̸= 0,

−ⁿ₂ log SSE₀( ˆβ₀)−∑

ilogy_i, falls λ= 0.(1.6) F¨ur λ = 1 resultiert z.B. pl(1) = −(n/2) log SSE₁( ˆβ₁) = nlog(∑n

i=1(y_i−x^⊤_i βˆ₁)² )₋1/2

. Obwohl SSE_λ(·) in λ = 0 unstetig ist, ist die Proﬁle Likelihood Funktion pl(λ) stetig.

Wegen

pl(λ) = −n 2 log

∑n i=1

(y_i^λ−x^⊤_i βˆ_λ)²

λ² + (λ−1)

∑n i=1

logy_i

= −n 2 log

∑n i=1

((y_i^λ−1)/λ−x^⊤_i β(λ))ˆ ²+ (λ−1)

∑n i=1

logy_i ist limλ→0pl(λ) =pl(0).

(7)

Wegen{2(pl(ˆλ)−pl(λ₀))∼χ²₁}beinhaltet ein approximatives Konfidenzintervall für den Parameter λ all jene Werte von λ₀, für die pl(λ₀) innerhalb von ¹₂χ²₁₋_α;1 Einheiten vom Funktionsmaximum entfernt ist (χ²_0.95;1= 3.8415,χ²_0.99;1 = 6.6349).

Ein wichtiger Aspekt dieser Variablentransformation ist, dass auf der transformierten Skala das Modell die Variation bezüglich des Erwartungswertes der (auf Normalverteilung) transformierten Variablen repräsentiert, während auf der Originalskala das Modell die Variation bezüglich des Medians der originalen Variablen darstellt. Dies sieht man am einfachsten für die Log-Transformation (λ= 0). Seien logy_i ∼N(x^⊤_i β, σ²), dann gilt

median(logy_i) = x^⊤_i β, E(logy_i) = x^⊤_i β, var(logy_i) = σ².

Die originalen Beobachtungen y_i unterliegen selbst einer Lognormalverteilung mit median(y_i) = exp(x^⊤_i β),

E(y_i) = exp(x^⊤_i β+σ²/2), var(y_i) = (

exp(σ²)−1)

exp(2x^⊤_i β+σ²).

Dies bedeutet, dass das additive Modell für den Erwartungswert (und daher auch für den Median) der logy_i ein multiplikatives Modell für die Mediane und Erwartungswerte dery_i ist. Für die Erwartungswerte von y_i wurde der Intercept Term um σ²/2 verschoben und die Varianzen sind nicht mehr konstant.

F¨ur eine Transformation y(λ) = y^λ mit λ̸= 0, also f¨ury_i^λ ∼N(µi, σ²), folgt median(y_i) = µ^1/λ_i ,

E(y_i) ≈ µ^1/λ_i (

1 +σ²(1−λ)/(2λ²µ²_i)) , var(y_i) ≈ µ^2/λ_i σ²/(λ²µ²_i).

Wiederum ist die oﬀensichtliche Unstetigkeit zwischenλ= 0 undλ̸= 0 in der Verwendung von y^λ anstelle von (y^λ−1)/λ begr¨undet.

1.3 Beispiel

n = 31 Black Cherry Trees; verwendbares Holzvolumen V in feet³ (1 foot = 30.48 cm), Baumhöhe H in feet, DurchmesserD in inches (1 inch = 2.54 cm) auf einer Höhe von 4.5 feet über dem Boden. Das Modell sollte das verwendbare Holzvolumen V aus den leicht zu messenden Größen H und D vorhersagen.

> %in%

> attach(trees); plot(H, V); plot(D, V) # Kr¨ummung (falsche Skala?)

> mod <- lm(V ~ D + H)

(8)

1.3. BEISPIEL 5

8 10 12 14 16 18 20

D

204060

V

65 70 75 80 85

H

204060

V

0 5 10 15 20 25 30

Index

0.050.100.150.20

leverages

8 10 12 14 16 18 20

D

-505

r

Abbildung 1.1: Oben: Volumen gegen Durchmesser (links) und Volumen gegen H¨ohe (rechts). Unten: Diagonalelemente der Hatmatrix (links) und Residuen gegen Durchmesser (rechts) unter dem linearen Modell f¨ur V.

Coefficients:

(Intercept) H D

-57.98766 0.3392512 4.708161

Degrees of freedom: 31 total; 28 residual Residual standard error: 3.881832

> plot(lm.influence(mod)$hat, ylab = "leverages") # 2 Hebelpunkte

> abline(2*mod$rank/length(V), 0); abline(v = c(20, 31))

> plot(residuals(mod), D); abline(0, 0) # Senke in der Mitte

> library(MASS, help=T)

> boxcox(V ~ D + H, lambda = seq(0.0, 0.6, len=18)) # max = -76.08

Das Maximum der Profile Likelihood Funktion tritt in der Nähe von λ = 0.3 auf. Das approximative 95% Konfidenzintervall ist sehr klein, etwa (0.12,0.49), und beinhaltet weder die Null noch die Eins. Es liegt in der Natur einer Volumenmessung, dass sich diese kubischbezüglich den linearen Größen Höhe und Durchmesser verhält. Daher erscheint es auch sinnvoll, die Kubikwurzel des Volumens als Response zu verwenden.

> plot(D, V**(1/3)) # Kr¨ummung beseitigt

(9)

0.0 0.1 0.2 0.3 0.4 0.5 0.6 lambda

-80-79-78-77-76

log-Likelihood

95%

8 10 12 14 16 18 20

D

2.53.03.54.0

V^(1/3)

20 40 60 80

mu^3

204060

V

20 40 60 80

mu^3 * (1 + (3 * s2)/mu^2)

20406080

mu^3

Abbildung 1.2: Oben: Profile Likelihood Schätzung vonλbei der Box-Cox Transformation (links) und Durchmesser gegen diey_i^1/3 (rechts). Unten: Gefittete Mediane gegen originale Beobachtungen V (links) und gefittete Mediane gegen gefittete Means (rechts).

> mod1 <- lm(V**(1/3) ~ H + D) Coefficients:

(Intercept) H D

-0.08538814 0.01447169 0.1515163

Die gefitteten Mediane von V unter diesem obigen Modell fürµ_1/3 = E(y^1/3) (mit festem λ = 1/3) sind ˆµ³_1/3, und die gefitteten Erwartungswerte sind etwa ˆµ³_1/3(1 + 3ˆσ²_1/3/µˆ²_1/3).

Diese k¨onnen mit den Beobachtungen y verglichen werden.

> mu <- fitted(mod1)

> s2 <- deviance(mod1)/mod1$df.residual

> plot(mu**3, V) # Median Modell

> plot(mu**3*(1+3*s2/mu**2), mu**3) # Median/Mean

Andere technische Überlegungen ergeben alternative Modelle. Die unerwünschte Krüm- mung in der Abbildung 1.1 kann auch durch logarithmische Transformation aller Variablen

(10)

1.3. BEISPIEL 7

2.2 2.4 2.6 2.8 3.0

log(D)

2.53.03.54.0

log(V)

-0.2 -0.1 0.0 0.1 0.2

lambda

-80-79-78-77-76-75

log-Likelihood

95%

Abbildung 1.3: Links: logV linear abhängig von logD. Rechts: Profile Likelihood Funktion für V modelliert durch logH+ logD.

entfernt werden. Dies legt eine Regression von auf log(D) und log(H) nahe. Soll man jetzt jedoch auf der log(V) Achse modellieren?

> plot(log(D), log(V)); lm(log(V) ~ log(D) + log(H)) Coefficients:

(Intercept) log(D) log(H) -6.631617 1.98265 1.117123

> boxcox(V ~ log(H) + log(D), lambda = seq(-0.25, 0.25, len = 10))

Die Profile Likelihood Schätzung bei der Box-Box Transformation liefert ein Maximum bei −0.07 und ein 95% Konfidenzintervall von (−0.24,0.11), welches zwar die Null (logarithmische Transformation), aber nicht mehr die Kubikwurzeltransformation λ = 1/3 oder die Identitätλ= 1 beinhaltet.

Beide Modelle liefern ann¨ahernd dieselben Maxima der Proﬁle Likelihood Funktionen.

Welches der beiden ist nun dasbessere? Wir k¨onnen sie durch einen Likelihood Quotienten Test miteinander vergleichen. Dazu werden die Modelle eingebettet in die Familie

V^∗ ∼ N(β₁+β₂H^∗+β₃D^∗, σ²) V^∗ = (V^λ^V −1)/λ_V

H^∗ = (H^λ^H −1)/λH

D^∗ = (D^λ^D−1)/λ_D

Wir vergleichen nun die maximale Likelihood bezüglich λ_V = 1/3, λ_H = λ_D = 1 mit der maximalen Likelihood des Modells λV = λH = λD = 0. Alle übrigen Parameter sind hierbei nuisance Parameter und die Likelihood Funktion wird über diese in beiden Modellen maximiert.

> boxcox(V ~ H + D, lambda = 1/3, plotit=FALSE)

(11)

$x: 0.3333333 $y: -76.12156

> boxcox(V ~ log(H) + log(D), lambda = 0, plotit=FALSE)

$x: 0 $y: -75.33877

Die doppelte Differenz der beiden Maxima beträgt nur −2∗(−76.122 + 75.339) = 1.566, was nicht signifikant verglichen mitχ²₃ist. Daher können wir auch nicht überzeugend eines der beiden Modelle wählen.

Bemerke aber, dass der Koeffizient für logH nahe bei Eins (1.117) liegt und der zu logD fast Zwei (1.983) ist. Nimmt man an, dass man einen Baum durch einen Zylinder oder Kegel beschreiben kann, so wäre sein Volumen πd²h/4 (Zylinder) oder πd²h/12 (Kegel).

In beiden F¨allen h¨atte man ein Modell der Form

logV =c+ logH+ 2 logD

mit c= log(π/4) (Zylinder) oder c= log(π/12) (Kegel). Jedoch beziehen sich die obigen Uberlegungen auf Messungen auf der gleichen Skala. Wir konvertieren daher zuerst¨ Dvon inches auf feet, d.h. wir betrachten D/12. Als Intercept resultiert −1.705.

> lm(log(V) ~ log(D/12) + log(H)) Coefficients:

(Intercept) log(D/12) log(H) -1.704917 1.98265 1.117123

Nun ﬁxiert man die beiden Koeﬃzienten und betrachtet das Modell logV =β₁+ logH+ 2 log(D/12)

Hierbei bezeichnet man 1×[logH+2 log(D/12)] alsoﬀset(ein Term mit festem Parameter Eins) und es muss nur noch β₁ gesch¨atzt werden.

> mod3 <- lm(log(V) ~ 1 + offset(2*log(D/12) + log(H)))

> summary(mod3) Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) -1.19935 0.01421 -84.42 <2e-16 ***

Residual standard error: 0.0791 on 30 degrees of freedom

> log(pi/4) [1] -0.2415645

> log(pi/12) [1] -1.340177

Das Holzvolumen eines Baums kann daher eher durch ein Kegelvolumen als durch das eines Zylinders beschrieben werden, hat jedoch ein etwas gr¨oßeres Volumen als ein Kegel.

(12)

Kapitel 2

Die Exponentialfamilie

Beim Linearen Modell (LM) wird angenommen, dass die abhängigen Variablen (Respon- ses) y_i stochastisch unabhängige, normalverteilte Größen sind mit Erwartungen µ_i =x^⊤_i β und konstanter Varianz σ². In manchen Situationen ist die Annahme einer Normalver- teilung sicherlich sehr künstlich und nur schwer zu vertreten. Man denke hierbei nur an Modelle für absolute Häufigkeiten oder relative Anteile. Weiters gibt es datengenerieren- de Mechanismen, die für größere Erwartungswerte auch größere Variabilität induzieren.

Hierzu zählen beispielsweise Modelle für konstante Variationskoeffizienten. Da bei einem LM alle Erwartungswerte auf einer Ebene imp-dimensionalen Raum liegen, ist ein solches LM für nicht-negative oder speziell für binäre Responses sicherlich unpassend.

Wir wollen nun wegen all dieser Schwachstellen die Klasse der Generalisierten Linea- ren Modelle (GLM) betrachten, die gerade bezüglich der oben angeführten Restriktionen flexible Verallgemeinerungen anbieten. So wird beim klassischen GLM statt der Normal- verteilung eine Verteilung aus der einparametrigen Exponentialfamilie angenommen und dadurch die Varianz als Funktion des Erwartungswertes modelliert (V ar(y_i) = V(µ_i)).

Dar¨uberhinaus wird der Erwartungswert nicht ausschließlich direkt linear modelliert, son- dern der lineare Pr¨adiktor η_i =x^⊤_i β entspricht einer bekannten Funktiong(µ_i), der Link- funktion.

Bevor nun die Exponentialfamilie eingef¨uhrt wird, soll an die wesentlichsten Eigenschaften von Score-Funktionen und Informations-Matrizen bei der Maximum-Likelihood Sch¨atzung erinnert werden. Dazu bezeichnef(y, θ) die Wahrscheinlichkeitsfunktion oder Dichtefunk- tion von y mit unbekannten Parametern θ.

Satz 2.1 F¨ur die Ableitung der Log-Likelihood Funktion, d.h. f¨ur die Score-Funktion gilt:

E

(∂logf(y, θ)

∂θ )

= 0, (2.1)

E

(∂logf(y, θ)

∂θ

)2

= E

(

−∂²logf(y, θ)

∂θ²

)

. (2.2)

9

(13)

Mit

∂logf(y, θ)

∂θ = 1

f(y, θ)

∂f(y, θ)

∂θ und

∫

R

f(y, θ)dy= 1 folgt (2.1), denn

E

(∂logf(y, θ)

∂θ

)

= E

(∂f(y, θ)

∂θ

1 f(y, θ)

)

=

∫

R

∂f(y, θ)

∂θ dy= ∂

∂θ

∫

R

f(y, θ)dy= 0. Weiters ergibt sich (2.2) mittels Kettenregel

E (

−∂²logf(y, θ)

∂θ²

)

= E (

−∂²f(y, θ)

∂θ²

1

f(y, θ) +∂f(y, θ)

∂θ

∂f(y, θ)

∂θ

1 f(y, θ)²

)

= −

∫

R

∂²f(y, θ)

∂θ² dy+

∫

R

∂logf(y, θ)

∂θ

∂logf(y, θ)

∂θ f(y, θ)dy

= E

(∂logf(y, θ)

∂θ

)2

.

Diese Likelihood-Eigenschaften gelten natürlich gerade auch für die Exponentialfamilie, die jetzt definiert wird.

Deﬁnition 2.1 Eine Zufallsvariable y sei aus einer Verteilung mit Dichte- oder Wahr- scheinlichkeitsfunktion

f(y, θ) = exp

(yθ−b(θ)

a(ϕ) +c(y, ϕ) )

f¨ur spezielle bekannte Funktionen a(·), b(·) und c(·) mita(ϕ)>0. Kann ϕ als feste Gr¨oße betrachtet werden, so bezeichnet man f(y, θ) als einparametrige Exponentialfamilie mit kanonischem Parameter θ.

Mit (2.1) resultiert bei der Exponentialfamilie E

(∂logf(y, θ)

∂θ

)

= 1

a(ϕ)E(

y−b^′(θ))

= 0 also E(y) = b^′(θ) , und mit (2.2) folgt

E

(∂²logf(y, θ)

∂θ²

) + E

(∂logf(y, θ)

∂θ

)2

=− 1

a(ϕ)b^′′(θ) + 1

a²(ϕ)var(y) = 0.

Diese zwei Identit¨aten ergeben unmittelbar die beiden ersten Momente (Kumulanten) der Exponentialfamilie

E(y) = b^′(θ) var(y) = a(ϕ)b^′′(θ).

(14)

11 Kumulanten h¨oherer Ordnung bestimmt man einfacher mit der Kumulantenerzeugenden Funktion K(t) = logM(t), wobei M(t) die Momentenerzeugende bezeichnet. Die k-te Kumulante κ_k ist nun gegeben durch K^(k)(t)|t=0 und steht mit dem k-ten Moment in einer einfachen Beziehung, denn

κ₁(y) = E(y) κ₂(y) = E(y−µ)² κ₃(y) = E(y−µ)³

κ₄(y) = E(y−µ)⁴−3var²(y). F¨ur die Exponentialfamilie gilt

1 =

∫

R

exp

(yθ−b(θ)

a(ϕ) +c(y, ϕ) )

dy = exp (

−b(θ) a(ϕ)

) ∫

R

exp ( y

a(ϕ) θ+c(y, ϕ) )

dy , woraus

exp (b(θ)

a(ϕ) )

=

∫

R

exp ( y

a(ϕ) θ+c(y, ϕ) )

dy folgt. Die Momentenerzeugende ist daher gegeben durch

M(t) = E(e^ty) = exp (

−b(θ) a(ϕ)

) ∫

R

exp ( y

a(ϕ) (

θ+a(ϕ)t )

+c(y, ϕ) )

dy

= exp (

−b(θ) a(ϕ)

) exp



b (

θ+a(ϕ)t )

a(ϕ)



= exp



b (

θ+a(ϕ)t

)−b(θ) a(ϕ)



,

und als Kumulantenerzeugende Funktion resultiert

K(t) = logM(t) = b

(

θ+a(ϕ)t

)−b(θ)

a(ϕ) .

Die k-te Kumulante von y, κ_k(y), ist somit κ_k(y) = K^(k)(t)|t=0 = a(ϕ)^k⁻¹b^(k)

(

θ+a(ϕ)t)

t=0

=a(ϕ)^k⁻¹b^(k)(θ). (2.3) Sei E(y) = b^′(θ) = µ und var(y) = a(ϕ)b^′′(θ) = a(ϕ)V(µ). Die Varianz von y ist also ein Produkt zweier Funktionen:

V(µ) h¨angt ausschließlich vom Erwartungswertµab unda(ϕ) ist vonµunabh¨angig.V(µ) nennt man Varianzfunktion, w¨ahrend ϕ als Dispersionsparameter bezeichnet wird.

Die Funktion b(θ) wirdKumulantenfunktion genannt.

(15)

2.1 Maximum Likelihood Sch¨ atzung

Liegt einen-elementige Zufallsstichprobe y₁, . . . , y_n aus der Exponentialfamilie vor, so ist der Maximum Likelihood Sch¨atzer von µdie L¨osung der Score-Funktion

∑n i=1

∂logf(y_i, θ)

∂µ =

∑n i=1

∂logf(y_i, θ)

∂θ

∂µ =

∑n i=1

y_i−b^′(θ) a(ϕ)

∂θ

∂µ. Mit b^′(θ) =µ und wegen

∂µ

∂θ = ∂b^′(θ)

∂θ =b^′′(θ) =V(µ) vereinfacht sich die obige Score-Funktion zu

∑n i=1

∂logf(y_i, θ)

∂µ =

∑n i=1

y_i −µ a(ϕ)V(µ) =

∑n i=1

y_i−µ

var(y). (2.4)

Diese recht simple Form resultiert bei der Exponentialfamilie nur bez¨uglich der Ableitung nachµ. Sie entspricht der Ableitung der Fehlerquadratsumme beim Linearen Modell mit var(y_i) = σ².

Generell könnten wir annehmen, dass beobachtungsspezifischea_i(·) vorliegen, es aber nur einen globalen Dispersionsparameter ϕ gibt. Falls speziell ai(ϕ) = ai ·ϕ mit bekannten Gewichten a_i, hängt der Maximum Likelihood Schätzer ˆµ nicht mehr von ϕ ab.

2.2 Mitglieder der Exponentialfamilie

Wir werden nun einige wichtige Mitglieder dieser Verteilungsfamilie kennen lernen. Dabei wird eine Parametrisierung verwendet, die den Erwartungswert immer durchµbezeichnet.

Die Varianzen sind dadurch oft proportional zu Potenzen von µ.

• Die Normalverteilungy∼N(µ, σ²):

f(y, µ, σ²) = 1

√2πσ² exp (

−(y−µ)² 2σ²

)

= exp

(yµ−µ²/2 σ² − y²

2σ² − 1

2log(2πσ²) )

, y ∈R. Setzen wir nunθ =µund ϕ =σ², so f¨uhrt dies zur Exponentialfamilie mit

a(ϕ) =ϕ , b(θ) = θ²/2, c(y, ϕ) =−y² 2ϕ − 1

2log(2πϕ), und mittels (2.3) zu

E(y) = b^′(θ) =θ =µ

var(y) = a(ϕ)b^′′(θ) = ϕ·1 =σ² κk(y) = 0 f¨urk >2.

(16)

2.2. MITGLIEDER DER EXPONENTIALFAMILIE 13

• Die Poissonverteilungy ∼P(µ):

f(y, µ) = µ^y

y!e⁻^µ = exp (ylogµ−µ−logy!) , y= 0,1,2, . . . . Mit θ= logµund festem ϕ = 1 f¨uhrt dies zur Exponentialfamilie mit

a(ϕ) = ϕ , b(θ) = exp(θ), c(y, ϕ) =−logy!, und mittels (2.3) zu den Kumulanten

E(y) = b^′(θ) = exp(θ) = µ var(y) = a(ϕ)b^′′(θ) = exp(θ) =µ

κ_k(y) = exp(θ) =µ f¨urk > 2.

• Die Gammaverteilungy ∼G(a, λ):

f(y, a, λ) = exp(−λy)λ^ay^a⁻¹ 1

Γ(a), a, λ, y >0.

Damit gilt E(y) = a/λ und V ar(y) = a/λ². Die Reparametrisierung µ= ν/λ mit ν =a liefert E(y) = µund V ar(y) = µ²/ν. Die entsprechende Dichtefunktion lautet damit

f(y, µ, ν) = exp (

−ν µy

) (ν µ

)ν

y^ν⁻¹ 1 Γ(ν)

= exp (

−ν

µy+νlogν−νlogµ+ (ν−1) logy−log Γ(ν) )

= exp



y (−¹_µ)

+ log_µ¹

1/ν +νlogν+ (ν−1) logy−log Γ(ν)



 , µ, ν, y >0.

Mit θ=−1/µ und ϕ= 1/ν f¨uhrt dies zur Exponentialfamilie mit a(ϕ) =ϕ , b(θ) = −log(−θ), c(y, ϕ) = 1

ϕlog 1 ϕ +

(1 ϕ −1

)

logy−log Γ (1

ϕ )

und mittels (2.3) zu den Kumulanten

E(y) = b^′(θ) =−1 θ =µ var(y) = a(ϕ)b^′′(θ) = ϕ1

θ² = 1 νµ² κ_k(y) = (k−1)!ν

(µ ν

)k

f¨urk > 2.

(17)

• Die Invers-Gaussverteilung y∼IG(µ, σ²):

f(y, µ, σ²) = 1

√2πσ²y³ exp (

− 1 2σ²y

(y−µ µ

)2)

= exp (

−y²−2yµ+µ² 2σ²yµ² − 1

2log(

2πσ²y³))

= exp



y

(−_2µ¹²) + ¹_µ

σ² − 1

2σ²y − 1 2log(

2πσ²y³)

, y >0.

Mit θ=−_2µ¹², (µ= (−2θ)⁻^1/2) und ϕ=σ² ergibt dies eine Exponentialfamilie mit a(ϕ) = ϕ, b(θ) = −(−2θ)^1/2, c(y, ϕ) = −1

2 ( 1

ϕy + log(

2πϕy³)) und mittels (2.3) zu den Kumulanten

E(y) = b^′(θ) = (−2θ)⁻^1/2 =µ,

var(y) = a(ϕ)b^′′(θ) =ϕ(−2θ)⁻^3/2 =σ²µ³, κ3(y) = 3σ⁴µ⁵, κ4(y) = 15σ⁶µ⁷.

• Die standardisierte Binomialverteilung my∼B(m, π):

f(y, m, π) = P(Y =y) = P(mY =my) = (m

my )

π^my(1−π)^m⁻^my

= exp (

log (m

my )

+mylogπ+m(1−y) log(1−π) )

= exp

(ylog₁₋^π_π −log₁₋¹_π

1/m + log

(m my

))

, y= 0, 1 m, 2

m, . . . ,1.

Mit θ= log_1−π^π , (π=e^θ/(1 +e^θ)) und ϕ= 1/m ist dies eine Exponentialfamilie mit a(ϕ) = ϕ, b(θ) = log 1

1−π = log(1 + exp(θ)), c(y, ϕ) = log (1

ϕ y ϕ

) , und mittels (2.3) zu den Kumulanten

E(y) = b^′(θ) = exp(θ)

1 + exp(θ) =π , var(y) = a(ϕ)b^′′(θ) =ϕ exp(θ)

(1 + exp(θ))² = 1

mπ(1−π), κ₃(y) = 1

m²(1−2π)π(1−π), κ₄(y) = 1

m³(1−6π(1−π))π(1−π). Das m-fache von y ist eine binomialverteilte Gr¨oße, y selbst ein relativer Anteil.

(18)

2.3. DIE QUASI-LIKELIHOOD FUNKTION 15

2.3 Die Quasi-Likelihood Funktion

Betrachtet man die Score-Funktion (2.4) zur Exponentialfamilie, so erkennt man, dass der Maximum-Likelihood Schätzer ˆµ nur von der zugrundeliegenden Varianzannahme abhängt. In diesem Abschnitt wird nun untersucht, welche Eigenschaften ein Schätzer für µaufweist, falls die Score-Funktion auch für Varianzannahmen verwendet wird, die keinen Mitgliedern aus der Exponentialfamilie gehören. Generell spricht man dann von einer Quasi-Score-Funktion. Ohne Verlust der Allgemeinheit wollen wir dazu annehmen, dass die Dispersion gegeben ist durcha(ϕ) =ϕ.

Definition 2.2 Für eine Zufallsvariable y mit E(y) = µ und var(y) = ϕV(µ) (mit a(ϕ) = ϕ) und bekannter Varianzfunktion V(·) ist die Quasi-Likelihood-Funktion q(y, µ) (eigentlich Log-Quasi-Likelihood-Funktion) definiert über die Beziehung

∂q(y, µ)

∂µ = y−µ

ϕV(µ), (2.5)

oder ¨aquivalent dazu durch q(y, µ) =

∫ _µ y−t

ϕV(t)dt+Funktion in y (und ϕ). (2.6) Die Ableitung ∂q/∂µ wird als Quasi-Score-Funktion bezeichnet. Verglichen mit (2.1) und (2.2) hat sie folgende Eigenschaften mit der Score-Funktion gemeinsam

E

(∂q(y, µ)

∂µ )

= 0, (2.7)

var

(∂q(y, µ)

∂µ )

= var(y)

ϕ²V²(µ) = 1

ϕV(µ) =−E

(∂²q(y, µ)

∂µ² )

. (2.8)

Satz 2.2 Wedderburn (1974): F¨ur eine Beobachtung y mit E(y) = µ und var(y) = ϕV(µ) hat die Log-Likelihood Funktion l(y, µ) = logf(y, µ) die Eigenschaft

∂l(y, µ)

∂µ = y−µ ϕV(µ),

dann und nur dann, wenn die Dichte bzw. Wahrscheinlichkeitsfunktion von yin der Form exp

(yθ−b(θ)

ϕ +c(y, ϕ) )

geschrieben werden kann, wobei θ eine Funktion von µ und ϕ unabh¨angig von µist.

(19)

⇒: Integration bez¨uglich µ liefert l(y, µ) =

∫ ∂l(y, µ)

∂µ dµ=

∫ y−µ ϕV(µ)dµ

= y

ϕ

∫ 1 V(µ)dµ

| {z }

θ

−1 ϕ

∫ µ V(µ)dµ

| {z }

b(θ)

= yθ−b(θ)

ϕ +c(y, ϕ).

⇐: Mit (2.3) folgt f¨ur die Kumulanten der einparametrigen Exponentialfamilie (a(ϕ) =ϕ) E(Y) =µ=b^′(θ) und var(Y) = ϕV(µ) =ϕb^′′(θ). Es gilt daher

dµ

dθ = db^′(θ)

dθ =b^′′(θ) = V(µ).

Da aber l(y, µ) = (yθ−b(θ))/ϕ+c(y, ϕ) und θ eine Funktion von µ ist, folgt

∂l(y, µ)

∂µ = y

ϕ dθ

dµ −b^′(θ) ϕ

dθ dµ

= y−µ ϕV(µ).

Mit dieser Konstruktionsidee wird f¨ur einige Varianzfunktionen der assoziierte Parameter θ hergeleitet, sowie die Quasi-Likelihood-Funktionen bestimmt.

• V(µ) = 1, ϕ =σ², y, µ∈R, (d.h. y∼N(µ, σ²)):

θ =

∫

dµ=µ, q(y, µ) =

∫ _µ y−t

σ² dt+ Funktion in y =−(y−µ)² 2σ² .

• V(µ) = µ, 0< µ, 0≤y, (d.h.y ∼P(µ)):

θ =

∫ 1

µdµ= logµ, q(y, µ) =

∫ µ

y−t

t dt=ylogµ−µ.

• V(µ) = µ², 0< µ, 0≤y, (d.h.y ∼G(µ,1)):

θ =

∫ 1

µ²dµ=−1 µ, q(y, µ) =

∫ µ

y−t

t² dt= y

µ−logµ.

(20)

2.3. DIE QUASI-LIKELIHOOD FUNKTION 17

• V(µ) = µ³, 0< µ, 0≤y, (d.h.y ∼IG(µ,1)):

θ =

∫ 1

µ³dµ=− 1 2µ², q(y, µ) =

∫ _µ y−t

t² dt =− y 2µ² + 1

µ.

• V(µ) = µ^k, 0< µ, 0≤y, k ≥3:

θ =

∫ 1

µ^kdµ=− 1 (k−1)µ^k−1, q(y, µ) =

∫ _µ y−t

t^k dt= 1 µ^k

( µ²

k−2 − yµ k−1

) .

• V(µ) = µ(1−µ), 0 < µ <1, 0≤y≤1, (d.h.mY ∼B(m, µ)):

θ =

∫ 1

µ(1−µ)dµ= log µ 1−µ, q(y, µ) =

∫ µ y−t

t(1−t)dt=ylog µ

1−µ + log(1−µ).

• V(µ) = µ²(1−µ)², 0< µ <1, 0≤y≤1:

θ =

∫ 1

µ²(1−µ)²dµ= 2 log µ

1−µ− 1

µ+ 1 1−µ, q(y, µ) =

∫ _µ

y−t

t²(1−t)²dt= (2y−1) log µ 1−µ− y

µ − 1−y 1−µ.

• V(µ) = µ+µ²/k, 0< µ, 0≤y, 0< k, (d.h. y∼N egB(k, µ)):

θ =

∫ 1

µ+µ²/kdµ= log µ k+µ, q(y, µ) =

∫ µ

y−t

t+t²/kdt=ylog µ

k+µ+klog 1 k+µ.

Während die ersten vier (Normal-, Poisson-, Gamma- und Inverse Gaußverteilung) und das sechste Beispiel (standardisierte Binomialverteilung) zu bereits bekannten Mitgliedern der Exponentialfamilie führen, stellen das fünfte sowie das siebente (speziell für Model- le für Prozentsätze) und achte Beispiel (Negativ-Binomialverteilung) neue (nicht in der Exponentialfamilie inkludierte) Varianzfunktionen dar. Hängt die Varianzfunktion von einem k ab, so muss diese Größe beim Quasi-Likelihood-Ansatz als fest betrachtet werden.

Es besteht (noch) keine M¨oglichkeit,k simultan mit µzu sch¨atzen.

(21)

2.3.1 Quasi-Dichten

Natürlich ist durch die Spezifikation einer Erwartungswert/Varianz-Beziehung auch eine Dichtefunktion spezifizierbar. Aus der (Log)-Quasi-Likelihood Funktion folgt mit der Normalisierungsfunktion

ω(µ) =

∫

R

exp(q(y, µ))dy alsQuasi-Dichte (siehe dazu Nelder & Lee (1992))

fq(y, µ) = exp(q(y, µ))

ω(µ) . (2.9)

ω(µ) ist ungleich 1, wenn die Varianz ϕV(µ) zu keiner Verteilung mit Dichte oder Wahr- scheinlichkeitsfunktion aus der Exponentialfamilie geh¨ort. Andererseits ist ω(µ) = 1, ∀µ, falls zur Varianz eine Exponentialfamilie existiert.

Zur Quasi-Dichte (2.9) korrespondiert nun die Log-Likelihood Funktion l_q(y, µ) = log(f_q(y, µ)) =q(y, µ)−log(ω(µ))

und

∂lq(y, µ)

∂µ = ∂q(y, µ)

∂µ −∂log(ω(µ))

∂µ .

Dieser Score unterscheidet sich vom Quasi-Score genau um

∂log(ω(µ))

∂µ = 1

ω(µ)

∂ω(µ)

∂µ = 1 ω(µ)

∫ ∂exp(q(y, µ))

∂µ dy

= 1

ω(µ)

∫ ∂q(y, µ)

∂µ exp(q(y, µ))dy=

∫ y−µ ϕV(µ)

exp(q(y, µ)) ω(µ) dy

=

∫ y−µ

ϕV(µ)f_q(y, µ)dy= E_q

(y−µ ϕV(µ)

)

= µ_q−µ ϕV(µ). Hierbei bezeichnet

µ_q=

∫

yf_q(y, µ)dy

den Quasi-Mean vony. Falls µ_q−µverglichen mity−µsehr klein ist, bedeutet dies, dass der Maximum-Quasi-Likelihood Schätzer sehr nahe dem Maximum-Likelihood Schätzer bezüglich der Quasi-Verteilung ist.

(22)

Kapitel 3

Das Generalisierte Lineare Modell

Unter Annahme der Existenz von E(y_i) und var(y_i) wird in der Klasse der Generalisierten Linearen Modelle (GLM) eine Parametrisierung der Form

stochastische Komponente: y_i ^ind∼ Exponentialfamilie(θ_i), E(y_i) =µ_i =µ(θ_i) systematische Komponente: η_i =x^⊤_i β

Linkfunktion: g(µ_i) =η_i

betrachtet, wobei der Zufallsvektor y = (y₁, . . . , y_n)^⊤ aus unabhängigen Komponenten y_i aufgebaut ist mit E(y_i) = µ_i und var(y_i) = ϕ_iV(µ_i). Die Dispersionsparameter ϕ_i sind gerade die Größen ai(ϕ) von zuvor. Es bezeichnet im weiteren xi = (xi1, . . . , xip)^⊤ den Vektor von bekannten erklärenden Variablen, zusammengefasst zu einer Designma- trix X = (x₁, . . . , x_n)^⊤, β = (β₁, . . . , β_p)^⊤ den Vektor mit den unbekannten Parametern, η= (η1, . . . , ηn)^⊤ den Vektor mit den Linearen Prädiktoren undg(·) eine bekannte Link- funktion.

Die wesentlichen Unterschiede zum herk¨ommlichen Linearen Modell sind:

• Es besteht keine allgemeine Additivit¨at bez¨uglich nicht-beobachtbarer Fehlerterme ϵ_i wie beim Linearen Modell,

• Eine Abh¨angigkeit der Varianzstruktur auch vom Erwartungswert ist m¨oglich,

• Eine Funktion des Erwartungswertes wird linear modelliert. Dies ist keinesfalls zu verwechseln mit einer einfachen Transformation der Responsevariablen.

Unser Hauptinteresse liegt nun in der Konstruktion eines Schätzers für den Parametervek- tor β, sowie an einem Maß f¨ur die Güte der Modellanpassung. Beides ist für Maximum- Likelihood-Schätzer besonders einfach und stellt nur eine Verallgemeinerung der Resultate für Lineare Modelle dar.

19

(23)

3.1 Maximum-Likelihood-Sch¨ atzungen

Falls y₁, . . . , y_n unabh¨angige Responses sind und die y_i aus derselben Exponentialfami- lie stammen mit Parameter (θ_i, ϕ_i), wobei der Vektor θ = (θ₁, . . . , θ_n)^⊤ die unbekannten Parameter beschreibt welche gesch¨atzt werden sollen, und ϕ = (ϕ₁, . . . , ϕ_n)^⊤ aus (vor- erst) bekannten (nuisance) Komponenten besteht, so ist die Log-Likelihood Funktion der Stichprobe gegeben durch

l(y, θ) =

∑n i=1

(y_iθ_i−b(θ_i)

ϕ_i +c(y_i, ϕ_i) )

.

Unter der recht allgemeinen Annahme µ=µ(β) folgt aus (2.4) die Score-Gleichung

∂l(y, θ(β))

∂β_j =

∑n i=1

y_i−µ_i ϕ_iV(µ_i)

∂µ_i

∂β_j = 0, j = 1, . . . , p.

Mit der Deﬁnition des linearen Pr¨adiktors gilt beim GLM

∂µ

∂β = ∂µ

∂η

∂β = ∂µ

∂g(µ)x= x g^′(µ) und deshalb

∂l(y, θ(β))

∂β_j =

∑n i=1

y_i−µ_i ϕ_iV(µ_i)

x_ij

g^′(µ_i) = 0, j = 1, . . . , p. (3.1) Den speziellen Link g(µ) = θ nennt man die kanonische Linkfunktion. Hierbei wird der Parameter θ direkt durch den linearen Pr¨adiktor η modelliert. In diesem Fall ist g(·) die Inverse von b^′(·) und wegen µ=b^′(θ) folgt

g^′(µ) = ∂g(µ)

∂µ = ∂θ

∂µ = 1

b^′′(θ) = 1 V(µ).

Die Score-Gleichung (3.1) vereinfacht sich f¨ur eine kanonische Linkfunktion zu

∂l(y, θ(β))

∂β_j =

∑n i=1

y_i−µ_i

ϕ_i xij = 0, j = 1, . . . , p. (3.2) F¨urϕ_i =ϕ(identische Dispersion) gilt hier bei Modellen mit einem Intercept (x_i1 = 1, ∀i)

∑n i=1

y_i =

∑n i=1

µ_i.

Der Maximum-Likelihood-Schätzer ˆβ ist also für den allgemeinen Fall als Lösung des Gleichungssystems (3.1) oder im kanonischen Fall durch (3.2) definiert. Beide Systeme

(24)

3.1. MAXIMUM-LIKELIHOOD-SCH ÄTZUNGEN 21 können nur iterativ gelöst werden. Einzige Ausnahme bildet das lineare Modell, in dem µ linear in β ist. Für alle anderen Situationen ist µ=g⁻¹(x^⊤β).

Die Newton-Raphson Methode liefert die Iterationsvorschrift β^(t+1) =β^(t)+

(

−∂²l(y, θ(β))

∂β∂β^⊤ )₋1

∂l(y, θ(β))

∂β , t= 0,1, . . . , (3.3) wobei beide Ableitungen der rechten Seite von (3.3) an der Stelleβ^(t) betrachtet werden.

In Matrixnotation folgt f¨ur den Scorevektor

∂l(y, θ(β))

∂β =X^⊤DW(y−µ), mit D=diag(d_i) und W =diag(w_i), wobei