Generalisierte Lineare Modelle: Eine Einfuhrung mit

Generalisierte Lineare Modelle:

Eine Einfuhrung mit

Herwig Friedl

Inhalt

Lineares Modell (LM): Wiederholung

Transformation auf Normalverteilung (Box-Cox) Lineare Exponentialfamilie (LEF)

Generalisiertes Lineares Modell (GLM) Quasi-Likelihood Schatzung (Uberdispersion) Logistische Regression (Binomiale Responses) Loglineare Modelle (Poisson Responses)

GLMs mit zufalligen Eekten (EM Algorithmus)

Linear Regression: Modell

Paare (x_i; y_i), i=1; : : : ; n;dim(x_i) =p;x_i fest Annahme: lineares Modell beschreibt Erwartungswert Modell: E(y_i) =_i =x^>_i ,var(y_i) =², cov(y_i; y_i⁰) =0

Linear Regression: Parameterschatzer

LSE ˆ= (X^>X) ¹X^>y

y_i ^ind Normal(x^>_i ; ²)) MLE Normalˆ (; ²(X^>X) ¹) MLE ˆ²= ¹_nSSE( ˆ) = _n¹P

i(y_i ˆ_i)²,E(ˆ²) = (1 ^p_n)² (biased)

S²= _{n p}¹ SSE( ˆ) (unbiased)

y_i ^ind Normal(x^>_i ; ²)) SSE( ˆ)=^{2 2}_{n p} ANOVA: SST=SSR( ˆ) +SSE( ˆ)

SST=P

i(y_i ¯y)²,SSR( ˆ) =P

i( ˆ_i ¯y)² R²= SSR( ˆ)

SST =1 SSE( ˆ)

SST 2(0; 1), R_adj² =1 SSE( ˆ)=(n p)

SST=(n 1) 62(0; 1),

Linear Regression: Grenzen

Probleme:

y_i 6 Normal(E(y_i); var(y_i)) E(y_i)6=x^>_i 2 R

var(y_i)6=² gleich (homoskedastisch)8i =1; : : : ; n

Losungsideen:

transformiere y_i so dass g(y_i)^ind Normal(x^>_i ; ²) verwende GLM mit y_i ^ind LEF(g ¹(x^>_i ); V(_i))

Box-Cox Transformation

Deniere fur positive Responses (y > 0)

y() =





 y 1

; falls 6=0;

log y; falls =0;

Fur ! 0strebt y()! log y, so dassy()eine stetige Funktion inist.

Annahme: es existiert ein Wert fur den

yi() Normal(i() =x^>_i (); ²()) Berechne MLE bzgl. marginaler Verteilung der Responsesy.

Box-Cox Transformation

Die Transformationsfamilie liefert dafur

f(yj; (); ²()) =

















 1

p2²()exp







_{(y 1)}

()2

2²()





y ¹; 6=0;

1

p2²()exp (log y ())² 2²()

!

y ¹; =0:

Falls 6=0und () =x^>(), dann ist f(yj; (); ²()) = 1

p2²²()exp y 1 x^>()2

2²²()

!

jjy ¹:

Box-Cox Transformation

Mit0=1+0() undj =j(),j =1; : : : ; p 1, sowie ²=²²() gilt

f(yj; (); ²()) = 1

p2²²()exp y 1 x^>()2

2²²()

! jjy ¹ f(yj; ; ²) =p 1

2²exp (y x^>)² 2²

jjy ¹:

Ist=0, so sei_j =_j(),j =0; : : : ; p 1, und²=²() f(yj0; ; ²) = p 1

2²exp

(log y x^>)² 2²

y ¹: Ware bekannt, dann ist MLE sehr einfach zu berechnen!

Box-Cox Transformation

Relevanter Teil der Log-Likelihood Funktion der Stichprobe:

6=0:

`(; ; ²jy) = n

2log ² 1 2²

Xn

i=1

y_i x^>_i 2

+n log jj+( 1) Xn

i=1

log y_i

=0:

`(0; ; ²jy) = n

2log ² 1 2²

n

X

i=1

log y_i x^>_i 2 n

X

i=1

log y_i

Box-Cox Transformation: MLE's

Fur festes losen MLE'sˆ und ˆ² (vgl. multiples LM) die Scoregleichung:

@`(; ; ²jy)

@ =













 1 ²

Xn

i=1

x_i(y_i x^>_i ) =0; 6=0;

1 ²

n

X

i=1

x_i(log y_i x^>_i ) =0; =0 ;

@`(; ; ²jy)

@² =













 n

2² + 1 2⁴

Xn

i=1

(y_i x^>_i )²=0; 6=0;

n

2² + 1 2⁴

n

X

i=1

(log yi x^>_i )²=0 ; =0 :

Box-Cox Transformation: MLE's

Wir erhalten ˆ=

