Bayesianische Inferenzstatistik: Eine Einführung in die Parameterschätzung in R mit dem Bayesian-First-Aid-Paket
Christoph Jörger & Johannes Gerwien
Libraries
library(BayesianFirstAid)Einleitung
Dieses Dokument ist als theoretische und praktische Einführung in bayesianische Verfahren zur statistischen Inferenz gedacht. Ausgehend von drei typischen statistischen Fragestellungen sollen die konzeptionellen Unterschiede zwischen frequentistischer und bayesianischer Inferenzstatistik sowie das praktische Vorgehen in entsprechenden Analysen vermittelt werden.
Wir betrachten die praktische Umsetzung bayesianischer Analysen in R am Bespiel des Bayesian-First-Aid-Pakets (für eine Einführung in die Arbeit mit R siehe hier). Dieses Paket implementiert einfache bayesianische Alternativen zu gängigen frequentistischen Inferenzverfahren, konkret zum t-Test, Korrelationstest und Proportions- bzw. Chi-Quadrat-Test, anhand derer sich das Prinzip bayesianischer Parameterschätzung gut veranschaulichen lässt. Es sei allerdings darauf hingewiesen, dass die Möglichkeiten bayesianischer Inferenzstatistik weit über die Anwendungsfälle und konkreten Lösungen hinausgehen, die wir hier besprechen. Wir bieten hier also nur einen einführenden Einblick (für weiterführende Ressourcen siehe hier). Zu den hier besprochenen Verfahren sei vorab angemerkt, dass sie stark datengetrieben, also mit relativ wenigen Vorannahmen arbeiten. Dieser Umstand sollte in Überlegungen zur Eignung der Verfahren für die eigene Datenanalyse berücksichtigt werden.
Das Bayesian-First-Aid-Paket wurde von Rasmus Bååth entwickelt. Die verschiedenen Verfahren werden in mehreren seiner Blog-Posts vorgestellt (t-test, Korrelationstest und Proportionstest).
Grundlagen
Wozu Inferenz?
Wir beginnen unsere Überlegungen mit einem einfachen Beispiel: Wir nehmen an, wir haben ein Experiment mit lexikalischer Entscheidungsaufgabe durchgeführt. 100 Versuchspersonen mussten bei Präsentation einer Buchstabensequenz per Tastendruck angeben, ob es sich bei der Sequenz um ein Wort handelt oder nicht. Wir nehmen weiterhin an, in dem Experiment fand nur eine kritische Manipulation statt: Die Vorkommenshäufigkeit der echten Wörter. Diese waren entweder hochfrequent (ca. ≥ 30 pro millionen Wörter) oder niedrigfrequent (ca. ≤ 5 pmw). Wir bezeichnen die Gruppe der hochfrequenten Wörter nachfolgend mit \(HF\) (high frequency), die der niedrigfrequenten Wörter mit \(LF\) (low frequency). Unsere Hypothese ist, dass die Worterkennung in der \(LF\)-Gruppe länger dauert, was in den Reaktionszeiten der Versuchspersonen reflektiert sein sollte.
Betrachten wir die Mittelwerte unserer (simulierten) Daten aus beiden Bedingungen, zeigt sich die erwartete Tendenz: Die Reaktionszeiten bei \(LF\) sind im Mittel länger als bei \(HF\). Die Differenz der Mittelwerte beträgt ca. 43 ms.
mean(LF) #Reaktionszeiten bei niedrigfrequenten Wörtern## [1] 714.3288mean(HF) #Reaktionszeiten bei hochfrequenten Wörtern## [1] 670.8887
Wir sind nun mit dem grundlegenden Problem statistischer Inferenz konfrontiert: Wie schlussfolgern wir, ob das vorliegende Ergebnis aus unserer Stichprobe einen Effekt in der Population reflektiert? Auf unser Beispiel bezogen: Liegt die Tendenz zur längeren Reaktionszeit bei niedrigfrequenten Wörtern nun zufällig bei den 100 Versuchspersonen vor, die wir getestet haben, oder auch bei den Individuen, die wir nicht getestet haben?
Die Frage setzt natürlich voraus, dass wir uns für mehr Personen interessieren als diejenigen, deren Daten uns vorliegen. Andernfalls würde sich die Frage nach statistischer Inferenz erübrigen, da wir die Population von Interesse vollständig erfasst hätten und es folglich keine unbekannte Größe gäbe, auf die man inferieren müsste. Wann immer unsere Daten aber nur einen Teil der Population reflektieren, auf die unsere Fragestellung abzielt (was fast immer der Fall ist), müssen wir einen Weg finden, anhand unserer Daten auf diese Population zu schließen. Um diesem Problem zu begegnen, liegen verschiedene Ansätze vor. Wir werden uns hier hauptsächlich mit dem bayesianischen Ansatz beschäftigen. Um dessen Motivation besser zu verstehen, betrachten wir aber zunächst den tradtionelleren, frequentistischen Ansatz.
Frequentistische Inferenz
Der frequentistische Ansatz liefert uns das wohl verbreitetste Verfahren statistischer Inferenz, den Nullhypothesen-Signifikanztest (NHST). Um zu evaluieren, ob das Muster in unseren Daten den Zustand der Population reflektiert, stellt dieses Verfahren die Frage, ob das gegenteilige Szenario mit hinreichend geringer Irrtumswahrscheinlichkeit ausgeschlossen werden kann. Letzteres Szenario bildet die Nullhypothese \(H_0\). Bei einer gerichteten Hypothese wie der unseren sagt die \(H_0\) aus, dass der Unterschied zwischen unseren Bedingungen in der Population entweder nicht existiert, oder anders gerichtet ist, als von uns vorhergesagt (\(H_0: LF \leq HF\)). Demgegenüber steht unsere Alternativhypothese \(H_1\), nach der ein Unterschied in der von uns angenommenen Richtung vorliegt (\(H_1: LF > HF\)).
Der konkrete Test für unser Szenario ist ein einseitiger t-Test für abhängige Stichproben, also für Messwertpaare derselben Merkmalsträger (alle Versuchspersonen haben sowohl auf hochfrequente als auch auf niedrigfrequente Wörter reagiert). Der Test erhält seinen Namen durch die Teststatistik, die wir bei seiner Durchführung ermitteln. In unerem Fall ist das ein t-Wert. Dieser Wert gibt an, wie groß ein beobachteter Unterschied in uneren Daten im Verhältnis zur Streuung der Daten ist. Jede Versuchsperson hat einen individuellen Mittelwert für hoch- und niedrigfrequente Wörter und eine daraus folgende Bedingungsdifferenz \(LF - HF\). Wir bilden den Mittelwert dieser Bedingungsdifferenzen und normalisieren diesen Wert an seinem Standartfehler. Auf diesem Weg erhalten wir die Teststatistik \(t\) unserer Daten.
Die t-Statistik folgt unter Annahme der \(H_0\) einer spezifischen Wahrscheinlichkeitsverteilung, einer t-Verteilung mit Erwartungswert bei \(0\) und \(N-1\) Freiheitsgraden (in unserem Fall also 99). Wir benötigen diese Verteilung als Referenz. Konkret müssen wir wissen, wie wahrscheinlich ein Wert größer oder gleich unserem vorliegenden t-Wert innerhalb dieser Verteilung ist. Liegt diese Wahrscheinlichkeit unterhalb des festgelegten Alpha-Niveaus (die maximal zu akzeptierende Irrtumswahrscheinlichkeit, konventionell \(\alpha = 0.05\)), verwerfen wir die Nullhypothese. Führen wir den entsprechenden Test in R durch, erhalten wir folgendes Ergebnis.
t.test(LF, HF, paired= TRUE, alternative = "greater")##
## Paired t-test
##
## data: LF and HF
## t = 2.5625, df = 99, p-value = 0.005949
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
## 15.29274 Inf
## sample estimates:
## mean difference
## 43.44005
Der t-Wert unserer Stichprobe (\(t = 2.56\) ) erhält einen sehr niedrigen p-Wert in der zugrundeliegenden t-Verteilung (\(p = 0.006\)). Ein Wert größer oder gleich dem vorliegenden ist in dieser Verteilung also unwahrscheinlicher als die für uns maximal zulässige Irrtumswahrscheinlichkeit \(\alpha = 0.05\). Ein entsprechendes Ergebnis führt nach der dichotomen Inferenzlogik des NHST zur Ablehnung der Nullhypothese. Wir gehen folglich davon aus, dass die Ergebnisse unserer Stichprobe auf einen analogen Zustand in der Population hinweisen. Die Inferenzlogik des NHST ist in der unteren Abbildung anhand unseres Beispiels visualisiert.
Die gestrichelte Linie in der Abbildung tangiert den kritischen Wert der t-Verteilung, bei dem eine Wahrscheinlichkeit von 0.05 für einen entsprechenden oder einen noch extremeren Wert gilt. Die Fläche unterhalb der Kurve rechts dieses Wertes bildet also die oberen 5 Prozent der Verteilung ab, die den Ablehnungsbereich der Nullhypothese darstellen. Wie wir sehen können, fällt der t-Wert aus unserer Stichprobe (rot) in diesen Bereich.
Von Frequentistischer zu Bayesianischer Inferenz
Um die Motivation für einen bayesianischen Inferenzansatz besser zu verstehen lohnt es sich, genauer zu betrachten, auf welcher Grundlage wir uns oben gegen die Nullhypothese (und damit implizit für die Alternativhypothese) entschieden haben:
Bei Betrachtung der oben dargestellten t-Verteilung fällt auf, dass diese symmetrisch um ein Maximum bei \(0\) organisiert ist. Zur Erinnerung: \(t\) quantifiziert in unserem Szenario die standardisierte Differenz zwischen unseren Bedingungen. In der oben gezeigten Verteilung liegt der Erwartungswert für dieses Maß bei \(0\), was keinem Unterschied zwischen den Bedingungen entspricht. Ferner sind kleinere Werte bzw. Werte nahe \(0\) immer wahrscheinlicher als größere Werte. Diese Verteilung ist ein Modell für die Datengenerierung in einer Nullhypothesen-Population. Wir gehen also von einer Population aus, in der kein Effekt vorliegt, und modellieren dann, wie sich das relevante Effektmaß verteilen würde, wenn wir wiederholt Stichproben aus dieser Population ziehen würden. Die Frage, die wir mit Hilfe unseres t-Werts beantworten wollen, ist also: Wie gut sind unsere Daten mit einem Populationsmodell vereinbar, das den Effekt von Interesse verneint? Im vorliegenden Fall reflektiert unser \(t\) Daten, die mit diesem Modell kaum vereinbar sind (\(p = 0.006\)). Folglich lehnen wir die Hypothese ab, die das Modell motiviert (Nullhypothese).
Ein zentraler Aspekt dieser Inferenzlogik ist, dass wir im gesamten Entscheidungsprozess nie evaluieren, wie wahrscheinlich die eine oder die andere Hypothese (\(H_0\) vs.\(H_1\)) tatsächlich ist. Wir ermitteln lediglich die Wahrscheinlichkeit, unsere Daten zu beobachten, wenn die Nullhypothese gilt.
Als Studierende oder Forschende interessiert uns häufig jedoch der umgekehrte Fall: Wie wahrscheinlich erscheint unsere Hypothese bei Betrachtung unserer Daten? Vielleicht interessiert uns zunächst auch nur, welcher konkrete Populationszustand ausgehend von unseren Daten plausibel erscheint. Auf unser Beispiel bezogen: Welches Ausmaß hat der Unterschied zwischen unseren experimentellen Bedingungen in der Population, wenn wir anhand unserer Daten urteilen? Entsprechende Urteile erlaubt die bayesianische Inferenz, deren Grundzüge nachfolgend erläutert werden. Die frequentistische Methode soll hierbei nicht als der “schlechtere” Ansatz verstanden werden. Es soll lediglich eine Sensibilisierung dafür erfolgen, welche konkreten inferenzstatistischen Fragen eine gegebene Methode überhaupt zu beantworten versucht, und wie geeignet die Methode dementsprechend für die Evaluierung unserer eigenen empirischen Fragen ist.
Bayes-Theorem und Parameterschätzung
In den vorangehenden Abschnitten haben wir gesehen, wie wir die Wahrscheinlichkeit unserer Daten gegen ein festes Modell der Population evaluieren, genauer gesagt gegen ein Modell der Art und Weise, wie sich Daten verteilen, die aus dieser Population gewonnen werden. Die Frage, die uns nachfolgend beschäftigt, mag ähnlich wirken, ist aber eine grundlegend andere: Gegeben sei ein Parameter \(\theta\), der die Population in Hinblick auf unsere Fragestellung beschreibt. In unserem Beispiel ist das die Differenz zwischen den Reaktionszeiten in unseren Bedingungen. Für einen beliebigen Wert \(i\) des Parameters \(\theta\) stellen wir nun die Frage: Mit welcher Wahrscheinlichkeit nimmt \(\theta\) in der Population einen Wert größer oder kleiner gleich \(i\) an? In unserem Fall beispielsweise: Mit welcher Wahrscheinlichkeit gilt in der Population, dass \(\theta ≤ 0\) (\(H_0\)) und komplementär: Mit welcher Wahrscheinlichkeit gilt, dass \(\theta > 0\) (\(H_1\))? Wir wollen diese Frage zwar anhand unserer Daten beantworten, doch bezieht sich die von uns gesuchte Wahrscheinlichkeit auf den Zustand der Population selbst, nicht auf die Daten!
Aber wie gelingt uns eine solche Inferenz? Wir stellen uns hierfür vor, unsere Daten werden nicht in einem epistemischen Vakuum interpretiert. Stattdessen aktualisieren unsere Daten den Wissensstand, der in Bezug auf unsere Frage vor Betrachtung der Daten vorlag. Die Wahrscheinlichkeit, mit der wir einen bestimmten Zustand in der Population annehmen, geht dann daraus hervor, wie stark unsere Daten die Wahrscheinlichkeit verändern, mit der wir diesen Zustand vor ihrer Betrachtung angenommen haben. Verstehen wir unsere Inferenz als eine solche Form der Wissensaktualisierung, ist das Verfahren mathematisch in Form des Bayes-Theorems beschreibbar.
\[p(a|b) = \frac{p(a)\cdot p(b|a)
}{p(b)}\]
Wenn wir die Formel versprachlichen, ergibt sich die Wahrscheinlichkeit von \(a\) gegeben \(b\) aus folgendem Bruch: Im Zähler steht das Produkt aus der Wahrscheinlichkeit von \(a\) und der Wahrscheinlichkeit von \(b\) gegeben \(a\), im Nenner die Wahrscheinlichkeit von \(b\). Wir werden nachfolgend versuchen, die einzelnen Elemente der Formel auf unser inferenzstatistisches Vorhaben abzubilden.
Posterior-Verteilung
Die gesuchte Wahrscheinlichkeit \(p(a|b)\) entspricht in unserem Fall der Wahrscheinlichkeit eines bestimmten Parameterwertes \(\theta_i\) vor Hintergrund unserer Daten \(D\). In der Praxis suchen wir allerdings keine einzelne Wahrscheinlichkeit \(p(\theta_i|D)\) für einen bestimmten Wert, sondern eine stetige Wahrscheinlichkeitsverteilung \(f(\theta|D)\), die verschiedenen Werten von \(\theta\) eine Wahrscheinlichkeitsdichte zuordnet. Dieser Verteilung können wir anschließend Wahrscheinlichkeiten für spezifische Werteintervalle entnehmen, beispielsweise \(p(\theta > 0)\). Die gesuchte Wahrscheinlichkeitsverteilung wird Posterior-Verteilung genannt, da sie ausdrückt, welche Wahrscheinlichkeit wir bestimmten Zuständen in der Population nach Betrachtung der Daten zuordnen. Diese Verteilung bildet das wesentliche Instrument unserer statistischen Inferenz.
Prior-Verteilung
Die Wahrscheinlichkeit \(p(a)\) entspricht in unserem Fall der Wahrscheinlichkeit, die wir einem gegebenen Wert \(i\) unseres Parameters \(\theta\) vor Betrachtung der Daten zuordnen. Entsprechend wird hierbei von der Prior-Wahrscheinlichkeit gesprochen. Erneut arbeiten wir nicht mit diskreten Wahrscheinlichkeiten \(p(\theta_i)\), sondern mit Wahrscheinlichkeitsdichten für einzelne Parameterwerte, die aus einer Prior-Verteilung \(f(\theta)\) hervorgehen. Diese Verteilung enkodiert unsere Vorannahme über den Zustand der Population in Bezug auf den Parameter, welche bei Betrachtung unserer Daten um die Posterior-Verteilung aktualisiert wird. Aber wie gelangen wir zu einer solchen Vorannahme?
Die Vorannahme ist zunächst als formale Notwendigkeit in der Anwendung der bayesianischen Methode zu verstehen: Um unser Wissen zu aktualisieren muss es schließlich etwas geben, das aktualisiert wird. Das schließt aber nicht aus, dass wir in unserer Vorannahme ausdrücken, praktisch wenig über den Parameter zu wissen. In diesem Fall läge eine relativ uninformative Prior vor. Diese kann beispielsweise als flache Normalverteilung oder Gleichverteilung über einem relativ großen Wertebereich enkodiert werden. Wir evaluieren dann im weiteren Verlauf, wie stark diese Annahme durch unsere Daten verändert wird. Ebenso möglich ist allerdings, dass wir aufgrund theoretischen oder empirischen Vorwissens spezifischere Annahmen über den Parameter machen können. In diesem Fall läge eine relativ informative Prior vor, die ebenso in einer entsprechenden Verteilungsform enkodiert werden könnte. Hierin zeigt sich, dass die Bestimmung der Prior-Verteilung kein objektives Verfahren ist. Wichtiger als die konkrete Form der Prior-Verteilung ist, wie wir diese Begründen!
Likelihood
Die bedingte Wahrscheinlichkeit \(p(b|a)\) entspricht der Wahrscheinlichkeit, unsere Daten zu beobachten, wenn \(\theta\) in der Population den Wert \(i\) annimmt. Konzeptuell bewegen wir uns hier also wieder im frequentistischen Rahmen. Entsprechend benötigen wir zunächst ein Modell über die Verteilung der Datenpunkte, die aus einer entsprechenden Population erhoben werden. Wir bestimmen hierfür eine Wahrscheinlichkeitsverteilung, deren Form durch den Parameter \(\theta\) mitbestimmt wird, und ordnen diesem Parameter dann den aktuellen Wert \(i\) zu (Beispielsweise eine Normalverteilung mit Mittelwert bei \(\theta_i\)). Wir bestimmen dann das Produkt der Wahrscheinlichkeitsdichten all unserer Datenpunkte innerhalb dieser Verteilung. Die resultierende Größe wird als Likelihood \(L\) bezeichnet und quantifiziert in unserem Fall, wie gut unsere Daten zu einem Modell der Datengenerierung passen, das von einer Population mit \(\theta_i\) ausgeht: Je höher die individuellen Dichtewerte unserer Datenpunkte, desto höher deren Produkt. Ein hoher Wert von \(L(D|\theta_i)\) drückt also aus, dass unsere Daten gut mit einem entsprechenden Modell vereinbar sind.
Marginal Likelihood
Die Wahrscheinlichkeit \(p(b)\) im Nenner des Bayes-Theorems entspricht in unserem Fall der grundsätzlichen Wahrscheinlichkeit, unsere Daten zu beobachten. Praktisch entspricht diese Wahrscheinlichkeit der Likelihood unserer Daten \(L(D)\) über den gesamten Parameterraum \(\theta\) hinweg. Dieser Wert wird als marginal Likelihood bezeichnet und fungiert im Bayes-Theorem als normalisierende Konstante, die das Produkt im Zähler auf einen echten Wahrscheinlichkeitsraum skaliert (0-1). Die relative Form der Posterior-Verteilung ändert sich hierdurch nicht. In der Praxis ist es oft schwer, die marginal Likelihood analytisch zu bestimmen, vor allem dann, wenn die bayesianische Methode mehr als einen Parameter gleichzeitig schätzen soll (was häufig der Fall ist). Unter anderem aus diesem Grund wird die gesuchte Posterior-Verteilung häufig nicht analytisch bestimmt.
Angepasstes Bayes-Theorem
Mit Hilfe der vorangehend erläuterten Konzepte können wir das Bayes-Theorem anwenden, um unseren Parameter zu schätzen.
\[f(\theta|D) = \frac{f(\theta)\cdot L(D|\theta) }{L(D)}\]
Die Formel ordnet einem gegebenen Wert \(\theta_i\) eine Posterior-Dichte zu. Diese ergibt sich aus dem Produkt der Prior-Dichte von \(\theta_i\) und der Likelihood unserer Daten \(D\) gegeben \(\theta_i\), geteilt durch die marginal Likelihood unserer Daten. Theoretisch können wir über dieses Verfahren also eine stetige Posterior-Verteilung für \(\theta\) ableiten, die wir für unsere Inferenz nutzen können. Praktisch findet ein solches analytisches Vorgehen aber selten Anwendung: Wie oben bereits angemerkt, ist die marginal Likelihood in vielen Fällen eine analytisch kaum zu berechnende Größe, sodass für einen gegebenen Parameterwert nur dessen unstandardisierte Posterior-Dichte (Prior*Likelihood) berechnet werden kann. Mit Hilfe dieses Wertes kann die Posterior-Verteilung des Parameters aber algorithmisch ermittelt werden. Hierbei wird eine Stichprobe generiert, deren Verteilungseigenschaften durch die gesuchte Posterior-Verteilung bestimmt sind. Zu diesem Zweck werden spezielle Sampling-Algorithmen, sogenannte Markov-Chain-Monte-Carlo-Verfahren eingesetzt (s.u.).
Bayesianischer Mittelwertsvergleich
Das erste Verfahren des Bayesian-First-Aid-Pakets, das wir betrachten werden, kann unmittelbar auf unser Datenbeispiel angewendet werden: Wir ermitteln darin auf Grundlage unserer Daten die Wahrscheinlichkeit, dass sich die Reaktionszeiten zwischen unseren Bedingungen unterscheiden, mit Richtung entsprechend unserer Vorhersage. Das Verfahren kann folglich als bayesianischer Ansatz für den zuvor demonstrierten t-Test für abhängige Stichproben verstanden werden. Diese Bezeichnung sollte allerdings nicht zu wörtlich genommen werden: Streng genommen findet hier kein formaler statistischer Test statt, wie wir ihn oben betrachtet haben. Wir schätzen lediglich eine Posterior-Verteilung der Bedingungsdifferenz auf Grundlage spezifischer Modellannahmen (s.u.) und leiten hieraus Wahrscheinlichkeiten für bestimmte Zustände in der Population ab.
Parameter und Priors
Das Verfahren schätzt insgesamt drei Parameter, die die Population in Bezug auf unsere Forschungsfrage beschreiben sollen. Die drei Parameter erhalten jeweils individuelle Priors, die nachfolgend kurz erläutert werden. Weitere Details sind Kruschke (2013 https://doi.org/10.1037/a0029146) zu entnehmen, auf dessen Arbeit das Verfahren zurückgeht.
\(\mu\) (mu) entspricht in unserem Fall der mittleren Bedingungsdifferenz in der Population. Dieser Wert ist für unsere Inferenz also der zentralste. Die Prior-Verteilung für \(\mu\) ist eine Normalverteilung mit Mittelwert entsprechend unserem Stichprobenmittelwert und einer Standardabweichung entsprechend unserer Stichprobenstandardabweichung \(s * 1000\). Durch diese extrem große Standardabweichung entsteht eine extrem flache und breite Verteilung, wie die untere Abbildung (obere Grafik) zeigt. Diese Verteilung enkodiert also eine relativ uninformative Prior, da die Wahrscheinlichkeitsdichte aller Werte für \(\mu\) gegen \(0\) strebt (vergleiche die Y-Achse in der Abbildung).
\(\sigma\) (sigma) beschreibt die Streuung um \(\mu\) in der Population. Es wird also geschätzt, wie stark Individuen tendenziell vom Populationsmittelwert abweichen. Die Prior-Verteilung für \(\sigma\) ist eine Gleichverteilung mit Minimum bei unserer Stichprobenstandardabweichung \(s / 1000\) und Maximum bei \(s * 1000\). Über einem sehr großen Wertebereich erhalten daher alle möglichen Werte für \(\sigma\) a-priori dieselbe Wahrscheinlichkeit (vgl. die Grafik links unten). Auch hier handelt es sich also um eine relativ uninformative Prior.
\(\nu\) (nu) beschreibt das Gewicht der Randbereiche einer Verteilung um \(\mu\) in der Population, also die relative Häufigkeit extremer Werte. Ein niedriger Wert von \(\nu\) zeigt an, dass die Randbereiche einer Verteilung um \(\mu\) relativ viel Gewicht erhalten und Extremwerte daher häufiger auftreten. Je höher \(\nu\), desto mehr nähert sich das Gewicht der Randbereiche der Form einer Normalverteilung an. Die Prior-Verteilung für \(\nu\) ist eine Exponentialverteilung mit einer Wachstumsrate von \(\frac{1 } {29}\). Wie der unteren Abbildung (rechts unten) zu entnehmen ist, erhalten schwere Randbereiche (~ \(\nu < 20\)) damit zwar höhere Wahrscheinlichkeitsdichten als leichte Randbereiche, kumulativ ist das Verhältnis zwischen schweren und leichten Randbereichen aber ausgeglichen, sodass keine starke Vorannahme über die relative Häufigkeit extremer Werte vorliegt.
Das Populationsmodell, das sich aus der Kombination dieser drei Parameter ergibt, entspricht einer t-Verteilung: Das Zentrum der Verteilung liegt bei \(\mu\), die Streuung der Verteilung wird durch \(\sigma\) beschrieben und die Anzahl der Freiheitsgrade der Verteilung entspricht \(\nu\). Für jeden der drei Parameter werden individuelle Posterior-Verteilungen ermittelt. Allerdings erfolgt das entsprechende Verfahren nicht unabhängig für jeden Parameter. Stattdessen werden Posterior-Dichten immer für eine Kombination der Werte \(\mu_i\), \(\sigma_i\) und \(\nu_i\) bestimmt. Die Prior-Dichte im Bayes-Theorem entspricht daher ebenfalls einer Kombination, konkret dem Produkt der individuellen Priors \(p(\mu_i)\), \(p(\sigma_i)\) und \(p(\nu_i)\), die jeweils aus den oben beschriebenen Prior-Verteilungen hervorgehen.
Likelihood
Analog zur kombinierten Prior-Dichte ergibt sich auch die Likelihood im vorliegenden Fall aus einer Kombination von \(\mu_i\), \(\sigma_i\) und \(\nu_i\). Wir ermitteln also die Passung unserer Daten zu einer Verteilung, die durch diese drei Parameter bestimmt ist. Wie oben bereits erwähnt, beschreiben die vorliegenden Parameter gemeinsam eine t-Verteilung. Wir nehmen entsprechend an, dass die Daten aus der Population t-verteilt sind. Für die Bestimmung der Likelihood unserer Daten gegeben die Parameterkombination \(\mu_i\), \(\sigma_i\), \(\nu_i\) verwenden wir daher eine t-Verteilung mit Zentrum bei \(\mu_i\), einer Streuung von \(\sigma_i\) und Freiheitsgraden entsprechend \(\nu_i\). Wir ermitteln dann das Produkt der Wahrscheinlichkeitsdichten unserer Datenpunkte innerhalb dieser Verteilung. Das Verfahren ist unten visualisiert für eine zufällige Parameterkombination \(\mu = 20\), \(\sigma = 100\) und \(\nu = 10\). Wir visualisieren die Dichten einzelner Messwerte unserer Stichprobe mittels roter Linien.
Bestimmung der Posterior-Verteilung
Wie oben bereits erwähnt, suchen wir im vorliegenden Verfahren Posterior-Dichten für Kombinationen aus drei Parametern. Dieser Umstand verdeutlicht das Problem der analytischen Bestimmung von Posterior-Dichten: Der Nenner des Bayes-Theorems, die marginal Likelihood, müsste in diesem Fall die kumulative Likelihood unserer Daten über alle denkbaren Wertekombinationen der drei Parameter ausdrücken, was rechnerisch kaum zu bewerkstelligen ist.
In der Praxis wird dieses Problem mit Hilfe von Stichproben umgangen,
deren Eigenschaften durch die gesuchte Posterior-Verteilung bestimmt
werden. Entsprechende Stichproben werden aus
Markov-Chain-Monte-Carlo-Verfahren gewonnen. Es handelt
sich hierbei um Sampling-Algorithmen, deren interne Evaluationskriterien
so beschaffen sind, dass die Häufigkeit einzelner Parameterwerte in der
Stichprobe langfristig aus der Form der Posterior-Verteilung hervorgeht.
Hierzu werden ausgehend von einem gegebenen Parameterwert wiederholt
lokale Vorschläge für neue Parameterwerte gezogen. Die
Wahrscheinlichkeit, mit der ein neuer lokaler Wert vom Sampler
übernommen wird, hängt von dessen unstandardisierter Posterior-Dichte
(Zähler des Bayes-Theorems: Prior * Likelihood) relativ zum vorangehend
gezogenen Wert ab.
Für das hier betrachtete Verfahren wird die gesuchte Posterior-Verteilung durch eine Stichprobe repräsentiert, die sich aus drei Markov-Ketten mit je 10.000 Einzelwerten zusammensetzt.
Zusammenfassung der Ergebnisse
Der R-Code für unsere Parameterschätzung ist fast identisch mit dem des klassischen t-Tests, der oben gezeigt wurde. Es wird lediglich bayes vorangestellt. Wir speichern das Ergebnis des Tests als Objekt, um seine Eigenschaften besser untersuchen zu können.
results_bayes = bayes.t.test(LF, HF, paired = TRUE)
results_bayes##
## Bayesian estimation supersedes the t test (BEST) - paired samples
##
## data: LF and HF, n = 100
##
## Estimates [95% credible interval]
## mean paired difference: 42 [8.8, 77]
## sd of the paired differences: 166 [140, 193]
##
## The mean difference is more than 0 by a probability of 0.993
## and less than 0 by a probability of 0.007
Aus dem Verfahren geht hervor, dass die mittlere Differenz zwischen den Reaktionszeiten mit einer Wahrscheinlichkeit \(p > 0.993\) über \(0\) liegt. Dieser Schluss beruht darauf, dass mehr als 99.3% der Masse der Posterior-Verteilung im Wertebereich über \(0\) liegen, und analog weniger als 0.7% darunter. An dieser Stelle lohnt sich ein Vergleich mit dem klassischen t-Test: Dieser hat ergeben, dass die Wahrscheinlichkeit unserer Daten bei Annahme einer Differenz von \(0\) in der Population sehr gering ist. Wir konnten hierdurch indirekt schließen, dass eine Population ohne Unterschied unplausibel ist. Mit der bayesianischen Parameterschätzung können wir - ausgehend von unseren spezifischen Modellannahmen und Priors - diesen Schluss unmittelbar ziehen, denn die Posterior-Verteilung sagt uns, mit welcher Wahrscheinlichkeit wir einen gegebenen Zustand für unseren Parameter in der Population annehmen können. Da 99.3% der Masse der Verteilung über \(0\) liegen, können wir inferieren, dass der Parameter in der Population mit größter Wahrscheinlichkeit über \(0\) liegt.
Wir können allerdings noch präzisere Schlüsse ziehen: Die beste Schätzung der mittleren Differenz unserer Bedingungen liegt bei 42 ms. Wie wir gleich sehen werden, entspricht das dem Mittelwert und dem Median der MCMC-Stichprobe. Zudem wird uns ein 95% Credible Interval angezeigt. Das Intervall entspricht hier einem 95% Highest Density Interval (HDI). Dieses Intervall bildet den Wertebereich, in dem ausgehend vom Median 95% der Masse der Posterior-Verteilung liegen. In unserem Fall liegen diese zwischen 8.8 ms und 77 ms. Ein HDI, das wie in unserem Fall vollständig oberhalb von \(0\) verläuft, kann als Anhaltspunkt für einen Effekt in der Population bei festgelegtem Sicherheitsniveau (hier 95%) betrachtet werden.
Wir erhalten zudem das Ergebnis der Schätzung für \(\sigma\), also für die mittlere Streuung um \(\mu\) in der Population: Die beste Schätzung (höchste Posterior-Dichte) liegt hierbei bei einer Streuung von 166 ms, mit einem 95% HDI zwischen 140 und 193 ms.
Details und Visualisierung
Unser Ergebnis-Objekt enthält eine Tabelle stats. Diese enthält Kennwerte der MCMC-Stichprobe, die im Rahmen des Verfahrens generiert wurde. Wir betrachten zunächst nur die erste Zeile der Tabelle. Diese enthält die Kennwerte für den Parameter \(\mu\), also die mittlere Differenz der Reaktionszeiten (wir runden die Kennwerte auf die fünfte Nachkommastelle).
round(results_bayes$stats[1,],5)## mean sd HDI% comp HDIlo HDIup
## 42.51268 17.31342 95.00000 0.00000 8.76734 76.84134
## %>comp %<comp q2.5% q25% median q75%
## 0.99253 0.00747 8.48736 30.92815 42.45031 54.06951
## q97.5% mcmc_se Rhat n_eff
## 76.64257 0.12803 1.00043 18329.00000
Der Output ist eine ausführlichere Darstellung der Ergebnisse: Wir sehen z.B. den Mittelwert unseres Parameters in der Stichprobe (42 ms), sowie dessen Standardabweichung (17.3). Wir sehen auch erneut die Grenzwerte des 95% HDI (HDIlo = 8.8 ms, HDIup = 76.9 ms) sowie den Anteil der Stichprobenwerte oberhalb und unterhalb \(0\) (%>comp und %<comp bei comp = \(0\)). Zusätzlich sehen wir weitere Quartilswerte sowie den Median und den Standardfehler der Stichprobenwerte (mcmc_se).
Wir erhalten dieselben Informationen über die Posterior-Verteilungen der Parameter \(\sigma\) und \(\nu\), wenn wir weitere Zeilen der Tabelle unseres Ergebnis-Objekts betrachten. Die erste Zeile des Outputs zeigt uns die Kennwerte der Posterior-Verteilung von \(\mu\), die wir oben bereits betrachtet haben. Die zweite Zeile zeigt dieselben Kennwerte für die Posterior-Verteilung von \(\sigma\), die dritte Zeile analog die Kennwerte für \(\nu\).
results_bayes$stats[1:3,]## mean sd HDI% comp HDIlo HDIup %>comp
## mu_diff 42.51268 17.31342 95 0 8.767342 76.84134 0.9925338
## sigma_diff 166.26562 13.34838 95 0 140.248392 192.73474 0.9999667
## nu 41.48160 30.47850 95 0 3.298512 102.50710 0.9999667
## %<comp q2.5% q25% median q75% q97.5%
## mu_diff 7.466169e-03 8.487362 30.92815 42.45031 54.06951 76.64257
## sigma_diff 3.333111e-05 141.291296 157.24179 165.70828 174.75011 193.95546
## nu 3.333111e-05 7.289781 19.48341 33.23225 54.83011 121.04700
## mcmc_se Rhat n_eff
## mu_diff 0.1280343 1.000427 18329
## sigma_diff 0.1149750 1.000517 13484
## nu 0.3726719 1.000369 6911
Das Paket bietet auch die Möglichkeit, die Ergebnisse der Parameterschätzung visuell zusammenzufassen. Wir verwenden hierfür den Befehl plot() mit unserem Ergebnis-Objekt als Argument. Der Output zeigt uns ein Histogramm der Posterior-Verteilungen für \(\mu\) (links oben) und \(\sigma\) (links unten). Diesen sind dieselben Informationen zu entnehmen, die wir oben bereits betrachtet haben: Die jeweiligen Mediane sowie die jeweiligen 95% HDIs. Aus beiden Parametern kann zudem eine Posterior-Verteilung für die Stärke des Effekts in der Population abgeleitet werden (Cohen’s \(d\) = \(\mu/\sigma\)). Die Verteilung ist im rechten oberen Histogramm zu sehen und impliziert einen kleinen (Median = 0.26) aber robusten Effekt, da das 95% HDI vollständig über \(0\) liegt bzw. weniger als 1% der Verteilungsmasse unter \(0\) liegen.
plot(results_bayes)
MCMC-Samples
Wir können auch eine noch tiefere Detailebene betrachten. Das Ergebnis-Objekt enthält die drei Markov-Ketten, die zur Parameterschätzung erstellt wurden. Diese enthalten jeweils 10.000 Einzelwerte für die drei geschätzten Parameter sowie für die Effektstärke, die aus einer gegebenen Schätzung hervorgeht. Wir betrachten unten beispielhaft die ersten 10 Zeilen aus der ersten Markov-Kette.
results_bayes$mcmc_samples[1][1:10,]## [[1]]
## mu_diff sigma_diff nu eff_size diff_pred
## [1,] 50.01712 164.0427 29.22199 0.304903014 31.95226
## [2,] 37.03222 161.0353 32.14734 0.229963334 236.64675
## [3,] 49.31920 157.4821 22.96745 0.313173376 182.74421
## [4,] 49.97199 154.8233 27.12930 0.322767890 55.26206
## [5,] 35.29876 152.9050 46.40909 0.230854112 109.84739
## [6,] 56.79979 177.0965 35.41169 0.320727917 -109.35174
## [7,] 27.50354 165.7889 43.66613 0.165894913 -193.57416
## [8,] 70.55631 164.9705 45.39869 0.427690389 -45.81030
## [9,] 1.13776 173.2478 33.99074 0.006567243 363.65183
## [10,] 64.00716 155.3912 74.08679 0.411909724 13.46671
Wir können die Werte der Markov-Ketten nutzen, um die Kennwerte der Posterior-Verteilung nachzuvollziehen, die wir oben betrachtet haben. Wir extrahieren hierfür die 10.000 Werte für \(\mu\) aus jeder Kette und fügen sie zu einem gemeinsamen Vektor zusammen. Wir lassen uns anschließend Kennwerte für diesen Vektor ausgeben, konkret dessen Mittelwert, Median und Standardabweichung. Wie wir sehen, sind diese Werte identisch mit den Kennwerten der Posterior-Verteilung von \(\mu\), die wir oben gesehen haben.
chain1 = results_bayes$mcmc_samples[[1]][1:10000]
chain2 = results_bayes$mcmc_samples[[2]][1:10000]
chain3 = results_bayes$mcmc_samples[[3]][1:10000]
all_chains = c(chain1,chain2,chain3)
mean(all_chains)## [1] 42.51268median(all_chains)## [1] 42.45031sd(all_chains)## [1] 17.31342Unabhängige Stichproben
Wir haben bislang ein Szenario betrachtet, in dem für dieselbe Gruppe von Versuchspersonen Messwerte aus zwei Bedingungen vorliegen. Gegenstand unserer Schätzung war hierbei die mittlere Differenz der gepaarten Messwerte. Wir können mit Hilfe derselben Modellannahmen auch eine Schätzung für das Szenario unabhängiger Stichproben vornehmen. In diesem Fall lägen also Messwerte aus zwei separaten Gruppen vor. Gegenstand unserer Schätzung wäre dann die Differenz der Mittelwerte beider Gruppen. Das Verfahren entspricht dann einer bayesianischen Alternative zum t-Test für unabhängige Stichproben.
Zur Veranschaulichung stellen wir uns kurz vor, die Messwerte unserer Bedingungen \(LF\) und \(HF\) stammen aus zwei verschiedenen Gruppen. Es könnte sich beispielsweise um Reaktionszeiten aus einer lexikalischen Entscheidungsaufgabe handeln, die mit zwei Gruppen von Deutsch-Lernenden durchgeführt wurde: Fortgeschrittene Lernende (\(HF\)) und Anfänger (\(LF\)). Wir würden dann erwarten, dass die Worterkennung bei fortgeschrittenen Lernenden schneller erfolgt als bei Anfängern. Damit läge eine gerichtete Hypothese analog zu der unseres vorherigen Szenarios vor, allerdings mit unabhängigen Stichproben. Ein solches Verfahren auf Grundlage des oben erläuterten Modells sähe folgendermaßen aus.
- Die Parameter \(\mu\) und \(\sigma\) repräsentieren jeweils Mittelwert und Streuung beider Gruppen \(HF\) und \(LF\) in der Population. Entsprechend werden für beide Gruppen Posterior-Verteilungen der Parameter ermittelt, also jeweils \(\mu_{HF}\) und \(\mu_{LF}\) bzw. \(\sigma_{HF}\) und \(\sigma_{LF}\). Für beide Gruppen wird ein gemeinsamer Parameter \(\nu\) angenommen.
- Die Formen der Prior-Verteilungen für beide \(\mu\) und \(\sigma\) sind identisch mit den oben erläuterten Formen für beide Parameter (s.o. Parameter und Priors).
- Referenzpunkte für die Prior-Verteilungen der Parameter in beiden Gruppen gehen aus den jeweiligen Stichproben hervor (z.B. die Prior-Verteilung von \(\mu_{HF}\) ist um den Stichprobenmittelwert von \(HF\) zentriert und hat eine Standardabweichung entsprechend der \(s_{HF} * 1000\). Gleiches gilt für die Prior-Verteilung von \(\mu_{LF}\) und die Messwerte der Gruppe \(LF\)).
- Die Posterior-Verteilung der Gruppendifferenz, also \(\mu_{LF} - \mu_{HF}\), ergibt sich dann aus den individuellen Differenzen der aktuellen Parameterwerte beider Gruppen in der MCMC-Stichprobe.
Für eine bayesianische Schätzung der Gruppendifferenz im oben gezeigten Modell müssen wir in R lediglich den Input des Arguments paired zu FALSE ändern.
results_ind = bayes.t.test(LF, HF, paired = FALSE)
results_ind##
## Bayesian estimation supersedes the t test (BEST) - two sample
##
## data: LF (n = 100) and HF (n = 100)
##
## Estimates [95% credible interval]
## mean of LF: 712 [684, 742]
## mean of HF: 671 [653, 689]
## difference of the means: 41 [7.3, 75]
## sd of LF: 141 [120, 164]
## sd of HF: 88 [74, 102]
##
## The difference of the means is greater than 0 by a probability of 0.991
## and less than 0 by a probability of 0.009
Der Output zeigt uns die Mittelwerte der Posterior-Verteilungen von \(\mu\) und \(\sigma\) in beiden Gruppen sowie deren jeweiliges 95% HDI. Dieselbe Information erhalten wir für die Posterior-Verteilung der Differenz der Gruppenmittelwerte. Unten wird zudem erneut der Anteil der Posterior-Verteilung (Gruppendifferenz) angegeben, der über bzw. unter \(0\) liegt.
Wenn wir die Ergebnisse visualisieren, erhalten wir Histogramme für die Posterior-Verteilungen von \(\mu\) in beiden Gruppen sowie ein Histogramm der Posterior-Verteilung der Gruppendifferenz (links). Rechts unten sehen wir zudem ein Histogramm der Posterior-Verteilung der Effektstärke (Cohen’s \(d\) für unabhängige Stichproben), die sich aus den individuellen Effektstärken auf Grundlage der aktuellen Parameterwerte der MCMC-Stichprobe zusammensetzt.
plot(results_ind)
Bayesianische Korrelationsanalyse
Vorangehend haben wir den Zusammenhang zwischen Wortfrequenz und Reaktionszeit in einem faktoriellen Design betrachtet. Wir haben also eine kontinuierliche Variable (Wortfrequenz) in zwei diskrete Gruppen (\(HF\) vs. \(LF\)) eingeteilt und untersucht, wie sich die Reaktionszeit in Abhängigkeit der Gruppenzugehörigkeit eines Wortes verändert. Ein anderer gängiger Analyseansatz betrachtet wiederum zwei kontinuierliche Variablen und untersucht die Stärke ihres Zusammenhangs. Solche Zusammenhänge werden häufig über die Pearson-Korrelation (auch Produkt-Moment-Korrelation) quantifiziert. Eine Voraussetzung hierbei ist, dass der Zusammenhang von Interesse linear gedacht wird, bei Veränderung in einer Variable also eine proportional konstante Veränderung in der anderen Variable erwartet wird. Wir setzen eine solche Linearität hier der Einfachheit wegen voraus. Wir betrachten zunächst die Pearson-Korrelation anhand eines Datenbeispiels und verdeutlichen den frequentistischen Inferenzansatz. Anschließend betrachten wir erneut eine bayesianische Alternative im Bayesian-First-Aid-Paket.
Wir nehmen erneut das eingangs geschilderte Szenario zum Ausgang: Wir haben mit 100 Versuchspersonen eine lexikalische Entscheidungsaufgabe durchgeführt, in der sowohl hoch- als auch niedrigfrequente Wörter enthalten waren. Wir haben eine erwartungskonforme Tendenz beobachtet, nach der die Reaktionszeit auf niedrigfrequente Wörter länger ausfällt als auf hochfrequente Wörter. Uns interessiert nachfolgend, wie stark dieser Effekt mit dem Leseverhalten unserer Versuchspersonen zusammenhängt. Die Überlegung ist folgende: Personen, die viel und regelmäßig lesen, werden häufiger mit niedrigfrequenten Wörtern in der Schriftsprache konfrontiert und erkennen diese daher schneller als Personen, die wenig lesen. Entsprechend sollte unser Wortfrequenzeffekt bei Versuchspersonen, die viel lesen, schwächer ausgeprägt sein als bei Versuchspersonen, die wenig lesen. Wir vermuten also einen negativen Zusammenhang (s.u.).
Wir nehmen an, wir haben neben der lexikalischen Entscheidungsaufgabe eine Reihe von Testverfahren durchgeführt, die das Leseverhalten unserer Versuchspersonen quantifizieren sollen (z.B. Tests zum Wortschatzumfang, Author Recognition, etc.). Aus den Ergebnissen dieser Verfahren haben wir einen gemeinsamen, standardisierten Score zum Leseverhalten ermittelt. Der unten gezeigte Datensatz enthält für jede Versuchsperson deren individuelle Reaktionszeitdifferenz in der lexikalischen Entscheidungsaufgabe (\(LF - HF\), nachfolgend \(DIFF\)) sowie ihren standardisierten Lese-Score.
head(data)## DIFF score
## 96 358.6458 -1.9815
## 52 28.1269 -1.7345
## 97 371.1762 -0.6852
## 49 20.3566 -0.0030
## 28 -41.5336 0.0996
## 12 -153.0463 1.1368Kovarianz und Korrelation
Um die Pearson-Korrelation für unsere Variablen zu bestimmen, berechnen wir zunächst deren Kovarianz. Wir berechnen hierzu für jedes Messwertpaar das Produkt der Differenzen beider Werte vom Mittelwert der jeweiligen Variable. Wir bilden dann die Summe über alle Produkte und teilen diese durch \(N - 1\) (in unserem Fall also 99). Wir wenden diese Berechnung unten auf unser Datenbeispiel an.
#Kovarianz
cov <- sum((data$DIFF - mean(data$DIFF)) * (data$score - mean(data$score)))/nrow(data)-1
cov## [1] -59.32186
Über die Kovarianz erfahren wir, welche Richtung ein möglicher
Zusammenhang zwischen unseren Variablen hat: Beispielsweise hat die von
uns berechnete Kovarianz ein negatives Vorzeichen. Das verweist darauf,
dass wir im Zählerterm unserer Berechnung überwiegend negative Produkte
aufsummiert haben. Die individuellen Werte unserer Messwertpaare weichen
also tendenziell in entgegengesetzte Richtungen vom
Mittelwert ihrer jeweiligen Variable ab, die entsprechenden Differenzen
weisen folglich unterschiedliche Vorzeichen auf. Es liegt daher ein
negativer Zusammenhang vor: Steigen Werte in einer
Variable, fallen sie tendenziell in der anderen Variable. Eine Kovarianz
mit positivem Vorzeichen liegt hingegen dann vor, wenn die individuellen
Werte unserer Messwertpaare tendenziell gleichgerichtet von ihrem
Mittelwert abweichen. Wir haben dann überwiegend Differenzen mit
gleichem Vorzeichen und folglich überwiegend positive Produkte. Das
entspräche einem positiven Zusammenhang: Steigen Werte
in einer Variable, steigen sie tendenziell auch in der anderen
Variable.
Wir können aus der Kovarianz allein allerdings nicht die Stärke des Zusammenhangs erkennen, da der konkrete numerische Wert mit der Skalierung der individuellen Variablen zusammenhängt. Wir normalisieren die Kovarianz daher am Produkt der Standardabweichungen beider Variablen und erhalten auf diesem Weg den Korrelationskoeffizienten \(r\). Dieser ist ein standardisiertes Maß zwischen \(-1\) und \(1\) für den Umfang, in dem zwei Variablen gemeinsam variieren. Ein \(r\) von \(0\) bedeutet, dass beide Variablen praktisch unabhängig voneinander variieren, also kein Zusammenhang besteht. Ein \(r\) von \(1\) bzw. \(-1\) entspräche einem perfekten (positiven oder negativen) linearen Zusammenhang. Wir berechnen \(r\) unten manuell und kontrollieren das Ergebnis mit Hilfe des R-Befehls für die Pearson-Korrelation.
#Manuelle Berechnung des Korrelationskoeffizenten
round(
cov/sd(data$DIFF)*sd(data$score),
2)## [1] -0.35#R-Befehl
round(
cor(data$DIFF, data$score),
2)## [1] -0.35
Wir erhalten für unsere Daten ein \(r\) von \(-0.35\), was nach gängigen Interpretationskonventionen einem schwachen negativen Zusammenhang beider Variablen entspricht (moderat \(\approx \pm 0.5\), stark \(\approx \pm 0.7\)). Wir stehen erneut vor dem Problem statistischer Inferenz: Was sagt uns die Zusammenhangsstärke, die wir in unserer Stichprobe feststellen, über den Zusammenhang der Variablen in der Population?
Frequentistsche Inferenz
Im frequentistischen Ansatz schließen wir erneut mittels eines Nullhypothesen-Signifikanztest auf den Populationszustand. Da wir theoretisch motiviert nach einem negativen Zusammenhang suchen, wenden wir ein einseitiges Testszenario an. \(H_0\) sagt dann aus, dass der Zusammenhang in der Population entweder nicht existent oder positiv gerichtet ist. \(H_1\) sagt aus, dass ein negativer Zusammenhang vorliegt. Für einen solchen Test berechnen wir aus unserem Korrelationskoeffizienten einen t-Wert, der unter Annahme von \(H_0\) einer t-Verteilung mit \(N-2\) Freiheitsgraden folgt (in unserem Fall also 98). Die Entscheidungslogik ist dann dieselbe, die oben demonstriert wurde, erfolgt diesmal allerdings am unteren Ende der Verteilung anstelle des oberen.
cor.test(data$DIFF, data$score, alternative = "less")##
## Pearson's product-moment correlation
##
## data: data$DIFF and data$score
## t = -3.6688, df = 98, p-value = 0.0001982
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
## -1.0000000 -0.1931435
## sample estimates:
## cor
## -0.347511
Der Output des Tests zeigt ein signifikantes Ergebnis: Ein t-Wert kleiner oder gleich demjenigen, der aus unserem \(r\) hervorgeht, ist um ein vielfaches unwahrscheinlicher als derjenige, der die unteren 5% der Verteilung definiert (\(p\) liegt also weit unter \(\alpha = 0.05\)).
Wir vergegenwärtigen uns erneut die Inferenzlogik des oben gezeigten Verfahrens: Die t-Verteilung modelliert die Verteilung der Teststatistik in einer Nullhypothesenpopulation. Die uns vorliegende Teststatistik ist in diesem Modell extrem unwahrscheinlich, sodass wir auf Grundlage dieser Beobachtung das Modell selbst ablehnen. Erneut attribuieren wir hierbei weder \(H_0\) noch \(H_1\) selbst eine bestimmte Wahrscheinlichkeit. Stattdessen ermitteln wir, wie wahrscheinlich die uns vorliegenden Beobachtungen unter Voraussetzung der Nullhypothese sind. Das nachfolgend gezeigte, bayesianische Verfahren ermöglicht wiederum Inferenzen über die Wahrscheinlichkeit verschiedener Populationszustände selbst. Diese Inferenzen gelten erneut unter Voraussetzung spezifischer Modellannahmen, die nachfolgend erläutert werden.
Bayes-Verfahren: Parameter und Priors
Wir haben vorangehend gesehen, wie wir über die Pearson-Korrelation den Zusammenhang zweier Variablen in unserer Stichprobe quantifizieren können. Ausgehend von diesem standardisierten Zusammenhangsmaß stellen wir uns nun die Frage, welche Wahrscheinlichkeit einer gegebenen Zusammenhangsstärke beider Variablen in der Population beigemessen werden kann. Unser Korrelationskoeffizient \(r\) wird somit zu einem Parameter, für den wir eine Posterior-Verteilung suchen. Aus dieser Verteilung können wir dann Inferenzen über die Wahrscheinlichkeit verschiedener Populationszustände ableiten, wie wir es vorangehend für die Bedingungsdifferenz \(LF-HF\) im bayesianischen t-Test gemacht haben.
Der Korrelationskoeffizient wird als Populationsparameter durch \(\rho\) (rho) denotiert. Die Prior-Verteilung für \(\rho\) entspricht einer Gleichverteilung über -1 bis 1. Sämtliche möglichen Ausprägungen der Korrelation erhalten also dieselbe Prior-Wahrscheinlichkeit. Entsprechend liegt hier erneut eine uninformative Prior vor. Die weiteren im Verfahren geschätzten Parameter beschreiben die Verteilungseigenschaften der individuellen Variablen, deren Zusammenhang wir untersuchen. Die Prior-Verteilungen dieser Parameter haben Formen, die wir zuvor bereits kennengelernt haben:
\(\mu_x\) und \(\mu_y\) denotieren die Populationsmittelwerte der beiden Variablen, deren Zusammenhang \(\rho\) ausdrückt. In unserem Fall sind das die Reaktionszeitdifferenz \(LF - HF\) (bzw. \(DIFF\)) sowie der standardisierte Lese-Score. Prior-Verteilungen für beide Parameter sind breite Normalverteilungen (s.o.) um den jeweiligen Stichprobenmittelwert.
\(\sigma_x\) und \(\sigma_y\) denotieren die Streuung um den Populationsmittelwert in beiden Variablen. Prior-Verteilungen für beide Parameter sind Gleichverteilungen (s.o.) über \(s / 1000\) bis \(s * 1000\) ausgehend von der jeweiligen Stichprobenstandardabweichung.
\(\nu\) denotiert die Anzahl der Freiheitsgrade (Gewicht der Randbereiche) für die Verteilungen um die Populationsmittelwerte beider Variablen. Für beide Variablen wird dasselbe \(\nu\) angenommen. Die Prior-Verteilung für \(\nu\) ist eine Exponentialverteilung (s.o.) mit einer Wachstumsrate von \(\frac{1}{29}\).
Das Verfahren ermittelt Posterior-Verteilungen für jeden der oben genannten Parameter, wobei das MCMC-Sampling erneut jeweils mit einer Parameterkombination \(\mu_{xi}\), \(\mu_{yi}\), \(\sigma_{xi}\), \(\sigma_{yi}\), \(\nu_i\), \(\rho_i\) arbeitet. Die Prior-Dichte einer Parameterkombination ist das Produkt der individuellen Prior-Dichten aller Parameterwerte, die aus deren jeweiligen Prior-Verteilungen hervorgehen.
Likelihood
Um zu bestimmen, wie gut unsere Daten zu einer gegebenen Kombination der oben genannten Parameter passen, wird im vorliegenden Verfahren eine bivariate t-Verteilung auf Grundlage der einzelnen Parameterwerte erstellt, in die unsere Daten anschließend eingesetzt werden. Wie der Name bereits aussagt, hat eine solche Verteilung keine einzelnen Werte \(x_i\) zum Gegenstand, sondern Wertepaare \(x_i, y_i\) . Die Datenverteilung für die Variablen \(x\) und \(y\) wird wie zuvor durch die Parameter \(\mu\), \(\sigma\) und \(\nu\) modelliert. Zusätzlich hierzu modelliert allerdings der Parameter \(\rho\) das Ausmaß der Korrelation beider Variablen. Die Wahrscheinlichkeitsdichte eines Messwertpaares innerhalb der bivariaten t-Verteilung ergibt sich also nicht nur daraus, wie stark \(x_i\) und \(y_i\) jeweils von \(\mu_x\) und \(\mu_y\) abweichen, sondern zusätzlich daraus, wie sehr die gepaarten Abweichungen \(x_i - \mu_x, y_i - \mu_y\) der Systematik entsprechen, die in \(\rho\) ausgedrückt ist.
Der Einfluss des Parameters \(\rho\) auf die Likelihood wird intuitiv verständlicher, wenn wir die bivariate t-Verteilung visualisieren. Da diese Verteilung zwei Variablen zum Gegenstand hat, ist ihr Massezentrum nicht als Dichtekurve im eindimensionalen Raum zu denken (wie im univariaten Fall, s.o.), sondern als Ellipse im zweidimensionalen Raum, der durch die Variablen \(x\) und \(y\) definiert ist. In der unteren Abbildung sehen wir solche Ellipsen für den zweidimensionalen Raum unserer Variablen (Bedingungsdifferenz und Lese-Score). Die Ellipse zeigt jeweils den Wertebereich in diesem Raum, in dem 95% der Masse einer bivariaten t-Verteilung liegen. Zusätzlich sehen wir jeweils die individuellen Messwertpaare unserer Stichprobe (schwarze Punkte).
Zwischen den Abbildungen variiert allein der Parameter \(\rho\), alle anderen Parameter der Verteilung bleiben konstant. Wie wir sehen, ändert sich die Orientierung der Ellipse mit dem Vorzeichen von \(\rho\), worin sich die Richtung der Korrelation ausdrückt (links oben vs. links unten). Liegt \(\rho\) bei \(0\) (rechts oben), hat die Ellipse keine spezifische Neigung, beide Variablen variieren unabhängig voneinander. Steigt der Betrag von \(\rho\) (rechts unten), steigt die Systematik gemeinsamer linearer Variation, sodass die Ellipse schmaler wird. Der Umfang unserer Datenpunkte, der von der Ellipse erfasst wird, liefert uns ein intuitives Verständnis davon, wie gut unsere Daten zu einer bivariaten t-Verteilung mit gegebenem Wert für \(\rho\) passen.
Die Likelihood unserer Daten für eine gegebene Parameterkombination ergibt sich dann, wie zuvor, aus dem Produkt der Wahrscheinlichkeitsdichten individueller Datenpunkte (hier Messwertpaare) innerhalb der Verteilung, die durch die aktuelle Parameterkombination bestimmt ist.
Ergebnisse
Ähnlich wie beim t-Test ist auch der Code für das Korrelationsverfahren fast identisch mit der frequentistischen Variante. Wir stellen erneut lediglich bayes. voran. Wir speichern das Ergebnis des Verfahrens als Objekt, um es besser untersuchen zu können.
results_cor <- bayes.cor.test(data$score, data$DIFF)
results_cor##
## Bayesian First Aid Pearson's Correlation Coefficient Test
##
## data: data$score and data$DIFF (n = 100)
## Estimated correlation:
## -0.34
## 95% credible interval:
## -0.52 -0.16
## The correlation is more than 0 by a probability of <0.001
## and less than 0 by a probability of >0.999
Der Output zeig uns die beste Schätzung der Korrelation in der Population auf Grundlage des Medians der Posterior-Verteilung von \(\rho\). Diese liegt bei \(-0.34\). Wir sehen zudem das 95%-HDI der Posterior-Verteilung, das vollständig unterhalb von \(0\) verläuft. Für unsere Inferenz auf den Populationszustand bedeutet das, dass wir mit hoher Posterior-Wahrscheinlichkeit einen negativen Zusammenhang in der Population annehmen können. Dem entspricht auch die Angabe zur Wahrscheinlichkeit einer Korrelation größer bzw. kleiner \(0\). Die Wahrscheinlichkeit, dass die Korrelation kleiner als \(0\), also negativ ist, strebt praktisch gegen 1.
Mit dem Plot-Befehl können wir die Ergebnisse der Parameterschätzung visualisieren. Die obere Abbildung zeigt uns ein Histogramm der Posterior-Verteilung von \(\rho\). Wie im Output bereits ersichtlich war, verläuft die Verteilung praktisch vollständig im negativen Bereich. Die untere Abbildung zeigt die individuellen Messwertpaare unserer Stichprobe gegen eine bivariate t-Verteilung, die jeweils auf den besten Schätzungen für alle untersuchten Parameter basiert (also \(\mu_{DIFF}\), \(\mu_{score}\), \(\sigma_{DIFF}\), \(\sigma_{score}\), \(\nu\) und \(\rho\)). Die dunkel gefärbte Fläche schließt den Bereich ein, in dem 50% der Masse dieser Verteilung liegen. Die helle Fläche schließt 95% der Verteilungsmasse ein.
plot(results_cor)
Wie zuvor im t-Test können wir uns die Kennwerte der MCMC-Stichprobe für alle Parameter ausgeben lassen.
results_cor$stats[1:6,]## mean sd HDI% comp HDIlo HDIup
## rho -0.341630709 0.09064594 95 0 -0.5162114 -0.1609878
## mu[1] 0.008276271 0.10266365 95 0 -0.1914469 0.2090729
## mu[2] 42.737675778 17.15443668 95 0 7.9662295 75.0918536
## sigma[1] 0.996127398 0.07738925 95 0 0.8436346 1.1467373
## sigma[2] 168.357500406 12.77636314 95 0 144.8806855 194.6247309
## nu 49.022830334 31.83494303 95 0 5.8932583 112.4754462
## %>comp %<comp q2.5% q25% median
## rho 0.0002666311 9.997334e-01 -0.5093942 -0.40371115 -0.344854127
## mu[1] 0.5337288362 4.662712e-01 -0.1936923 -0.06030162 0.008110564
## mu[2] 0.9918677510 8.132249e-03 9.0015942 31.32473164 42.700417778
## sigma[1] 0.9999333422 6.665778e-05 0.8534219 0.94244236 0.991675038
## sigma[2] 0.9999333422 6.665778e-05 144.8836073 159.35174602 167.913190997
## nu 0.9999333422 6.665778e-05 11.2047000 25.99079636 41.278304889
## q75% q97.5% mcmc_se Rhat n_eff
## rho -0.28302116 -0.1504762 0.0010381556 1.0012926 7626
## mu[1] 0.07707568 0.2075141 0.0012446357 1.0007288 6822
## mu[2] 54.19890931 76.2265979 0.2033184840 1.0005963 7119
## sigma[1] 1.04591892 1.1588742 0.0009335485 1.0008328 6965
## sigma[2] 176.59286138 194.6260708 0.1492801601 0.9999192 7334
## nu 63.39081688 130.7136136 0.4906359696 1.0002385 4262
Ebenso können wir erneut auf die einzelnen MCMC-Ketten zugreifen. Unten beispielhaft auf die ersten 10 Werte der ersten Kette.
results_cor$mcmc_samples[1][1:10,1:6]## [[1]]
## rho mu[1] mu[2] sigma[1] sigma[2] nu
## [1,] -0.2950235 -0.037701426 59.23503 1.0650165 157.7701 99.15356
## [2,] -0.3669400 -0.046588729 57.59017 1.0828447 168.8477 128.43758
## [3,] -0.3249030 0.077921443 45.43521 1.1109512 173.6796 36.16859
## [4,] -0.3239248 0.045250242 55.71033 1.0624964 165.1642 48.24208
## [5,] -0.4846882 0.098121236 55.40229 0.9408826 165.4641 60.30894
## [6,] -0.3490578 0.009413685 57.91542 1.0242562 196.8899 67.68546
## [7,] -0.3431469 -0.097688851 58.28678 0.9789553 180.5275 78.92540
## [8,] -0.4533610 -0.026760872 52.09857 0.9872667 168.0727 23.87391
## [9,] -0.3368892 -0.001529006 33.29660 0.9488904 170.2038 50.45260
## [10,] -0.3435254 0.024620926 46.05560 0.9553054 157.5370 31.39889Bayesianische Proportionsanalyse
Wir werden abschließend betrachten, wie wir mit Hilfe bayesianischer Parameterschätzung auf Unterschiede in der Auftretenswahrscheinlichkeit eines bestimmten Ereignisses zwischen zwei Gruppen schließen können. Wir gehen hierbei vom Ergebnis eines Bernoulli-Versuchs aus (Erfolg vs. Misserfolg) und beobachten, dass sich die Erfolgsrate zwischen zwei Gruppen unserer Stichprobe unterscheidet. Ausgehend von dieser Beobachtung wollen wir nun auf den Populationszustand bezüglich der Erfolgsrate in beiden Gruppen schließen. Wir bleiben in unserem Beispiel bei der Frage nach dem Zusammenhang von Leseverhalten und Worterkennung. Vorangehend haben wir festgestellt, dass der Effekt der Wortfrequenz auf die Geschwindigkeit, mit der ein Wort erkannt wird, negativ mit dem Leseverhalten unserer Versuchspersonen korreliert. Uns interessiert nun zusätzlich, ob das Leseverhalten auch einen Einfluss darauf hat, ob niedrigfrequente Wörter überhaupt als Wörter erkannt werden.
Zur Erinnerung: Wir haben eine lexikalische Entscheidungsaufgabe durchgeführt. Versuchspersonen müssen hierbei per Tastendruck entscheiden, ob eine gegebene Buchstabensequenz ein deutsches Wort ergibt oder nicht. Von primärem Interesse ist hierbei die Geschwindigkeit, mit der ein echtes Wort erkannt wird. Fehlentscheidungen, in denen ein echtes Wort fälschlicherweise als Nicht-Wort eingeordnet wird, sind zwar prinzipiell möglich, aber unwahrscheinlich und daher weniger informativ. Aber wie verhält es sich mit einem wirklich seltenen Wort, beispielsweise Dachtel (= Ohrfeige)? Wir nehmen dieses Extrembeispiel als Testdomäne für unsere Frage: Angenommen, das Wort befand sich in der Gruppe niedrigfrequenter Wörter unserer Aufgabe. Wurde das Wort von Personen, die viel lesen, häufiger als deutsches Wort erkannt?
Um Erfolgsraten für die Erkennung des Wortes in zwei Gruppen bestimmen zu können, teilen wir unsere Stichprobe kategorial in häufige und seltene Leser ein, ausgehend vom Median des Lese-Scores. Der unten gezeigte Datenausschnitt zeigt für mehrere Versuchspersonen deren Gruppenzuordnung sowie die Kategorisierung ihrer Antwort auf das Wort Dachtel (richtig oder falsch = r vs. f ).
head(data[,3:5])## Leser Wort Antwort
## 96 wenig Dachtel f
## 52 wenig Dachtel f
## 97 wenig Dachtel f
## 49 wenig Dachtel f
## 28 viel Dachtel r
## 12 viel Dachtel f
Wir erstellen aus unseren Daten eine Kontingenztabelle, die für beide Gruppen die Anzahl beider Antworttypen enthält. Wie wir sehen können, wurde das Wort in beiden Gruppen mehrheitlich nicht erkannt. Allerdings liegen unter den häufigen Lesern mehr korrekte Antworten vor, als unter den seltenen Lesern. Können wir hiervon ausgehend auf einen entsprechenden Unterschied in der Population schließen?
contingencies <- table(data$Leser, data$Antwort)
contingencies##
## r f
## viel 16 34
## wenig 7 43Frequentistische Inferenz
Der frequentistische Proportionstest verfährt zunächst wie ein Chi-Quadrat-Test auf stochastische Unabhängigkeit. Die Chi-Quadrat-Statistik \(\chi^2\) quantifiziert, wie stark die gruppenweisen Häufigkeiten eines Ereignisses (Zellenwerte der Kontingenztabelle) von den erwarteten Häufigkeiten unter Annahme stochastischer Unabhängigkeit abweichen. Die erwartete Häufigkeit entspricht dem Produkt der Randhäufigkeiten (Zeilensumme*Spaltensumme) normalisiert an der Anzahl an Beobachtungen. Die Teststatistik drückt nur die Größe der Abweichung aus, nicht deren Richtung. \(\chi^2\) bildet daher immer die Grundlage für einen ungerichteten Test. Wir könnten in unserem Fall also testen, ob sich die Erfolgsrate zwischen den Gruppen unterscheidet, aber nicht explizit, ob das kritische Wort von häufigen Lesern verlässlicher erkannt wurde.
Werden allerdings nur zwei Gruppen verglichen, wie im vorliegenden Fall, ist ein gerichteter Test möglich: Die Quadratwurzel aus \(\chi^2\) folgt einer Standartnormalverteilung (\(z\)-Verteilung). Wir können aus \(\chi^2\) also einen z-Wert ableiten, mit Vorzeichen entsprechend dem der Differenz in der Erfolgsrate beider Gruppen. Die Teststatistik ist dann also symmetrisch um \(0\) verteilt, sodass sowohl ein unterer als auch ein oberer Ablehnungsbereich für \(H_0\) definiert werden kann. Wir wenden den Proportionstest unten auf unsere Kontingenztabelle an.
prop.test(contingencies, alternative = "greater")##
## 2-sample test for equality of proportions with continuity correction
##
## data: contingencies
## X-squared = 3.6138, df = 1, p-value = 0.02865
## alternative hypothesis: greater
## 95 percent confidence interval:
## 0.02476142 1.00000000
## sample estimates:
## prop 1 prop 2
## 0.32 0.14
Der Output des Verfahrens zeigt uns die Chi-Quadrat-Statistik unserer Daten. Der daneben berichtete p-Wert bezieht sich allerdings auf den z-Wert, der aus \(\chi^2\) hervorgeht (\(z = -1.9\)). Unten sehen wir zudem die Erfolgsraten beider Gruppen in unserer Stichprobe. Wie wir sehen können, liegt ein signifikantes Ergebnis vor. Unter Annahme der Nullhypothese, dass die Erfolgsraten in beiden Gruppen nicht oder anders gerichtet divergieren (seltene Leser > häufige Leser), ist die vorliegende Teststatistik hinreichend unwahrscheinlich (\(\alpha = 0.05\)), um das Nullhypothesen-Populationsmodell zu verwerfen. Wie in den vorherigen Beispielen attribuieren wir erneut keinem Populationszustand (\(H_0\) vs. \(H_1\)) eine bestimmte Wahrscheinlichkeit, sondern unseren Daten unter Voraussetzung eines bestimmten Populationszustands (\(H_0\))
Bayes-Verfahren: Parameter und Priors
Ziel des bayesianischen Verfahrens ist eine Inferenz auf die Wahrscheinlichkeit, dass sich die Erfolgsraten in beiden Gruppen in der von uns angenommenen Richtung unterscheiden. Wir suchen also eine Posterior-Verteilung für die Differenz der Erfolgsraten beider Gruppen. Wir erhalten diese Verteilung, indem wir separate Posterior-Verteilungen für die Erfolgsraten in beiden Gruppen ermitteln. Die Posterior-Verteilung der Differenz zwischen den Erfolgsraten beider Gruppen wir dann ermittelt, indem für jedes MCMC-Sample die Differenz zwischen den geschätzten Erfolgsraten berechnet wird.
Anders als in den vorherigen Verfahren schätzen wir pro Gruppe \(x\) und \(y\) nur einen einzigen Populationsparameter, nämlich die Erfolgsrate \(\theta_x\) bzw. \(\theta_y\). Die Prior-Verteilung für \(\theta\) ist eine beta-Verteilung, eine Wahrscheinlichkeitsverteilung über \([0,1]\) mit den Parametern \(\alpha\) und \(\beta\). Wir können uns die beta-Verteilung als stetige Verteilung über einen Wahrscheinlichkeitsraum vorstellen, in der das Verhältnis der Parameter \(\alpha\) und \(\beta\) das Verhältnis von Erfolgen zu Misserfolgen enkodiert. Hieraus ergibt sich die konkrete Form der Verteilung, wie unten anhand einiger Beispiele veranschaulicht:
Bei \(\alpha = \beta\) liegt eine symmetrische Verteilung mit Erwartungswert bei \(0.5\) vor (links unten), bei Ungleichheit verlagert sich die Verteilungsmasse schief in Richtung \(0\) (\(\alpha < \beta\) = mehr Misserfolge als Erfolge, links oben) bzw. \(1\) (\(\alpha > \beta\) = mehr Erfolge als Misserfolge, rechts oben). Das Gewicht der Randbereiche ergibt sich aus der größe von \(\alpha\) und \(\beta\): Je kleiner beide Parameter, desto flachgipfliger die Verteilung (vgl. links unten mit den oberen Abbildungen). Hierbei bildet der Fall \(\alpha = \beta = 1\) (rechts unten) den Extremfall, in dem eine Gleichverteilung über dem gesamten Wahrscheinlichkeitsraum vorliegt. Diese spezielle Form der beta-Verteilung wird im vorliegenden Verfahren als Prior-Verteilung für die Erfolgsrate \(\theta\) verwendet. Sämtliche denkbaren Erfolgsraten erhalten a-priori also dieselbe Wahrscheinlichkeit, weshalb erneut eine uninformative Prior vorliegt.
Likelihood
Die Bestimmung der Likelihood ist im vorliegenden Fall weniger komplex als in den vorangehenden Verfahren. Unsere Datengrundlage bilden Erfolgshäufigkeiten bei einer festen Anzahl an Beobachtungen (hier jeweils 50 pro Gruppe). Wir wollen ermitteln, wie gut diese Daten zu einer gegebenen Erfolgsrate \(\theta_i\) passen. Die Datenverteilung gegeben \(\theta_i\) ist über eine Binomialverteilung modellierbar, also eine diskrete Wahrscheinlichkeitsverteilung über die Anzahl an Erfolgen \(k\) bei einer festen Anzahl an Beobachtungen \(n\) und einer Erfolgswahrscheinlichkeit \(p\). Letztere wird durch \(\theta_i\) repräsentiert, \(k\) und \(n\) entnehmen wir unseren Daten.
Wir betrachten zur Veranschaulichung die Daten aus der Gruppe der häufigen Leser: Das kritische Wort wurde von 16 der 50 Versuchspersonen erkannt. Es liegen also \(k = 16\) Erfolge bei \(n = 50\) Beobachtungen vor. Die Likelihood dieser Daten gegeben eine zufällige Erfolgsrate \(\theta_i = 0.25\) entspräche
\[P(k = 16| p = 0.25) = \binom{50}{16} (0.25)^{16} (0.75)^{34} = 0.0647\]
also ungefähr einer Wahrscheinlichkeit von \(0.065\). Anders als zuvor arbeiten wir hier also mit echten Wahrscheinlichkeiten, nicht mit Wahrscheinlichkeitsdichten, da wir unsere Daten mit einer diskreten Wahrscheinlichkeitsverteilung modellieren.
Die unstandardisierte Posterior-Dichte (Prior*Likelihood) einer Parameterkombination \(\theta_{xi}\), \(\theta_{yi}\), die im MCMC-Verfahren berechnet wird, ergibt sich aus dem Produkt der Prior-Dichten von \(\theta_{xi}\) und \(\theta_{yi}\) multipliziert mit dem Produkt der Likelihoods \(L(D_x|\theta_{xi})\) und \(L(D_y|\theta_{yi})\).
Ergebnisse
Der Code für die bayesianische Proportionsanalyse entspricht dem des frequentistischen Tests, es wird lediglich bayes. vorangestellt. Wir speichern das Ergebnis als Objekt, um seine Eigenschaften besser untersuchen zu können.
results_prop <- bayes.prop.test(contingencies)
results_prop##
## Bayesian First Aid proportion test
##
## data: contingencies
## number of successes: 16, 7
## number of trials: 50, 50
## Estimated relative frequency of success [95% credible interval]:
## Group 1: 0.33 [0.20, 0.46]
## Group 2: 0.15 [0.065, 0.25]
## Estimated group difference (Group 1 - Group 2):
## 0.18 [0.012, 0.33]
## The relative frequency of success is larger for Group 1 by a probability
## of 0.983 and larger for Group 2 by a probability of 0.017 .
Der Output zeigt uns für beide Gruppen den Median sowie das 95%-HDI der Posterior-Verteilung der Erfolgsrate. Auf Grundlage dieser Verteilungen erhalten wir die Posterior-Verteilung der Gruppendifferenz. Deren beste Schätzung (Median) liegt bei \(0.18\), mit einem 95%-HDI oberhalb \(0\). Wir können also inferieren, dass die Erfolgsrate in Gruppe 1, in unserem Fall sind das die häufigen Leser, in der Population mit hoher Wahrscheinlichkeit größer ist als die Erfolgsrate in Gruppe 2, also in der Gruppe der niedrigen Leser. Dem entspricht die untere Angabe zum Anteil der Posterior-Verteilung über bzw. unter \(0\).
Wenn wir die Ergebnisse über den Plot-Befehl visualisieren, erhalten wir Histogramme der Posterior-Verteilungen von \(\theta\) in beiden Gruppen (grün) sowie Histogramme der Posterior-Verteilung der Gruppendifferenz, jeweils in beide Richtungen (blau).
plot(results_prop)
Wir können erneut das Stats-Objekt unserer Ergebnisse aufrufen. Darin sehen wir die Kennwerte der Posterior-Stichproben für die Erfolgsraten in beiden Gruppen sowie deren Differenz.
results_prop$stats[c(1,2,5),]## mean sd HDI% comp HDIlo HDIup
## theta[1] 0.3283562 0.06526781 95 0.5 0.20194818 0.4577822
## theta[2] 0.1534410 0.04949590 95 0.5 0.06456821 0.2517431
## theta_diff[1,2] 0.1749151 0.08160787 95 0.0 0.01225492 0.3318220
## %>comp %<comp q2.5% q25% median
## theta[1] 5.799227e-03 0.99420077 0.20615729 0.2827092 0.3255341
## theta[2] 6.665778e-05 0.99993334 0.07063092 0.1174064 0.1488301
## theta_diff[1,2] 9.831356e-01 0.01686442 0.01379394 0.1200673 0.1750130
## q75% q97.5% mcmc_se Rhat n_eff
## theta[1] 0.3715331 0.4634410 0.0005222148 1.001803 15800
## theta[2] 0.1841665 0.2618128 0.0004358997 1.000861 12974
## theta_diff[1,2] 0.2303647 0.3348201 NA NA NA
Ebenso können wir auf die einzelnen MCMC-Ketten zugreifen, hier erneut exemplarisch auf die ersten 10 Zeile der ersten Kette. Darin sehen wir einzelne MCMC-Samples für die Erfolgsraten in beiden Gruppen.
results_prop$mcmc_samples[1][1:10,1:2]## [[1]]
## theta[1] theta[2]
## [1,] 0.2946595 0.16093760
## [2,] 0.2406945 0.13510567
## [3,] 0.3124596 0.21719481
## [4,] 0.3178707 0.22875436
## [5,] 0.2628023 0.15192963
## [6,] 0.2654938 0.21096604
## [7,] 0.3412009 0.23758146
## [8,] 0.4548585 0.08168685
## [9,] 0.2876865 0.16085222
## [10,] 0.2620113 0.14809148Weiterführende Ressourcen
Bayesianische Multi-level Regressionsmodelle
Aufbauend auf den hier eingeführten Konzepten bietet sich als nächster Schritt die Beschäftigung mit bayesianischen Mehrebenen-Regressionsmodellen an. Die exzellente brms R-library bietet dafür einen sehr flexiblen Zugang. Am besten gelingt dieser Schritt wohl über die Website des Autors von brms, Paul Buerkner, auf der Anleitungen zur Installation und erste Schritte zu finden sind.
Ebenfalls empfehlenswert sind die frei verfügbaren Online-Kurse, die auf Richard McElreaths einflussreichem Buch “Statistical Rethinking” aufbauen, darunter etwa die auf YouTube verfügbaren Vorlesungen des Autors selbst. Beispielsweise hier.