11 Optimierung, Logit-, Probit-, Poisson-Regression in:

Andreas Behr, Ulrich Pötter

Einführung in die Statistik mit R, page 180 - 197

2. Edition 2010, ISBN print: 978-3-8006-3599-3, ISBN online: 978-3-8006-4878-8, https://doi.org/10.15358/9783800648788_180

Series: Vahlens Kurzlehrbücher

Bibliographic information
N O R B - Z In Rahmen statistischer Analysen treten Optimierungsprobleme insbesondere im Zusammenhang mit der Maximum-Likelihood-Schätzung auf. Im Folgenden wird der Newton-Raphson-Algorithmus dargestellt und an einem einfachen Beispiel veranschaulicht. In der Funktion optim() sind mehrere wichtige Optimierungsverfahren implementiert. Zusätzlich besprechen wir in diesem Kapitel die Logit- und Probit-Regression mit einer binären abhängigen Variablen und die Poisson-Regression für Zählvariable. 11.1 Numerische Optimierung 11.1.1 Der Newton-Raphson-Algorithmus 11.1.2 Der Befehl optim() 11.1.3 Der Befehl maxLik() 11.2 Verallgemeinerte Lineare Modelle 11.3 Logit- und Probit-Regression 11.3.1 Maximum-Likelihood-Schätzung 11.3.2 Logit- und Probit-Regression mit R 11.3.3 Probit- und Logit-Regression mit glm() 11.4 Poisson-Regression 11.4.1 Poisson-Regression mit R 11.4.2 Poisson-Regression mit glm() 11.5 Übungsaufgaben . Numerische Optimierung Wir betrachten zunächst ein einfaches Beispiel. Ausgehend von Beobachtungen soll eine nichtlineare Produktionsfunktion geschätzt werden. Die Optimierung, Logit-, Probit-, Poisson-Regression Funktion habe folgende, in den Parametern nichtlineare Form: Y = f (X, e) = θ Xθ + e Der Störterm e wird dabei als standardnormalverteilt angenommen. Wir erzeugen entsprechende Daten und stellen die Log-Likelihood graphisch dar. > set.seed(123) > theta < c(1,0.8) > n < 1000 > x < 1:n/1000 > e < rnorm(n) > y < theta[1]∗x^theta[2]+e > loglik < function(theta)sum(dnorm(y theta[1]∗x^theta[2],log=T)) > t0 < seq(0.5,1.5,0.05) > t1 < seq(0,3,0.05) > t01 < expand.grid(t0,t1) > lik.t01 < matrix(apply(t01,1,function(z) + loglik(z)),length(t0),length(t1)) > persp(t0,t1,lik.t01,phi=40,theta=40, + xlab="theta0",ylab="theta1",zlab="log L") > contour(t0,t1,lik.t01,levels=c(seq( 1480, 1420,20), 1415, + 1412, 1410.5),xlab=expression(theta[0]), + ylab=expression(theta[1]),bty="l") (a) θ θ (b) Abbildung 11.1: Log-Likelihoodfunktion (a) und derenHöhenlinien (b) Betrachten wir die Graphik der Höhenlinien, dann kann man sich ein numerischesOptimierungsverfahren folgendermaßen veranschaulichen:Durch einen gewählten Startpunkt (θ , θ ) geht genau eineHöhenlinie. Zu dieser . Numerische Optimierung Höhenlinie wird die Tangente an diesem Punkt betrachtet. Läuft man in Richtung der Senkrechten zur Tangente, so bewegt man sich auf das Maximum zu. Läuftman in diese Richtung bis zum nächsten lokalenMaximum, so erreicht man wieder eine Tangente einer Höhenlinie, an der erneut senkrecht zu dieser Tangente in Richtung des Maximums weitergelaufen werden muss. Mit diesem Vorgehen kann man beliebig nahe an das Maximum gelangen und bricht das Verfahren einem gewählten Abbruchkriterium zufolge irgendwann ab. Diese intuitive Vorstellung ist jedoch noch kein implementierbares Verfahren, da noch nicht festgelegt ist, wie das lokale Maximum in Richtung der Senkrechten zur Tangente bestimmt werden soll. Wir können aber festhalten: Da die Senkrechte zur Tangente an einem Punkt einer Höhenlinie gerade der Gradient der Funktion in diesem Punkt ist, sind die Ableitungen der Log-Likelihoodfunktion log L(.) nach den gesuchten Parametern zu bilden. . . Der Newton-Raphson-Algorithmus Der Newton-Raphson-Algorithmus benutzt zusätzlich zu dem Gradienten die Krümmung der Funktion, um zu einem numerischen Verfahren zu gelangen. Er beruht auf der Taylor-Reihen-Approximation . Ordnung. Die Log- Likelihoodfunktion log L(.) an der Stelle θt = (θ , θ )t lässt sich approximieren durch die um die Stelle θ entwickelte Funktion: log L (θ) log L θ + ∂ log L (θ) ∂θ t θ θ θ + θ θ t ∂ log L (θ) ∂θ∂θt θ θ θ Nullsetzung der Ableitung ergibt: ∂ log L (θ) ∂θ θ + ∂ log L (θ) ∂θ∂θ θ θ θ = Nach Umstellung resultiert die Iterationsvorschrift θ = θ ∂ log L (θ) ∂θ∂θt θ ∂ log L (θ) ∂θ θ Der Algorithmus konvergiert in der Regel sehr schnell, jedoch kann eine nicht positiv definite und damit nicht invertierbare Hesse-Matrix (die ausgewertete Matrix der zweiten Ableitungen) resultieren. Zudem sind die erste und zweite Ableitung der Zielfunktion nach dem Parametervektor vorzugeben. Die analytische Bestimmung der Ableitungen kann aufwendig sein. Zudem gibt es statistische Optimierungsprobleme, bei denen entweder die zweite oder sogar die erste Ableitung nicht existieren. Optimierung, Logit-, Probit-, Poisson-Regression Führen wir das oben begonnene Beispiel fort, dann benötigen wir zur Implementierung des Newton-Raphson-Algorithmus die Ableitungen der Log-Likelihoodfunktion nach dem Parametervektor θ: ∂ log L ∂θ ∂ log L ∂θ = xθ (y xθ θ ) θ log(x) ∂ log L ∂θ ∂ log L ∂θ ∂θ ∂ log L ∂θ ∂θ ∂ log L ∂θ = xθ xθ log(x)(y xθ θ ) log(x)(y xθ θ ) θ log (x)(y xθ θ ) Neben der Log-Likelihoodfunktion benötigen wir noch zwei Funktionen, die den Spaltenvektor der . Ableitungen (z) bzw. die Matrix der . Ableitungen (h) berechnen. > z < function(theta) matrix(c( + 2∗sum(x^theta[2]∗(y x^theta[2]∗theta[1])), + 2∗sum(x^theta[2]∗theta[1]∗log(x)∗(y x^theta[2]∗theta[1])) + )) > h < function(theta){ + m11 < 2∗sum(x^(2∗theta[2])) + m22 < 2∗sum(x^theta[2]∗theta[1]∗log(x)^2∗ + (y 2∗x^theta[2]∗theta[1])) + m12 < 2∗sum(x^theta[2]∗log(x)∗(y 2∗x^theta[2]∗theta[1])) + matrix(c(m11,m12,m12,m22),2,2)} Ermitteln wir nun die ersten drei Iterationen: > b0 < matrix(c(1.5,0.5)) > b.1 < b0 solve(h(b0))%∗%z(b0);t(b.1) [,1] [,2] [1,] 0.587 0.225 > b.2 < b.1 solve(h(b.1))%∗%z(b.1);t(b.2) [,1] [,2] [1,] 0.895 0.590 > b.3 < b.2 solve(h(b.2))%∗%z(b.2);t(b.3) [,1] [,2] [1,] 1.004 0.754 Der Newton-Raphson Algorithmus soll nun in einer Schleife berechnet werden, bis ein Abbruchkriterium erfüllt ist: . Numerische Optimierung > b0 < matrix(c(1.5,0.5)) > crit < 1;i < 0 > L0 < loglik(b0);L0 > while (crit>0.00001&i<20){ + if (i==0) cat(Šiteration=Š,i,Šb=Š,b0,ŠL=Š,L0,Šcrit=Š,crit,Š\nŠ) + i < i+1 + b1 < b0 solve(h(b0))%∗%z(b0) + crit < sum(abs(b0 b1)) + L0 < loglik(b1) + cat(Šiteration=Š,i,Šb=Š,b1,ŠL=Š,L0,Šcrit=Š,crit,Š\nŠ) + b0 < b1 } Für unser Beispiel erhalten wir folgenden Output: iteration= 0 b = 1.500 0.500 L= 1506 crit= 1 iteration= 1 b = 0.587 0.225 L= 1433 crit= 1.19 iteration= 2 b = 0.895 0.590 L= 1412 crit= 0.674 .... iteration= 6 b = 1.020 0.782 L= 1410 crit= 1.32e 06 Nach der sechsten Iteration ergibt sich für die sechste Nachkommastelle keine Veränderung mehr, so dass entsprechend des gesetzten Kriteriums die Iteration abgebrochen wird. . . Der Befehl optim() Mit optim()steht in R ein Befehl zur numerischen Optimierung zur Verfügung. optim() minimiert eine zu übergebende Zielfunktion. Bei einem gesuchten Parametervektor ist als erstes Argument ein Vektor mit Startwerten zu übergeben. Das zweite Argument ist der Name einer Zielfunktion. Der Befehl optim() stellt verschiedene Algorithmen zur Verfügung, die mit dem method Argument ausgewählt werden können. – Nelder-Mead: Dieser Algorithmus nutzt lediglich die Funktionswerte. Er ist damit leicht anwendbar, weil keine Ableitungen vorzugeben sind, jedoch relativ langsam. – BFGS: Dieser Algorithmus nach Broyden, Fletcher, Goldfarb und Shanno verwendet auch den Gradienten. – CG:Dieser Algorithmus beruht auf der konjugiertenGradientenmethode und gilt als weniger robust, ist jedoch für größere Optimierungsprobleme geeignet. – L-BFGS-B: Dieser Algorithmus entspricht dem BFGS-Algorithmus, jedoch können hier Ungleichheitsrestriktionen angegeben werden. – SANN: Dieser Algorithmus (Simulated Annealing) gehört zur Klasse der stochastischen globalen Optimierungsmethoden. Er benutzt lediglich die Funktionswerte und ist relativ langsam. Optimierung, Logit-, Probit-, Poisson-Regression Das bereits eingeführte nicht-lineare Optimierungsproblem lösen wir mit Hilfe der verschiedenen in optim() implementierten Algorithmen. Bei dem einfachen Beispiel konvergieren die Algorithmen alle zur gleichen Lösung. Unterschiede sind erst in den Nachkommastellen der geschätzten Parameter zu finden. Da die Voreinstellung von optim() dieMinimierung der Zielfunktion ist, schreiben wir die Log-Likelihoodfunktion mit negativem Vorzeichen (mloglik). > mloglik < function(theta){ + sum(dnorm(y theta[1]∗x^theta[2],log=T))} > b0 < c(1.5,0.5) > #Nelder Mead > Ąt1 < optim(b0,mloglik,method="Nelder Mead");t(Ąt1$par) [,1] [,2] [1,] 1.02 0.782 > #BFGS, ohne Ableitung > Ąt2 < optim(b0,mloglik,method="BFGS");t(Ąt2$par) [,1] [,2] [1,] 1.02 0.782 > #BFGS, mit Ableitung > Ąt3 < optim(b0,mloglik,z,method="BFGS");t(Ąt3$par) [,1] [,2] [1,] 1.02 0.782 > #CG, ohne Ableitung > Ąt4 < optim(b0,mloglik,method="CG");t(Ąt4$par) [,1] [,2] [1,] 1.02 0.782 > #CG, mit Ableitung > Ąt5 < optim(b0,mloglik,z,method="CG");t(Ąt5$par) [,1] [,2] [1,] 1.02 0.782 > #SANN > set.seed(123) > Ąt6 < optim(b0,mloglik,method="SANN");t(Ąt6$par) [,1] [,2] [1,] 1.01 0.779 . . Der Befehl maxLik() Das Paket maxLik stellt den Befehl maxLik() zur Verfügung, der die Arbeit mit Maximum-Likelihood-Schätzern etwas vereinfacht. So werden zusätzliche Informationen wie der Wert der Zielfunktion, die Standardfehler, t-Werte und p-Werte der geschätzten Parameter zur Verfügung gestellt. Für unser Beispiel erhält man den folgenden Output: > library(maxLik) . Verallgemeinerte Lineare Modelle > ml < maxLik(loglik,start=b0) > summary(ml) Maximum Likelihood estimation Newton Raphson maximisation, 6 iterations Return code 1: gradient close to zero Log Likelihood: 1410 2 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] 1.0163 0.0719 14.13 < 2e 16 ∗∗∗ [2,] 0.7822 0.1289 6.07 1.3e 09 ∗∗∗ . Verallgemeinerte Lineare Modelle In diesem Kapitel betrachten wir Regressionsmodelle für binäre Variable und für Zählvariable. Für binäre Variable, d.h. für Variable, die lediglich zwei Ausprägungen annehmen können, verwenden wir Logit- und Probitmodelle. Für Zählvariable, d.h. Variable, die nur ganzzahlige positive Werte (einschließlich der ) annehmen, verwenden wir ein Poisson-Modell. Diese Modelle lassen sich unter die Klasse der verallgemeinerten linearen Modelle subsumieren. Diese Modellklasse verallgemeinert das einfache Regressionsmodell, behält aber viele Eigenschaften der linearen Regression bei. Verallgemeinerte lineare Modelle haben als wesentliche Bestandteile neben einer Verteilungsannahme einen linearen Prädiktor η, der eine Linearkombination der erklärenden Variablen x ist und eine Linkfunktion g(µx), die den Erwartungswert µx vonY, also µx := E(Y x), in den linearen Prädiktor η transformiert: g(µx) = η := b + b x + . . . + bmxm Die Linkfunktion g(.) ist dabei eine strikt monotone (invertierbare) und stetig differenzierbare Funktion. Wir setzen noch h(η) := g (η). Offenkundig resultiert das lineare Regressionsmodell aus der Wahl der Identität als Linkfunktion. Die in diesem Kapitel behandelten Probit-/Logit- und Poisson-Regressionen können als Spezialfälle in diesem Modellrahmen mit unterschiedlichen Verteilungsannahmen und Linkfunktionen betrachtet werden. . Logit- und Probit-Regression Logit- und Probit-Regressionen sind spezielle Regressionsmodelle für eine dichotome abhängige Variable Y , d.h. eine Variable, die lediglich zwei Ausprägungen annehmen kann. Die Ausprägungen kodieren wir mit und . Die Verteilung der zugehörigen ZufallsvariablenY ist durch Pr(Y = x) = E(Y x) = µx gegeben. Optimierung, Logit-, Probit-, Poisson-Regression Eine einfache Möglichkeit der Modellierung besteht nun darin, den Erwartungswert direkt als Linearkombination der erklärenden Variablen zu betrachten E(Y x) = b + b x + . . . + bmxm =: xb In der Notation der verallgemeinerten Modelle würde dies der Wahl der Identität als Linkfunktion entsprechen. Dieses lineare Wahrscheinlichkeitsmodell hat jedoch Nachteile, insbesondere den, dass die Regressionswerte nicht auf das Intervall [ , ] beschränkt sind. Dies ist ungünstig, da in diesem Fall Wahrscheinlichkeiten kleiner als oder größer als resultieren könnten. Dieser Nachteil kann vermieden werden, wenn nicht die Linearkombination xb sondern eine Transformation von xb auf das Intervall [ , ] verwandt wird. Da wir eine strikt monotone Transformation in das Intervall [ , ] benutzen möchten, bietet sich die Klasse der strikt monotonen und differenzierbaren Verteilungsfunktionen an, deren Quantilsfunktion dann als Linkfunktion dienen kann. Prinzipiell können alle möglichen Verteilungsfunktionen gewählt werden, üblich ist die Wahl der Normalverteilung (Probit) Pr(Y = x, b) = xb ϕ(t) dt = Φ(xb) =: hΦ(xb) oder der logistischen Verteilung (Logit) Pr(Y = x, b) = exb + exb =: hΛ(xb) Betrachten wir die marginalen Effekte in diesem nichtlinearen Modell: ∂E(Y x, b) ∂x = ∂h(xb) ∂x = ∂h(η) ∂η ∂η ∂x =: h (xb)b Hier soll h (u) für ∂h(η)/∂η η=u stehen. Offenkundig hängen die geschätzten Wahrscheinlichkeiten E(Y x, b) = h(xb) nun nichtlinear von x ab. Der Einfluss einer partiellen Änderung von x auf h(xb) ist proportional zumWert der Dichtefunktion (genauer: der Ableitung der inversen Linkfunktion nach η) am linearen Prädiktor xb. Für das Probitmodell ergibt sich ∂E(Y x, b) ∂x = ∂Φ(xb) ∂x = ϕ(xb)b und für das Logitmodell ∂E(Y x, b) ∂x = ∂hΛ(xb) ∂x = hΛ(xb) ( hΛ(xb)) b D.h. der Einfluss einer Veränderung von x ist am höchsten, wenn die Dichte am höchsten ist. . Logit- und Probit-Regression . . Maximum-Likelihood-Schätzung Laut Modell ist jede Beobachtung ( oder ) das Ergebnis eines Bernoulli- Experiments mit Erfolgswahrscheinlichkeit Pr(Y = x) = h(xb) Da die Zufallsvariablen stochastisch unabhängig sind, ergibt sich die Wahrscheinlichkeit für das Auftreten der Werte als Produkt der Einzelwahrscheinlichkeiten Pr(Y = y ,Y = y , . . . ,Yn = yn X) = n i= Pr(Yi = yi xi) Wir betrachten nun die Beobachtungen als gegeben und den Parametervektor b als variabel und interpretieren die Wahrscheinlichkeit für die Realisationen als Likelihood des Parametervektors. Umsortieren der Wahrscheinlichkeiten nach den Ausgängen und ergibt L(b Y ,X) = i yi= [ h(xib)] i yi= h(xib) Wegen x = und x = x kann die Likelihood in einem Produkt geschrieben werden: L(b Y ,X) = n i= [ h(xib)] yi n i= h(xib) yi = n i= h(xib) yi [ h(xib)] yi Logarithmierung führt zur Log-Likelihood log L(b) = n i= yi log h(xib) + ( yi) log [ h(xib)] Die Ableitung nach dem gesuchten Parametervektor b, die Scorefunktion, ist U(b) = n i= xti yi hi hi ( yi) hi hi = n i= xti hi hi( hi) yi hi wobei hi abkürzend für h(xib) steht, hi für ∂h(η)/∂η xib. Bis auf den Faktor hi/(hi( hi)) entspricht das der Schätzgleichung der linearen Regression. Bisher haben wir eine allgemeine inverse Linkfunktion h(.) verwandt. Nun betrachten wir die Verteilungsfunktion hΛ(.) der logistischen Verteilung. Die Scorefunktion sowie die Fisher-Information und die beobachtete Information lauten dann: U(b) = n i= xti yi hΛ(xib) I(b) = Î(b) = n i= hΛ(xib) [ hΛ(xib)] x t ixi Optimierung, Logit-, Probit-, Poisson-Regression Fisher-Information und beobachtete Fisher-Information sind gleich, da die beobachtete Fisher-Information offenbar nicht mehr von den Zufallsvariablen Yi abhängt. Wählen wir die Verteilungsfunktion der Normalverteilung Φ(.), dann lautet die Scorefunktion U(b) = n i= xti ϕ(xib) Φ(xib)( Φ(xib)) yi Φ(xib) Vergleicht man die Scorefunktion von Logit- und Probitmodell, dann sieht man, dass die Scorefunktion im Logitmodell sich von der Schätzgleichung im linearen Regressionsmodell nur durch die Transformation des linearen Prädiktors in der Definition der Residuen yi hΛ(xib) unterscheidet. Dagegen hat die Scorefunktion im Probitmodell einen zusätzlichen Term. In allen verallgemeinerten linearen Modellen lässt sich immer eine Linkfunktion finden, die die entsprechende Scorefunktion auf die Form der Scorefunktion einer linearen Regression reduziert. Diese Linkfunktion nennt man kanonische Linkfunktion. Für die beobachtete Fisher-Information ergibt sich Î(b) = n i= xti ∂ ∂b ϕi Φi( Φi) (yi Φi) + x t i ϕi Φi( Φi) ∂Φi ∂b Bildet man den Erwartungswert, dann ist der erste Term , denn E(Y x) = Φ(xb). Daher ist die erwartete Fisher-Information: I(b) = n i= ϕ(xib) Φ(xib)( Φ(xib)) xtixi . . Logit- und Probit-Regression mit R Wir erzeugen uns zunächst einen Datensatz entsprechend den Modellannahmen des Probitmodells: > set.seed(1) > x < 1:10/10 0.6;n < length(x) > X < cbind(1,x) > b < matrix(1:2) > eta < X%∗%b > y < rbinom(n,1,pnorm(eta)) Der lineare Prädiktor η geht als Argument in die Verteilungsfunktion der (Standard-) Normalverteilung ein. Die so ermittelten Wahrscheinlichkeiten liegen den Realisationen der unabhängigen Binomialexperimente zugrunde, die mit Hilfe des Befehls rbinom() erzeugt werden. . Logit- und Probit-Regression Die (negative) Log-Likelihood können wir mitHilfe der implementierten Normalverteilung als Funktion schreiben, die für einen übergebenen Parametervektor den Wert der Log-Likelihood zurückgibt: > l < function(b) sum(y∗log(pnorm(X%∗%b))+ + (1 y)∗log(1 pnorm(X%∗%b))) MitHilfe des Befehls optim() finden wir mittels numerischer Optimierung das Maximum der Log-Likelihood. Da der Befehl optim() als Voreinstellung minimiert, wurde einfach das Vorzeichen der Log-Likelihood verändert. Statt die Log-Likelihood zu maximieren, minimieren wir also die negative Log- Likelihood. Dies führt zu demselben Ergebnis. Wir benennen das Ergebnis der Optimierung und können so die interessierenden Unterobjekte (Parameter, Hesse-Matrix, usw.) später abrufen. Um die geschätzten Standardfehler der Parameter ausgehend von der Hesse-Matrix ermitteln zu können, muss deren Berechnung im Aufruf von optim() mit der Option hessian=T angefordert werden. Um die Standardfehler der Parameter zu berechnen, greifen wir die Diagonalelemente der Inversen der Hesse-Matrix heraus und berechnen die Quadratwurzeln dieser Varianzen: > o < optim(1:2,l,hessian=T) > o$par [1] 0.3469283 1.3702499 > o$hessian [,1] [,2] [1,] 5.8829034 0.4719711 [2,] 0.4719711 0.4685582 > sqrt(diag(solve(o$hessian))) [1] 0.4300332 1.5237581 Den geschätzten partiellen Effekt der erklärenden Variablen x auf den Erwartungswert von Y ermitteln wir entsprechend der oben gefundenen Ableitung: > as.vector(dnorm(X%∗%o$par)∗o$par[2]) [1] 0.5162656 0.5357003 0.5455271 0.5452009 0.5347398 [6] 0.5147237 0.4862411 0.4507906 0.4101510 0.3662339 Je nach demWert der Dichtefunktion an der betrachteten Stelle x ergibt sich ein anderer, hier zunächst ansteigender dann fallender, Einfluss einer Veränderung von x auf den Erwartungswert vonY. Auf die nach dem Probitmodell generierten Daten wenden wir zur Illustration nun auch das Logitmodell an. Die Vorgehensweise ist ganz analog, nur dass wir anstelle der Verteilungsfunktion der Normalverteilung (pnorm()) nun die Verteilungsfunktion der logistischen Verteilung (plogis()) in der Definition der Log-Likelihood verwenden: > l < function(b) sum(y∗log(plogis(X%∗%b))+ + (1 y)∗log(1 plogis(X%∗%b))) Optimierung, Logit-, Probit-, Poisson-Regression > o < optim(1:2,l,hessian=T) > o$par [1] 0.5489294 2.1309884 > o$hessian [,1] [,2] [1,] 2.2054632 0.1851227 [2,] 0.1851227 0.1840546 > sqrt(diag(solve(o$hessian))) [1] 0.7037254 2.4360143 . . Probit- und Logit-Regression mit glm() Der Befehl glm() („generalized linear models“) erweitert den lm() Befehl für lineare Regressionen und erlaubt die Schätzung beliebiger verallgemeinerter linearer Modelle, behält aber gleichzeitig die Argumente zur Angabe von Formeln und die Optionen zur Behandlung von fehlenden Werten, die Erstellung von Diagnostiken, Graphiken und Zusammenfassungen bei. Im Aufruf des Befehls muss der lineare Prädiktor in der gleichen Notation wie bei dem Befehl lm() mit einer Formel spezifiziert werden (z.B. y~x1+x2). Daneben muss die Verteilungsfamilie und die Linkfunktion benannt werden. In unserem Beispiel der dichotomen abhängigen Variable, deren Realisationen aus Bernoulli-Experimenten stammen, muss dementsprechend die Option family=binomial gewählt werden. Im Falle der Probit-Schätzung, also derWahl der Verteilungsfunktion der Normalverteilung, muss zusätzlich die Linkfunktion genannt werden: family=binomial(link="probit"), da als Voreinstellung immer die kanonische Linkfunktion gewählt wird. Der Aufruf des Befehls führt zur Berechnung eines glm() Objektes, das neben den Regressionsergebnissen weitere Unterobjekte und die Regression charakterisierende Maßzahlen enthält. Um auf das Ergebnisobjekt zugreifen zu können, benennen wir es: > reg < glm(y~x,family=binomial(link="probit"));reg Call: glm(formula=y~x,family=binomial(link="probit")) Coeicients: (Intercept) x 0.347 1.370 Degrees of Freedom: 9 Total (i.e. Null); 8 Residual Null Deviance: 13.46 Residual Deviance: 12.6 AIC: 16.6 Der Aufruf des Ergebnisobjektes führt zu einer knappen Übersicht der Regressionsergebnisse. Mit dem Befehl summary(reg) lassen sich mehr Details der Rechenergebnisse anzeigen: > summary(reg) Call: . Logit- und Probit-Regression glm(formula=y~x,family=binomial(link="probit")) Deviance Residuals: Min 1Q Median 3Q Max 1.5217 1.1598 0.6764 0.9761 1.3166 Coeicients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.3470 0.4262 0.814 0.416 x 1.3703 1.4724 0.931 0.352 Null deviance: 13.460 on 9 degrees of freedom Residual deviance: 12.601 on 8 degrees of freedom AIC: 16.601 Number of Fisher Scoring iterations: 4 Die Devianz ist die zweifache Differenz der Log-Likelihood des geschätzten Modells und der Log-Likelihood, bei der alle Beobachtungen ihren eigenen Parameter haben, also dem Modell mit perfekter Anpassung (Logarithmus des Likelihoodquotienten zwischen den beiden Modellen). In einemModell für binäre Variable vereinfacht sich die Devianz zu (log L(Y) log L(b̂)) = log L(b̂) Für jede Beobachtung resultiert eine Residual-Devianz, deren Verteilung mit fünf Maßzahlen charakterisiert wird. Zudem wird das Akaike-Informationskriterium (AIC) ausgewiesen, das den zweifachen negativen Wert der Log-Likelihood ergänzt um das Zweifache der Parameterzahl angibt. Für Modellvergleiche gilt, dass das Modell mit dem geringeren Wert des Informationskriteriums vorzuziehen ist. Die geschätzten Parameterwerte und geschätzten Standardabweichungen entsprechen den Ergebnissen, die wir auch bei der direkten numerischen Optimierung erhalten haben. Nur können mit dem glm() Befehl Regressionen sehr einfach symbolisch spezifiziert werden. Zudem werden alle Diagnostiken bereitgestellt, die für den Befehl lm() schon demonstriert wurden. Und die Behandlung fehlender Werte erfolgt genau so wie in der linearen Regression. Um eine Logit-Regression zu berechnen, ist im glm() Befehl lediglich die Linkfunktion anders zu wählen. Da aber die Logit-Regression einer kanonischen Linkfunktion entspricht, braucht die Linkfunktion gar nicht explizit angegeben werden (link=ŠlogitŠ als Option geht natürlich auch): > reg < glm(y~x,family=binomial) > summary(reg) Call: glm(formula=y~x,family=binomial) Deviance Residuals: Min 1Q Median 3Q Max 1.5133 1.1642 0.6965 0.9789 1.3088 Optimierung, Logit-, Probit-, Poisson-Regression Coeicients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.5489 0.7037 0.780 0.435 x 2.1314 2.4361 0.875 0.382 Null deviance: 13.460 on 9 degrees of freedom Residual deviance: 12.631 on 8 degrees of freedom AIC: 16.631 Number of Fisher Scoring iterations: 4 Auch hier entsprechen die Parameterwerte und die geschätzten Standardfehler der Parameter den Ergebnissen der Optimierung unserer selbst geschriebenen Zielfunktion (Log-Likelihood) mitHilfe des Befehls optim(). . Poisson-Regression Im vorherigen Abschnitt haben wir uns mit einer dichotomen abhängigen Variablen beschäftigt. Nun wollen wir den Fall einer Zählvariablen betrachten. D.h. die Zufallsvariable nimmt nur ganzzahlige nicht-negative Werte an. Dabei gehen wir davon aus, dass die Zufallsvariable der Poisson-Verteilung folgt: Pr(Y = y) = e λ λy y! Dann ist der Erwartungswert geradeE(Y) =: λ. Zusätzlich betrachten wir Kovariable x, deren Ausprägungen den Erwartungswert E(Y x) =: λx beeinflussen. Da λx positiv sein muss, ist eine naheliegende Linkfunktion h(xb) = λx = e xb bzw. g(λx) = log(λx) = xb Die Linkfunktion ist also der Logarithmus. Das Modell kann dann geschrieben werden als Pr(Y = y x) = e λx λ y x y! = e exp(xb) exp(xb)y y! Die Wahrscheinlichkeit für das Auftreten der vorliegenden Beobachtungen ergibt sich wegen der angenommenen Unabhängigkeit der Beobachtungen als n i= Pr(Yi = yi xi, b) = n i= e exp(x b i ) exp(xib) yi yi! Die Höhe der Wahrscheinlichkeit hängt von dem unbekannten Parametervektor b ab. Betrachten wir nun die Daten als fix und b als variabel, dann geben verschiedene Vektoren b den vorliegenden Daten unterschiedliche Wahrscheinlichkeiten. Diese Wahrscheinlichkeiten betrachten wir als Likelihood für den . Poisson-Regression Vektor b. Als Schätzwert für den unbekannten Vektor b soll derjenige Vektor gewählt werden, der die höchste Likelihood hat, d.h. den beobachteten Daten die höchste Wahrscheinlichkeit des Auftretens gibt. Wir bilden die Log-Likelihoodfunktion log L(b Y ,X) = n i= exp(xib) + n i= yixib n i= log(yi!) Die Scorefunktion ist dann: U(b) = n i= xti yi exp(xib) Daher ist der Logarithmus die kanonische Linkfunktion für die Poisson- Verteilung. Betrachten wir die zweiten Ableitungen, dann erhalten wir die Hesse-Matrix (beobachtete Information), die wiederum mit der Fisher-Information übereinstimmt, da die beobachtete Information nicht von den ZufallsvariablenY abhängt: I(b) = Î(b) = ∂ log L(b Y ,X) ∂b∂bt = n i= exp(xib)x t ixi Da die Erwartungswerte λx unbekannt sind, schätzen wir die Varianz mit Hilfe des geschätzten Parametervektors b̂ V̂(b̂) = n i=i exp(xib̂)x t ixi Betrachten wir nun die Veränderung des Erwartungswertes vonY bei einer Veränderung von xj. Dabei sei xj die j-te Komponente von x: ∂E(Y x)/∂xj = bj exp(xib). Damit hängt die Wirkung der Veränderung von xj vom Erwartungswert, d.h. von exp(xb) ab. . . Poisson-Regression mit R Wir konstruieren zunächst ein Zahlenbeispiel, das den Annahmen des Poisson-Regressionsmodells entspricht: > set.seed(1) > n < 10;x < 1:10/10;b < 1:2 > eta < b[1]+b[2]∗x > y < rpois(n,exp(eta)) Optimierung, Logit-, Probit-, Poisson-Regression Nun wollen wir versuchen, den Parametervektor bmitHilfe der Maximum- Likelihood-Methode zu schätzen. Dafür schreiben wir eine Funktion, die für jeden übergebenen Parametervektor den Wert der Log-Likelihoodfunktion berechnet. Eine einfache Möglichkeit bietet die Verwendung der in R implementierten Poisson-Verteilung. Der Befehl dpois() liefert für einen Vektor von Realisationen Y mit Erwartungswertvektor λ die Wahrscheinlichkeiten für diese Realisationen. Die Summe der logarithmierten Wahrscheinlichkeiten entspricht demWert der Log-Likelihood. Die Erwartungswerte werden in Abhängigkeit des übergebenen Parametervektors b (in R das Objekt b) berechnet. Zudem schreiben wir noch Funktionen für die Scorefunktion und die Hesse-Matrix (jeweils mit negativem Vorzeichen). > l < function(b) sum(dpois(y,exp(X%∗%b),log=T)) > S < function(b) as.vector(t(X)%∗%(y exp(X%∗%b))) > H < function(b) t(X)%∗%(matrix(exp(X%∗%b),n,ncol(X))∗X) Verwenden wir unsere Funktion, um die Log-Likelihood an dem Parametervektor bt = ( , ) zu berechnen, erhalten wir: > l(1:2) [1] 23.43374 MitHilfe des Befehls optim() ermitteln wir nun numerisch den Maximum- Likelihood-Schätzer. Um auch Schätzer der Standardabweichungen der Parameter zu erhalten, ermitteln wir numerisch dieHesse-Matrix amMaximum- Likelihood-Schätzer. Wir benennen das Ergebnis der Optimierung und lassen uns die ermittelten Parameter und für diese Parameter die Hesse-Matrix ausgeben: > o < optim(1:2,l,S,method="BFGS") > o$par [1] 1.103164 1.885072 > H(o$par) 98.00002 68.30001 x 68.30001 54.41686 Um die Standardfehler der Parameter zu berechnen, greifen wir die Diagonalelemente der Inversen der Hesse-Matrix heraus und berechnen die Quadratwurzeln dieser Varianzen: > sqrt(diag(solve(H(o$par)))) 0.2854242 0.3830338 Mit den Maximum-Likelihood-Schätzern können wir den Effekt einer Veränderung von x um eine Einheit für die Beobachtungen aufE(Y x) = λx ermitteln: > (o$par[1]+o$par[2]∗x)∗o$par[2] [1] 2.434892 2.790241 3.145591 3.500940 3.856290 4.211639 [7] 4.566989 4.922338 5.277688 5.633037 . Übungsaufgaben Eine Veränderung von x führt zu einer umso größeren Veränderung von λx, je größer der Wert von xb vor der betrachteten Veränderung ist. . . Poisson-Regression mit glm() Die Poisson-Regression kann ebenfalls mit der Funktion glm() berechnet werden. Als Verteilung ist nun die Poisson-Verteilung im Funktionsaufruf anzugeben. Als Linkfunktion ist der Logarithmus als Voreinstellung (die kanaonische Linkfunktion) festgelegt, so dass die Linkfunktion für unser Beispiel nicht spezifiziert werden muss. > reg < glm(y~x,family=poisson) > summary(reg) Call: glm(formula=y~x,family=poisson(link="log")) Deviance Residuals: Min 1Q Median 3Q Max 1.0519 0.8372 0.2779 0.7500 1.3224 Coeicients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.1032 0.2854 3.865 0.000111 ∗∗∗ x 1.8851 0.3830 4.921 8.59e 07 ∗∗∗ Null deviance: 33.8033 on 9 degrees of freedom Residual deviance: 7.4294 on 8 degrees of freedom AIC: 50.728 Number of Fisher Scoring iterations: 4 Die geschätzten Parameter und Standardfehler der Parameter entsprechen denen, die wir durch die Optimierung der Log-Likelihood mittels der Funktion optim() erhalten haben. . Übungsaufgaben 1) Erzeugen Sie wie in Abschnitt . . Poisson-verteilte Pseudozufallszahlen (n = ) mit Erwartungswerten exp( 1+x x^2) (kanonische Linkfunktion). Dabei sei x standardnormalverteilt. Berechnen Sie eine Poisson-Regression. Welche Bilder werden vom Befehl plot() für Ihr Regressionsergebnis erzeugt? Können Sie Splines wie im Kapitel . . benutzen? 2) Erzeugen Sie zwei Vektoren x1 und x2 mit jeweils Realisationen einer im Intervall [ , ] gleichverteilten Zufallsvariablen und einen Vektor e mit Realisationen einer normalverteilten Zufallsvariablen mit Erwartungswert und Standardabweichung . . Erzeugen Sie davon Optimierung, Logit-, Probit-, Poisson-Regression ausgehend einen Vektor y mittels der folgenden Funktion und den Parametern θ = . , θ = . , θ = . , θ = : y = θ + θ log θ xθ + ( θ ) xθ + e a) Ermitteln Sie die ersten Ableitungen der Log-Likelihoodfunktion nach den Parametern. b) Ermitteln Sie die Schätzwerte der Parameter nach der Maximum- Likelihood-Methode. 3) Simulieren Sie Weibull-verteilte Pseudozufallszahlen mit den Parametern shape=2 und scale=1.5. Schreiben Sie eine Funktion, die die Log-Likelihood zurückgibt und verwenden Sie die Funktion maxLik um ML-Schätzwerte und deren Standardabweichungen zu berechnen. 4) Erzeugen Sie einen Vektor x der Länge n = der Werte . , . , ..., . a) Erzeugen Sie den linearen Prädiktor η (Vektor) für den Parametervektor bt = , . b) Erzeugen Sie die Werte der Verteilungsfunktion der Standardnormalverteilung an den Stellen η. c) Erzeugen Sie Realisationen (set.seed( )) der binomialverteilten ZufallsvariablenY mit den Wahrscheinlichkeiten Φ(η). d) Schreiben Sie eine Funktion, die für einen übergebenen Parametervektor den Wert der (negativen) Log-Likelihoodfunktion zurückgibt. e) Schreiben Sie eine Funktion, die für einen übergebenen Parametervektor den Wert der Ableitung der (negativen) Log-Likelihoodfunktion nach dem Parametervektor b zurückgibt. f) Ermitteln Sie numerische Schätzwerte für bmit Hilfe der Funktion optim(). Verwenden Sie dabei einmal das Verfahren nach Nelder und Mead und einmal den ’BFGS’-Algorithmus. g) Vergleichen Sie die Ergebnisse mit den Ergebnissen, die Sie bei Verwendung der Funktion glm() erhalten. 5) Seien (X , . . . ,Xr) r unabhängige und identisch verteilte Zufallsvariable, wobei X binomial-verteilt mit den Parametern n und p = . sei. Schreiben Sie eine R Funktion, die die Log-Likelihoodfunktion für den Parameter n berechnet. Wie würden Sie diese Funktion maximieren? Probieren Sie Ihr Verfahren in einer Simulation aus. Setzen Sie r = und n , . . . , . Erzeugen Sie dazu binomialverteilte Pseudozufallszahlen und berechnen das Maximum der Log-Likelihood. Stellen Sie die Verteilung der Schätzergebnisse für jedes n als Boxplot dar.