((X^>X) ¹X^>y; 6=0;

(X^>X) ¹X^>log y; =0;

ˆ²= 1

nSSE( ˆ) =













 1 n

Xn

i=1

(y_i x^>_i ˆ)²; 6=0;

1 n

n

X

i=1

(log y_i x^>_i ˆ)²; =0:

y (bzw.log y) sind elementeweise gerechnet. SSE( ˆ) ist die Fehlerquadratsumme vony (bzw.log y) in ˆ fur fest.

Wegen dieser Parameterisierung istSSE( ˆ) in =0unstetig.

Box-Cox Transformation: Prole-Likelihood Schatzer

Prole (Log-) Likelihoodfunktionp`(jy) =`(;ˆ;ˆ²jy) =

=













 n

2log SSE( ˆ) +n log jj+ ( 1)

n

X

i=1

log y_i; 6=0;

n

2log SSE₀( ˆ₀) Xn

i=1

log y_i; =0:

Fur 6=0folgt p`(jy) = n

2log

n

X

i=1

(y_i x^>_i ˆ)²

² + ( 1)

n

X

i=1

log yi

= n

2logXⁿ

i=1

(y_i 1)= x^>_i ˆ()₂

+ ( 1) Xn

i=1

log y_i; und somitlim_!0p`(jy) =p`(0jy) (stetig).

Box-Cox Transformation: Prole-Likelihood Schatzer

Ignorieren wir auch P_n

i=1log y_i, so folgt p`(jy) = n

2log

n

X

i=1

(y_i 1)= x^>_i ˆ() ₂

+

n

X

i=1

log y_i

= n

2logXⁿ

i=1

r_i²+log exp

"

n1 n

Xn

i=1

log y_i

#!

= n

2log

n

X

i=1

r_i²+log cⁿ

= n

2log

n

X

i=1

r_i²+n

2log(c)²

= n

2logXⁿ

i=1

r_i c

₂

mit Residuenri zum Modell x^>_i ()fur Responses (y_i 1)=, und mit der Konstantenc=exp(¹_nP

ilog y_i).

Box-Cox Transformation: Prole-Likelihood Schatzer

LRT:H₀:=₀ gegenH₁: 6=₀. LRT Statistik deniert als

Λ(y) =

2supΘ₀L(; ; ²jy) sup2ΘL(; ; ²jy); Unter gewissen Regularitatsbedingung gilt

2 log(Λ(y)) = 2

`(<sub>0</sub>;ˆ<sub>0</sub>;ˆ<sup>2</sup><sub>0</sub>jy) `(ˆ;;ˆ ˆ²jy)

= 2

p`(0jy) p`(ˆjy) _D

! ²₁:

Wegen 2(p`(<sub>0</sub>jy) p`(ˆjy)) ²₁ beinhaltet ein(1 ) Kondenzintervall alle Werte0, fur diep`(0jy) maximal

12²_1;1 vonp`(ˆjy) entfernt ist (²_1;0:95=3:841, ²₁;0:99=6:635).

Box-Cox Transformation: Eigenschaften

Log-Transformation (=0): furlog y_i Normal(x^>_i ; ²) gilt median(log y_i) =x^>_i ;

E(log y_i) =x^>_i ; var(log y_i) =²: Die originaleny_i sind lognormalverteilt mit

median(y_i) =exp(x^>_i );

E(yi) =exp(x^>_i +²=2) =exp(x^>_i )exp(²=2); var(yi) = exp(²) 1

exp(2x^>_i +²):

Additives Modell fur Erwartungswert (und fur Median) vonlog y_i ist multiplikativ fur Median und Erwartungswert vony_i. E(y_i)ist das1 < exp(²=2)-fache des Medians und die Varianz ist nicht mehr konstant furi =1; : : : ; n.

Box-Cox Transformation: Beispiel

n=31 Black Cherry Baume. Zusammenhang zwischen

HolzvolumenV in feet³, BaumhoheH in feet und Durchmesser D in inches (1 inch = 2.54 cm, 12 inches = 1 foot).

> trees <- read.table("trees.dat", header=TRUE); attach(trees)

> plot(D, V); lines(lowess(D, V)) # curvature (wrong scale?)

> plot(H, V) # increasing variance?

8 10 12 14 16 18 20

10203040506070

D

V

65 70 75 80 85

10203040506070

H

V

Box-Cox Transformation: Beispiel

> (mod <- lm(V ~ H + D)) # still fit a linear model for volume Coefficients:

(Intercept) H D

-57.9877 0.3393 4.7082

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

> h.crit <- 2*mod$rank/length(V)

> abline(h.crit, 0) # 2 leverage points

> plot(D, residuals(mod), ylab="residuals"); abline(0, 0)

> lines(lowess(D, residuals(mod))) # sink in the middle

Box-Cox Transformation: Beispiel

0 5 10 15 20 25 30

0.050.100.150.20

Index

leverages

8 10 12 14 16 18 20

−505

D

residuals

> library(MASS)

> bc<-boxcox(V~H+D,lambda=seq(0.0,0.6,length=100),plotit=FALSE))

> ml.index <- which(bc$y == max(bc$y))

> bc$x[ml.index]

[1] 0.3090909

> boxcox(V~H+D, lambda = seq(0.0, 0.6,len = 18)) # plot it now

> # directly calculate pl(lambda|y) - doesn't work if lambda=0!

> require(MASS)

> bc.trafo <- function(y, lambda) (y^lambda - 1)/lambda

> n <- length(V)

> lambda <- seq(0.01, 0.4, len=20) # avoide lambda=0

> res <- matrix(0, nrow = length(lambda), 2)

> C <- exp(mean(log(V))) # scaling constant

> for(i in seq_along(lambda)) {

+ r <- resid(lm(bc.trafo(V, lambda[i]) ~ H + D)) + pl <- -(n/2) * log(sum((r/(C^lambda[i]))^2)) + res[i, ] <- c(lambda[i], pl)

> }

> boxcox(V~H+D, lambda = lambda) # compare with box cox

> points(res[,1], res[,2], pch=16) # add points to verify match

Box-Cox Transformation: Beispiel

0.0 0.1 0.2 0.3 0.4 0.5 0.6

2122232425

λ

log−Likelihood

95%

Box-Cox Transformation: Beispiel

Volumenmessung verhalt sich kubisch in Hohe und Durchmesser?

> plot(D, V^(1/3), ylab=expression(V^{1/3}))

> lines(lowess(D, V^(1/3))) # curvature almost removed

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

(Intercept) H D

-0.08539 0.01447 0.15152

Fur=1=3fest ist median\ (V) = ˆ³₁₌₃mit E(V¹⁼³) =₁₌₃. Eˆ(V) = ˆ³₁₌₃(1+3ˆ²₁₌₃=ˆ²₁₌₃). Vergleiche Responses mit Medianen.

> mu <- fitted(mod1)

> plot(mu^3, V) # fitted median modell

Box-Cox Transformation: Beispiel

8 10 12 14 16 18 20

2.53.03.54.0

D

V13

10 20 30 40 50 60 70 80

10203040506070

μ1 3 3

V

Box-Cox Transformation: Beispiel

Krummung kann vielleicht durch log-Transformation aller Variablen reduziert werden (Regression auflog(D) und log(H)).

Soll man jetzt jedoch auf derlog(V) Achse modellieren?

> plot(log(D), log(V)) # shows nice linear relationship

> lm(log(V) ~ log(H) + log(D)) # response log(V) or still V?

Coefficients:

(Intercept) log(H) log(D)

-6.632 1.117 1.983

> boxcox(V~log(H)+log(D), lambda=seq(-0.35,0.25,length=100))

Box-Cox Transformation: Beispiel

2.2 2.4 2.6 2.8 3.0

2.53.03.54.0

log(D)

log(V)

−0.3 −0.2 −0.1 0.0 0.1 0.2

212223242526

λ

log−Likelihood

95%

Box-Cox Transformation: Beispiel

Welches der beiden Modelle ist besser? Vergleich mittels LRT.

Bette beide Modelle ein in die Modellfamilie V Normal(0+1H+2D; ²) V = (V^V 1)=_V

H = (H^H 1)=_H D = (D^D 1)=_D

Vergleiche Prole-Likelihoodfunktion in_V =1=3,_H =_D =1 (E(V¹⁼³) =₀+₁H+₂D), mit jenem in _V =_H =_D =0 (E(log(V)) =0+1log(H) +2log(D)).

Box-Cox Transformation: Beispiel

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

> bc1$y [1] 25.33313

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

> bc2$y [1] 26.11592

LRT Statistik: 2(25:333 26:116) =1:566(nicht signikant).

Box-Cox Transformation: Beispiel

Bemerkung: Schatzer zulog(H) nahe bei Eins (ˆ1=1:117) und Schatzer zulog(D) nahe bei Zwei (ˆ₂=1:983).

Baum durch Zylinder oder Kegel beschreibbar. Volumenhd²=4 (Zylinder) oderhd²=12 (Kegel), also

log(V) =c+1 log(H) +2 log(D)

mitc =log(=4) (Zylinder) oderc =log(=12) (Kegel).

Vorsicht: D von inches auf feet konvertieren ) D=12als Pradiktor.

Box-Cox Transformation: Beispiel

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

(Intercept) log(H) log(D/12)

-1.705 1.117 1.983

Konvertierung beeinusst nur Interceptwert!

Fixiere Slopes(₁; ₂) auf (1, 2) und schatze nur Intercept₀, d.h. betrachte das Modell

E(log(V)) =₀+1 log(H) +2 log(D=12):