Chapter Preview

References

Zusammenfassung

Vorteile

- Einführung in die statistische Analyse mit R für Wirtschafts- und Sozialwissenschaftler

- Inklusive hilfreicher Tipps wie "Ansprechende Grafiken mit R gestalten"

Zum Thema

R ist ein Statistikprogramm, das kostenlos über das Internet verbreitet wird und dessen Source Codes frei zugänglich sind.

Aufgrund dieses kostenlosen Angebots gehen immer mehr Dozenten dazu über, neben SPSS auch R zu lehren bzw. SPSS durch R zu ersetzen.

In R steht dem Nutzer die gesamte Bandbreite statistischer Verfahren zur Verfügung. Durch die eigenständige Programmierumgebung ist die Software sehr flexibel und erlaubt notwendige Modifikationen und Erweiterungen verfügbarer Prozeduren.

Zum Werk

Dieses Buch führt leicht verständlich in die statistische Analyse mit R ein. Anhand von Beispielen wird die Umsetzung der wichtigsten Methoden der Statistik, wie sie üblicherweise in den Grundkursen gelehrt werden, mit R vorgestellt.

Das Buch verfolgt entsprechend zwei Ziele:

1. Vorstellung der statistischen Methoden,

2. Benutzung des Werkzeuges R zur Analyse von Daten.

Inhalt

- Grundlagen von R

- Datenbehandlung und graphische Darstellungen mit R

- Datenbeschreibungen (deskriptive Statistik)

- Wahrscheinlichkeitsverteilungen

- Regressionsanalysen

- Optimierungsverfahren

- Simulationen mit R

Neben vielen neuen, wirtschaftsorientierten Beispielen wird nun auch in die Paneldatenanalyse und Stichprobentheorie eingeführt.

Zu den Autoren

Dr. Andreas Behr ist wissenschaftlicher Mitarbeiter am Institut für Statistik und Ökonometrie der Universität Münster.

Dr. Ulrich Pötter ist wissenschaftlicher Mitarbeiter am Institut für Statistik der Universität Bochum.

Zielgruppe

Für Studierende und Dozenten der Wirtschaftswissenschaften im Bachelor an Universitäten und Fachhochschulen.