Den Term1 log H+2 log(D=12) nennt man oset (Pradiktor mit festem Parameter).

Box-Cox Transformation: Beispiel

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

(Intercept) -1.199

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

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

Holzvolumen kann eher durch ein Kegelvolumen als durch das eines Zylinders beschrieben werden, hat jedoch ein etwas groeres Volumen als ein Kegel.

Lineare Exponentialfamilie (LEF): Denition

Denition: Eine Zufallsvariabley sei aus einer Verteilung mit Dichte- oder Wahrscheinlichkeitsfunktion

f(yj) =exp

y b()

a() +c(y; )

fur bekannte Funktionena(), b() und c() mita()> 0.

Istfest, nennt man f(yj) einparametrige, lineare Exponentialfamilie in kanonischer Form (LEF) mit kanonischem Parameter.

Lineare Exponentialfamilie (LEF)

Bemerkung: Allgemeine Exponentialfamilie (AEF) deniert durch

f(yj) =h(y)p()exp







k

X

j=1

t_j(y)w_j()







;

mit reellen Funktionenh(y) 0und t1(y); : : : ; t_k(y) in y sowie p() 0und w₁(); : : : ; w_k() in .

Um darin LEF zu erkennen, schreiben wir AEF um zu

f(yj) =exp





 Xk

j=1

t_j(y)w_j() +log(p()) +log(h(y))





 : Setzek =1(einparametrig), t(y) =y (linear),w() = (kanonische Parametrisierung), sowielog(p()) = b() und log(h(y)) =c(y; ). Bis auf Skalierung ist dies eine LEF.

Lineare Exponentialfamilie: Momente

Bemerkung: Fur die AEF (log(p()) = b()) gilt E





k

X

j=1

@w_j()

@_l t_j(y)



= @

@_l log(p()); var





k

X

j=1

@wj()

@_l t_j(y)



= @²

@_l²log(p()) E





k

X

j=1

@²wj()

@²_l t_j(y)



:

Istk =1,t(y) =y und w() =, liefert diesE(y) =b⁰()und var(y) =b⁰⁰().

Lineare Exponentialfamilie: Score & Information

Bemerkung: Fur Scorefunktion und Informationszahl gilt:

E

@ log f(yj)

@

=0;

var

@ log f(yj)

@

=E

@ log f(yj)

@

₂

=E

@²log f(yj)

@@^>

: Hier:

E

@ log f(yj)

@

= 1

a()E y b⁰()

=0;

alsoE(y) =b⁰() , und E

@²log f(yj)

@²

+E

@ log f(yj)

@

₂

= 1

a()b⁰⁰()+ 1

a²()var(y) =0 : Alsovar(y) =a()b⁰⁰().

Lineare Exponentialfamilie: Kumulanten

SeiE(y) =b⁰() =und somit var(y) =a()b⁰⁰() =a()V(). Varianz vony ist Produkt zweier Funktionen: V() und a().

VarianzfunktionV(), Dispersionsparameter. b() nennt man Kumulantenfunktion.

Kumulantenerzeugenden FunktionK(t) =log M(t). k-te Kumulante_k gegeben durchK^(k)(t)j_t=0. Zusammenhang mit Momenten

1(y) =E(y) 2(y) =E(y )² ₃(y) =E(y )³

₄(y) =E(y )⁴ 3var²(y):

Lineare Exponentialfamilie: Kumulantenerzeugende

Fur die LEF gilt 1=

Z

Rexp

y b()

a() +c(y; )

dy

=exp

b() a()

Z

Rexp y

a() +c(y; )

dy ; also

exp b()

a()

= Z

Rexp y

a() +c(y; )

dy :

Lineare Exponentialfamilie: Kumulantenerzeugende

Wegen exp

b() a()

= Z

Rexp y

a() +c(y; )

dy folgt als Momentenerzeugende FunktionM(t)

E(e^ty) =exp

b() a()

Z

Rexp y

a()

+a()t

+c(y; )

dy

=exp

b() a()

exp



 b

+a()t a()





=exp



 b

+a()t b() a()



:

Lineare Exponentialfamilie: Kumulantenerzeugende

Wegen

E(e^ty) =exp



 b

+a()t b() a()



 folgt als Kumulantenerzeugende Funktion

K(t) =log M(t) = b

+a()t b() a() : Diek-te Kumulante vony ist somit

_k(y) =K^(k)(t)jt=0= a()^{k 1}b^(k)

+a()t

t=0=a()^{k 1}b^(k)():

Lineare Exponentialfamilie: MLE fur Erwartung

Annahme: y₁; : : : ; y_n (iid) Zufallsstichprobe aus LEF().

MLE furist Nullstelle der Scorefunktion Xn

i=1

@ log f(y_ij)

@ =

Xn

i=1

@ log f(y_ij)

@

@

@ = Xn

i=1

y_i b⁰() a()

@

@: Mitb⁰() = und wegen (Ableitung der inversen Funktion)

@

@ = @b⁰()

@ =b⁰⁰() =V() ist die Scorefunktion

Xn

i=1

@ log f(y_ij)

@ =

Xn

i=1

y_i a()V() =

Xn

i=1

y_i var(y_i):

Lineare Exponentialfamilie: MLE fur Erwartung

MLEˆ lost Xn

i=1

@ log f(y_ij)

@ =

Xn

i=1

y_i

a()V() =0 :

Entspricht der Ableitung der Fehlerquadratsumme beim LM mit var(y_i) =² (Normalverteilung).

Bekannte Losung fur gesamte LEF:

ˆ= 1 n

Xn

i=1

y_i:

Lineare Exponentialfamilie: Dispersion

Von nun an Annahme: beobachtungsspezische Funktionen a_i() hangen nur von einem globalen Dispersionsparameterab!

(sonst ist die Anzahl der Dispersionsparameter gleichn).

Beispiel: N Mittel ¯yk von Stichproben mit Umfangenn1; : : : ; nN. Nur diese Mittely¯_k = _n¹

k

P_nk

i=1y_ki sind beobachtet.

Fallsy_ki (iid) Zufallsstichprobe mit E(y_ki) = und var(y_ki) =², dannE(¯y_k) =und var(¯y_k) =²=n_k =a_k mita_k =1=n_k bekannt und unbekannter Dispersion=².

Daher werden wir uns im Folgenden ausschlielich auf den Fall a_i() =a_i mit bekannten Gewichtena_i beschranken. Unter diesem Modell hangt der MLEˆ nicht mehr von ab.

Lineare Exponentialfamilie: Mitglieder

Normalverteilung y Normal(; ²):

f(yj; ²) = p 1 2²exp

(y )² 2²

=exp

y ²=2 ²

y² 2²

1

2log(2²)

; y 2 R:

Setze= und =², so fuhrt dies zur LEF mit a=1 ; b() =²=2 ; c(y; ) = y²

2 1

2log(2); wofur gilt

E(y) =b⁰() == var(y) =b⁰⁰() = 1=²

_k(y) =0 fur k > 2 :

Lineare Exponentialfamilie: Mitglieder

Poissonverteilungy Poisson(): f(yj) = ^y

y!e =exp(y log log y!); y =0; 1; 2; : : : : Setze=log und =1, so fuhrt dies zur LEF mit

a=1 ; b() =exp(); c(y; ) = log y!; wofur gilt

E(y) =b⁰() =exp() = var(y) =b⁰⁰() =exp() =

_k(y) =exp() = fur k > 2 :

Die Dispersion ist bei der Poissonverteilung bekannt Eins und somit wirklich kein freier Parameter.

Lineare Exponentialfamilie: Mitglieder

Gammaverteilung y Gamma(a; ): f(yja; ) =exp( y)^ay^{a 1} 1

Γ(a); a; ; y > 0 : Dafur istE(y) =a=undvar(y) =a=².

Reparametrisierung==, =a gibtE(y) =,var(y) =²= und Dichtefunktion fur; ; y > 0

f(yj; ) =exp

y

y ¹ 1 Γ()

=exp

y+ log log + ( 1)log y logΓ()

=exp



 y

1

+log¹

1= + log + ( 1)log y logΓ()



:

Lineare Exponentialfamilie: Mitglieder

Dichtefunktion fur; ; y > 0 f(yj; ) =exp



 y

1

+log¹

1= + log + ( 1)log y logΓ()



:

Setze= 1=und =1=, so fuhrt dies zur LEF mit a=1 ; b() = log( ); c(y; ) = 1

log1 +

1

1

log y logΓ 1

;

wofur gilt

E(y) =b⁰() = 1 = var(y) =b⁰⁰() =1

² = 1 ² _k(y) = (k 1)!

_k

fur k > 2 : Varianz ist proportional zum Quadrat des Erwartungswertes.

Lineare Exponentialfamilie: Mitglieder

Inverse Gaussverteilung y InvGauss(; ²):

f(yj; ²) = 1

p2²y³exp 1 2²y

y

2!

=exp



 y

21²

+^{1 2}

1 2²y

1

2log 2²y³



; y > 0:

Mit= ₂¹₂,(= ( 2) ¹⁼²) und=² gibt dies LEF mit a=1; b() = ( 2)¹⁼²; c(y; ) = 1

2 1

y +log 2y³

E(y) =b⁰() = ( 2) ¹⁼²=; var(y) =b⁰⁰() =( 2) ³⁼²=²³;

3(y) =3⁴⁵; 4(y) =15⁶⁷: Varianz wachst proportional zu³.

Lineare Exponentialfamilie: Mitglieder

Standardisierte Binomialverteilung my Binomial(m; ): f(yjm; ) =Pr(Y =y) =Pr(mY =my) =

m my

^my(1 )^{m my}

=exp

log m

my

+my log +m(1 y)log(1 )

=exp y log₁ log₁ ¹ 1=m +log

m my

!

; y =0; 1 m; 2

m; : : : ; 1:

Mit=log₁ ,( =e=(1+e)) und=1ist dies LEF mit a= 1

m; b() =log 1

1 =log(1+exp()); c(y; ) =log ₁

y

:

Lineare Exponentialfamilie: Mitglieder

Standardisierte Binomialverteilung my Binomial(m; ):

E(y) =b⁰() = exp()

1+exp() = ; var(y) =a b⁰⁰() = 1

m1 exp()

(1+exp())² = 1

m(1 ); ₃(y) = 1

m²(1 2)(1 ); ₄(y) = 1

m³(1 6(1 ))(1 ):

Relative Haugkeity. Absolute Haugkeit my ist binomialverteilt.

Bemerke, dass Dispersion bekannt Eins ist undm reziprok als Gewicht eingeht.

Lineare Exponentialfamilie: Quasi-Likelihood

Scorefunktion fur bei LEF hangt nur von der Varianz ab.

Idee: verwende auch Varianz zu der gar kein LEF Mitglied existiert. Man spricht von Quasi-Scorefunktion.

Annahme: Dispersion seia =, also Gewichta sei Eins.

Denition: Fur eine Zufallsvariabley mitE(y) =und var(y) =V() mit bekannter Varianzfunktion V()ist die (Log-)Quasi-Likelihoodfunktionq(jy) deniert uber die Beziehung

@q(jy)

@ = y V(); oder aquivalent dazu durch

q(jy) =

Z y t

V(t)dt+Funktion in y (und ) :

Lineare Exponentialfamilie: Quasi-Likelihood

Die Ableitung@q=@nennt man Quasi-Scorefunktion.

Verglichen mit der Scorefunktion hat sie folgende Eigenschaften gemeinsam

E

@q(jy)

@

=0;

var

@q(jy)

@

= var(y)

²V²() = 1

V() = E

@²q(jy)

@²

:

Lineare Exponentialfamilie: Quasi-Likelihood

Satz (Wedderburn, 1974) Fur eine Beobachtungy mit E(y) =und var(y) =V() hat die Log-Likelihoodfunktion

`(jy) =log f(yj) die Eigenschaft

@`(jy)

@ = y V();

genau dann, wenn die Dichte- bzw. Wahrscheinlichkeitsfunktion vony in der Form

exp

y b()

+c(y; )

geschrieben werden kann, wobeieine Funktion von, und unabhangig von ist.

Lineare Exponentialfamilie: Quasi-Likelihood

Beweis:

): Integration bezuglich liefert

`(jy) =

Z @`(jy)

@ d=

Z y V()d

= y

Z 1 V()d

| {z }

1

Z V()d

| {z }

b()

=y b()

+c(y; ):

Lineare Exponentialfamilie: Quasi-Likelihood

Beweis:

(: Mit den Kumulanten der LEF gilt E(y) ==b⁰() und var(y) =V() =b⁰⁰(). Daher ist

d

d = db⁰()

d =b⁰⁰() =V():

Da aber`(jy) = (y b())=+c(y; )und eine Funktion vonist, folgt

@`(jy)

@ = y

d d

b⁰()

d d

= y

V():

Lineare Exponentialfamilie: QL Modelle

V() =1,=²,y; 2 R, (vgl. mity Normal(; ²)):

= Z

d= ; q(jy) =

Z y t

² dt+Funktion in y = (y )² 2² : V() =,0 < ,0 y, (vgl. mity Poisson()):

= Z 1

d=log ; q(jy) =

Z y t

t dt =y log :

Lineare Exponentialfamilie: QL Modelle

V() =², 0 < ,0 y, (vgl. mit y Gamma(; 1)):

= Z 1

²d= 1 ; q(jy) =

Z y t

t² dt = y

log : V() =³, 0 < ,0 y, (vgl. mit y InvGauss(; 1)):

= Z 1

³d= 1 2²; q(jy) =

Z y t

t² dt = y 2² + 1

:

Lineare Exponentialfamilie: QL Modelle

V() =^k,0 < ,0 y,k 3:

= Z 1

^kd= 1 (k 1)^{k 1}; q(jy) =

Z y t

t^k dt= 1 ^k

² k 2

y k 1

:

V() =(1 ),0 < < 1,0 y 1, (vgl. mit my Binomial(m; )):

=

Z 1

(1 )d=log 1 ; q(jy) =

Z y t

t(1 t)dt=y log

1 +log(1 ):

Lineare Exponentialfamilie: QL Modelle

V() =²(1 )², 0 < < 1,0 y 1:

=

Z 1

²(1 )²d=2 log 1

1 + 1

1 ; q(jy) =

Z y t

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

y

1 y 1 : V() =+²=k,0 < ,0 y,0 < k, (vgl. mit

y NegBinomial(k; )):

=

Z 1

+²=kd=log k+; q(jy) =

Z y t

t+t²=kdt=y log

k+ +k log 1 k+:

Lineare Exponentialfamilie: Quasi Dichte

Durch Spezikation der Erwartungswert/Varianz-Beziehung ist auch Dichtefunktion spezizierbar.

Aus der (Log)-QL Funktion folgt mit

!() = Z

Rexp(q(jy))dy als Quasi-Dichte (vgl. Nelder & Lee, 1992)

fq(yj) = exp(q(jy))

!() :

Normalisierungsfunktion!()6=1, wenn Varianz V() zu keinem LEF Mitglied gehort. Andererseits ist!() =1,8, falls zur Varianz eine LEF existiert.

Lineare Exponentialfamilie: Quasi Dichte

Zur Quasi-Dichtef_q(yj) korrespondiert Log-Likelihoodfunktion

`_q(jy) =log(f_q(yj)) =q(jy) log(!()) und Scorefunktion

@`_q(jy)

@ = @q(jy)

@

@ log(!())

@ :

Dieser Score unterscheidet sich vom Quasi-Score um

@ log(!())

@ = 1

!()

@!()

@ = 1

!()

Z @ exp(q(jy))

@ dy

= 1

!()

Z @q(jy)

@ exp(q(jy))dy

=

Z y

V()f_q(yj)dy =E_q

y V()

= _q

V() :

Lineare Exponentialfamilie: Quasi Dichte

Unterschied liegt in

@ log(!())

@ = _q

V() : Hierbei bezeichnet

q = Z

yfq(yj)dy

den Quasi-Mean vony. Falls _q verglichen mit y klein ist, ist der Maximum-Quasi-Likelihood Schatzer nahe dem Maximum-Likelihood Schatzer bezuglich der Quasi-Verteilung.

Generalisiertes Lineares Modell (GLM)

Klasse der Generalisierten Linearen Modelle (GLM):

Parametrisierung der Form

stochastische Komponente: y_i ^ind LEF(_i); E(y_i) =_i =(_i) systematische Komponente: _i =x^>_i

Linkfunktion: g(_i) =_i:

Responsevektory= (y₁; : : : ; y_n)^> aus unabhangigeny_i mit E(yi) =i und var(yi) =aiV(i). Dispersion ist Produkt ai. x_i = (x_i0; x_i1; : : : ; x_{i;p 1})^> sei p 1Vektor bekannter Pradiktoren zusammengefasst zurn p Designmatrix X= (x₁; : : : ; x_n)^>, = (₀; ₁; : : : ; _{p 1})^> sei p 1Vektor unbekannter Parameter, = (₁; : : : ; _n)^> sei n 1 Vektor mit den Linearen Pradiktoren, g()sei eine bekannte (monoton und 2x stetig dierenzierbare) Linkfunktion.

Generalisiertes Lineares Modell (GLM)

Wesentlichen Unterschiede zum herkommlichen LM sind:

keine allgemeine Additivitat bzgl. nicht-beobachtbarer Fehlerterme_i wie beim LM,

Varianzstruktur kann auch vom Erwartungswert abhangen, Funktion des Erwartungswertes wird linear modelliert (nicht

zu verwechseln mit einer Transformation der Response).

Interessiert in der Konstruktion eines Schatzers von und an einem Ma fur die Gute der Anpassung. Beides ist fur MLE sehr einfach und stellt nur Generalisierung der Resultate beim LM dar.

Generalisiertes Lineares Modell: MLE

Seieny₁; : : : ; y_n unabhangige Responses aus derselben LEF mit Parameter(_i; ), so ist die Log-Likelihood der Stichprobe

`(jy) = Xn

i=1

y_ii b(_i)

a_i +c(y_i; )

: Unter der Annahme=()folgt die Scorefunktion

@`(()jy)

@_j =

n

X

i=1

y_{i i} a_iV(_i)

@_i

@_j ; j =0; 1; : : : ; p 1 : Mit der Denition des linearen Pradiktors=x^> gilt beim GLM

@

@ = @

@

@

@ = @

@g()x= x g⁰() und deshalb

@`(()jy)

@_j =

n

X

i=1

y_{i i} a_iV(_i)

x_ij

g⁰(_i); j =0; 1; : : : ; p 1 :

Generalisiertes Lineares Modell: MLE

Dab⁰(_i) =_i und wegeng(_i) =x^>_i gilt _i =b^{0 1}(_i) =b^{0 1}(g ¹(x^>_i )):

Den Linkg() =b^{0 1}(), also g(_i) =_i =x^>_i , nennt man kanonisch. Damit wird direkt durch modelliert und es folgt

g⁰() = @g()

@ = @

@ = 1

b⁰⁰() = 1 V(): Scorefunktion vereinfacht sich fur kanonische Links zu

@`(()jy)

@_j =

n

X

i=1

y_{i i}

a_i x_ij; j =0; 1; : : : ; p : Fallsa_i =1gilt bei Modellen mit Intercept (x_i0=1; 8i)

n

X

i=1

y_i =

n

X

i=1

ˆ_i:

Generalisiertes Lineares Modell: MLE

Sind allea_i =1und die Dispersionbekannt, so folgt bei kanonischen Links

f(yj()) =Yⁿ

i=1

exp

y_ii() b(_i())

n

Y

i=1

exp(c(y_i; ))

=exp 1

n

X

i=1

y_ix^>_i b(x^>_i )

! _n Y

i=1

exp(c(y_i; ))

=g(T(y)j)h(y);

undT(y) =X^>y ist somit suziente Statistik fur. Bemerke, dassdim(T(y)) =dim() =p gilt, und dass eine suziente Statistik ausschlielich bei einer kanonischen Linkfunktion existiert.

Generalisiertes Lineares Modell: MLE

MLEˆist Nullstelle der Scorefunktion. Gleichungssystem nur numerisch (iterativ) losbar. Einzige Ausnahme ist das LM, in dem linear in ist. Fur alle anderen Situationen ist=g ¹(x^>).

Newton-Raphson Methode liefert Iteration (t =0; 1; : : :) ^(t+1) =^(t)+

@²`(()jy)

@@^>

₁ =^(t)

@`(()jy)

@ =^(t): Ableitungen auf der rechten Seite sind in^(t) betrachtet. In Matrixnotation folgt fur den Scorevektor

@`(()jy)

@ = 1

X^>DW(y ); mitD=diag(d_i) und W=diag(w_i), wobei

d_i =g⁰(_i);

wi = aiV(i)g⁰²(i) ₁ :

Generalisiertes Lineares Modell: MLE

Als negative Hessematrix der Log-Likelihoodfunktion resultiert

@²`(()jy)

@@^> = 1 X^>

@DW

@^> (y ) DW @

@^>

X

= 1 X^>

W @DW

@^> (y )

X ;

wegen@=@^>=D ¹. Weiters istd_iw_i = a_iV(_i)g⁰(_i) ₁

@d_iw_i

@_i = a_iV⁰(_i)^@_@ⁱ

ig⁰(_i) +a_iV(_i)g⁰⁰(_i)^@_@ⁱ (a_iV(_i)g⁰(_i))² i

= V⁰(_i)g⁰(_i) +V(_i)g⁰⁰(_i) a_iV²(_i)g⁰³(_i) :

Generalisiertes Lineares Modell: MLE

Zusammenfassen von

w_i =w_i @d_iw_i

@_i (y_{i i})

zuW (mitE(W) =W) liefert Newton-Raphson Vorschrift

^(t+1) =^(t)+(X^>W^(t)X) ¹X^>D^(t)W^(t)(y ^(t)); t =0; 1; : : : :

Bemerke, dass das Produkt Scorevektor mal inverse Hessematrix hier unabhangig vom Dispersionsparameter ist.

Generalisiertes Lineares Modell: MLE

Mittels Pseudobeobachtungen (adjusted dependent variates) z=X+W ¹DW(y )

die Newton-Raphson Vorschrift umschreiben in eine Iterative (Re)Weighted Least Squares Prozedur (IWLS/IRLS)

^(t+1)= (X^>W^(t)X) ¹X^>W^(t)z^(t):

Generalisiertes Lineares Modell: MLE

Fur kanonische Links giltg⁰() =1=V() und

g⁰⁰() = V⁰()=V²() und es verschwinden die Ableitungen

@d_iw_i

@_i = V⁰(_i)=V(_i) V⁰(_i)=V(_i) a_i=V(_i) =0 ) W=W, und Pseudobeobachtungen vereinfachen sich zu

z=X+D(y ) =X+V ¹(y ); mitV=diag(V(_i)). Dies liefert als Iterationsvorschrift

^(t+1) = (X^>W^(t)X) ¹X^>W^(t)z^(t):

Generalisiertes Lineares Modell: MLE

Erinnerung an LM: LS Schatzer deniert als ˆ= (X^>X) ¹X^>y

(keine iterative Losung notwendig daW=I und z=y).

Generell wird gerne mit erwarteter Information gearbeiten. Wegen E(X^>WX) =X^>WXfolgt fur diese Fisher Scoring Technik wieder

^(t+1)= (X^>W^(t)X) ¹X^>W^(t)z^(t) mit den Pseudobeobachtungen

z=X+D(y ):

Dafur istE(z) =X und wegenvar(y) =(DWD) ¹ folgt var(z) =D var(y)D=W ¹.