Fallstudie ‘Filmbewertungen’
ses
2021-12-09 17:31:41
- 1 Vorbereitung
- 2 Hintergrund
- 3 Explorative Analyse
- 4 Modelldefinition
- 5 Modell in R
- 5.1 Modell 1: Standard-Priors (
post1) - 5.2 Modell 2: Informierte (?) Priors
- 5.2.1 Modellefinition in R (rstanarm)
- 5.2.2 Priors
- 5.2.3 Prior-Prädiktiv-Verteilung
- 5.2.4 Posteriori-Verteilung
- 5.2.5 Visualisieren von
post2 - 5.2.6 Quantile
- 5.2.7 Wahrscheinlichkeiten für Parameterwerte
- 5.2.8 \(R^2\)
- 5.2.9 Verhältnis von Prior zu Posteriori
- 5.2.10 PPV
- 5.2.11 90%-Vorhersage-Intervalle
- 5.3 Modell 3
- 5.1 Modell 1: Standard-Priors (
- 6 Modellvergleich
- 7 Fazit
- 8 SessionInfo
1 Vorbereitung
1.1 R-Pakete
library(rstanarm) # Bayes-Modelle
library(tidyverse) # Datenjudo
library(bayesplot) # Plotting
library(gt) # Tabellen
library(parallel) # Turbo
library(rstatix) # Deskriptive Statistiken
library(bayestestR) # Vernachlässigbare Unterschiede/Zusammenhänge
#library(see) # Visualisierung
library(tictoc) # Zeit messen, wie lange das Modell rechnetTurbo einschalten:
options(mc.cores = parallel::detectCores())1.2 Daten: Filmbeurteilung
library(ggplot2movies)
data(movies)Hilfe zu den Daten gibt es hier:
help(movies)2 Hintergrund
2.1 Bezug zum Studiengang AWM
Nach dem Studium haben Sie bei einem Online-Händler angeheuert und zwar in der Medien-Abteilung. Der Hinweis, dass Sie etwas mit Medien studiert haben genügte. Vielleicht war auch die Behauptung, dass Sie der absolute R-Checker seien nützlich, um den Job zu bekommen … Jedenfalls müssen Sie den Behauptungen Taten folgen lassen, Schluck!
2.2 Forschungsfrage
Wie (stark) ist der Zusammenhang von logarithmierten Budget und Bewertung eines Films?
Unser Kunde – ein reicher Mäzen mit extremen Lebenswandel, der gerne ein paar Millionen investieren möchte – geht von einem positivem Zusammenhang, \(\beta\) aus. Entsprechend sei unsere Hypothese \(H_1: \beta > 0\).
3 Explorative Analyse
3.1 Datensatz vorverarbeiten
Es macht für Analysen, die in Stan laufen, grundsätzlich Sinn:
- nur die relevanten Variablen in die Analyse aufzunehmen
- fehlende Werte zu entfernen
- vorab (vor dem Modellieren) etwaige Transformationen (wie Quadrieren, logarithmieren) vorzunehmen
- UV und AV zu z-standardisieren
movies <-
movies %>%
mutate(budget_log10 = log10(budget))
movies2 <-
movies %>%
select(budget_log10, rating) %>%
drop_na() %>%
filter(budget_log10 != -Inf)Einige Filme haben ein Budget von 0. Das Logarithmus von 0 ist aber minus Unendlich. Mit dieser “Zahl” kann man nicht rechnen. Daher filtern wir alle Zeilen mit -Inf heraus.
Eine Standardisierung der Prädiktoren könnte nützlich sein, ist hier aber nicht weiter ausgeführt.
3.2 Logarithmus
budget_log10 fasst die Größenordnung des Budgets:
- 1000: 10^3 -> log10(10^3) = 3
- 10000: 10^4 -> log10(10^4) = 4
- 100000: 10^5 -> log10(10^5) = 5
Wir modellieren also die Größenordnung, nicht den Betrag selber.
3.3 Deskriptive Statistiken
movies2 %>%
get_summary_stats() %>%
gt()| variable | n | min | max | median | q1 | q3 | iqr | mad | mean | sd | se | ci |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| budget_log10 | 5183 | 3 | 8.301 | 6.477 | 5.439 | 7.176 | 1.737 | 1.222 | 6.232 | 1.204 | 0.017 | 0.033 |
| rating | 5183 | 1 | 10.000 | 6.300 | 5.200 | 7.200 | 2.000 | 1.483 | 6.137 | 1.544 | 0.021 | 0.042 |
Die Umrechung (im Kopf) von logarithmierten Werten ist nicht ganz einfach, daher hier nochmal die ursprünglichen, nicht-transformierten Werte:
movies %>%
select(budget, rating) %>%
get_summary_stats(type = "mean_sd")## # A tibble: 2 × 4
## variable n mean sd
## <chr> <dbl> <dbl> <dbl>
## 1 budget 5215 13412513. 23350085.
## 2 rating 58788 5.93 1.553.4 Zusammenhang im Datensatz visualisieren
plot0 <- movies %>%
ggplot(aes(x = budget_log10, y = rating)) +
geom_point(alpha = .2)
plot0 +
geom_smooth(method = "lm")
Es scheint einen sehr schwachen, negativen Zusammenhang zu geben.
Das gefällt unserem Auftraggeber nicht.
Aber, sagt er, und das zu Recht, entscheidend ist ja nicht die Stichprobe, sondern die Population. Weswegen wir uns bitte schön sofort an die Inferenzstatistik bewegen sollten.
Aha, unser Mäzen kennt sich also sogar mit Statistik aus.
4 Modelldefinition
4.1 Likelihood
- \(r_i \sim N(\mu_i, \sigma)\)
mit \(r_i\): Rating für Film \(i\)
4.2 Lineares Modell
- \(\mu_i = \alpha + \beta_1 b\)
4.3 Prioris
Die Prioriwerte sind etwas weiter unten aufgeführt.
5 Modell in R
5.1 Modell 1: Standard-Priors (post1)
5.1.1 Modellefinition in R (rstanarm)
tic()
post1 <- stan_glm(rating ~ budget_log10,
data = movies2,
refresh = 0 # Mit `refresh = 0` bekommt man nicht so viel Ausgabe
)
toc()## 1.766 sec elapsedEinige Infos zu den Priori-Werten bei stan_glm() findet sich hier.
Einen Überblick zu allen Funktionen (R-Befehlen) aus rstanarm findet sich übrigens hier.
tic() und toc() messen Anfang und Ende der dazwischen liegenden Zeit. Da die Modelle ein paar Sekunden brauchen, ist es ganz interessant zu wissen, wie lange wir warten mussten.
5.1.2 Priors
Welche Priori-Werte wurden (per Standard) gewählt?
prior_summary(post1)## Priors for model 'post1'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 6.1, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 6.1, scale = 3.9)
##
## Coefficients
## Specified prior:
## ~ normal(location = 0, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 0, scale = 3.2)
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.65)
## ------
## See help('prior_summary.stanreg') for more detailsMit coefficients ist das Regressionsgewicht \(\beta\) gemeint.
Das sind keineswegs besonders schlaue Prioris. Könnten wir unseren Mäzen befragen, er scheint ja Experte zu sein, kämen wir vielleicht zu klügeren Prioris (allerdings ist der Mäzen gerade auf einer “Musen-Reise” und nicht zu sprechen).
5.1.3 Posteriori-Verteilung
Überblick über die Parameter:
print(post1)## stan_glm
## family: gaussian [identity]
## formula: rating ~ budget_log10
## observations: 5183
## predictors: 2
## ------
## Median MAD_SD
## (Intercept) 6.8 0.1
## budget_log10 -0.1 0.0
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 1.5 0.0
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregLangfassung:
summary(post1)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: rating ~ budget_log10
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 5183
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 6.8 0.1 6.6 6.8 6.9
## budget_log10 -0.1 0.0 -0.1 -0.1 -0.1
## sigma 1.5 0.0 1.5 1.5 1.6
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 6.1 0.0 6.1 6.1 6.2
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3907
## budget_log10 0.0 1.0 3940
## sigma 0.0 1.0 3201
## mean_PPD 0.0 1.0 3668
## log-posterior 0.0 1.0 1688
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Nur die mittleren Schätzwerte für die Regression:
coef(post1)## (Intercept) budget_log10
## 6.7643587 -0.1009524Man kann sich die Posteriori-Intervalle so ausgeben lassen:
posterior_interval(post1) %>%
round(2)## 5% 95%
## (Intercept) 6.58 6.95
## budget_log10 -0.13 -0.07
## sigma 1.51 1.56Wir wissen jetzt also schon das Wesentliche: Mit einer Wahrscheinlichkeit von 90% liegt der Zusammenhang (das Beta-Gewicht) für Log-Budget zwischen -0.13 und -0.07. Es gibt also einen gewissen, negativen Zusammenhang zwischen Budget und Bewertung eines Films, laut unserem Golem zumindest.
5.1.4 Visualisieren von post1
5.1.4.1 Priori-Verteilung
post1_prior_pred <- stan_glm(rating ~ budget_log10,
data = movies2,
prior_PD = TRUE # DIESER Schalter gibt uns die Prior-Prädiktiv-Verteilung
#, refresh = 0
) # Mit `refresh = 0` bekommt man nicht so viel AusgabeAus der Hilfeseite:
prior_PD A logical scalar (defaulting to FALSE) indicating whether to draw from the prior predictive distribution instead of conditioning on the outcome.
post1_prior_pred## stan_glm
## family: gaussian [identity]
## formula: rating ~ budget_log10
## observations: 5183
## predictors: 2
## ------
## Median MAD_SD
## (Intercept) 5.6 20.2
## budget_log10 0.2 3.2
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 1.1 1.1
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregDie Koeffizienten aus post1_prior_pred sind also rein durch die Priori-Werte definiert; die Daten sind nicht eingeflossen.
Wenn wir das Objekt mit as_tibble() in eine Tabelle umwandeln, bekommen wir eine Tabelle mit den Stichproben:
post1_prior_pred_draws <-
post1_prior_pred %>%
as_tibble() %>%
rename(a = `(Intercept)`, # schönere, einfachere Namen
b = budget_log10) %>%
slice_sample(n = 100)movies2 %>%
ggplot() +
geom_point(aes(x = budget_log10, y = rating)) +
geom_abline(data = post1_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2)
Puh, die Priori-Werte sind … vogelwild 🐦.
5.1.4.2 Posteriori-Verteilung: Regressionsgerade
plot1 <- plot0 +
geom_abline(intercept = coef(post1)[1],
slope = coef(post1)[2],
color = "blue")
plot1
col_names <- c("a", "b", "sigma")
draws_m1 <-
post1 %>%
as_tibble()
names(draws_m1) <- col_namesEin Blick in die ersten paar Zeilen der Post-Stichproben:
draws_m1 %>%
slice_head(n=10) %>%
gt() %>%
fmt_number(everything(), decimals = 1)| a | b | sigma |
|---|---|---|
| 6.9 | −0.1 | 1.5 |
| 6.8 | −0.1 | 1.5 |
| 6.7 | −0.1 | 1.6 |
| 6.7 | −0.1 | 1.5 |
| 7.0 | −0.1 | 1.6 |
| 6.6 | −0.1 | 1.5 |
| 6.9 | −0.1 | 1.6 |
| 6.7 | −0.1 | 1.5 |
| 6.9 | −0.1 | 1.5 |
| 6.7 | −0.1 | 1.5 |
Und hier die Posteriori-Verteilung für \(\alpha\) und \(\beta\) visualisiert:
plot0 +
geom_abline(data = draws_m1,
aes(intercept = a,
slope = b),
color = "skyblue1",
alpha = .1) +
geom_abline(intercept = coef(post1)[1],
slope = coef(post1)[2],
color = "blue")
5.1.4.3 Verteilung von \(\beta\)
Hier die Verteilung für die Steigung (Regressionsgewicht \(\beta\)):
draws_m1 %>%
ggplot(aes(x = b)) +
geom_density()
5.1.4.4 Posterior-Intervalle
Die Posteriori-Intervalle kann man sich schnöde mit plot() ausgeben lassen, wobei man als Parameter den Namen des Modells übergibt.
plot(post1)
Es gibt aber auch andere Darstellungsarten, z.B. als Dichtediagramme:
mcmc_areas(post1) +
labs(title = "Posteriori-Verteilung",
caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")
S. Hilfe hier
Man kann sich auch angeben, welchen Parameter man visualisieren möchte, und zwar mit dem Argument pars (wie parameters).
Oder nur das Posteriori-Interval für den Regressionskoeffizienten:
mcmc_areas(post1,
pars = "budget_log10")
5.1.5 Fazit
Die Wahrscheinlichkeit, dass der Zusammenhang in Wirklichkeit zwischen -0.15 und -0.05 ist, ist sehr hoch.
Aber was bedeutet das, wie interpretiert man den Befund?
Rufen wir uns dazu ins Gedächtnis, was log10 bedeutet: Es bedeutet, dass man das Budget mit 10 multipliziert, also verzehnfacht, sozusagen eine Null hintendran schreibt.
Also: Verzehnfacht man das Budget, so verringert sich mittlere Bewertung um etwa 0.15 bis 0.05 Punkte.
Ob das viel ist, muss ein Medienexperte beantworten. Vielleicht sind Sie ja einer!
5.2 Modell 2: Informierte (?) Priors
Kramen wir all unser Wissen über Filme und die Filmindustrie zusammen! Wie wichtig sind die Kohlen für die Güte eines Films? Haben die einen großen Einfluss?
(Wenn wir von “Einfluss” sprechen, denken wir sicher automatisch an kausalen Einfluss. Jedenfalls geht es mir so und ich glaube, es ist schwierig, nicht an kausalen Einfluss, sondern nur an blanke Assoization, zu denken).
5.2.1 Modellefinition in R (rstanarm)
Außerdem z-standardisieren wir jetzt UV und AV.
movies3 <-
movies2 %>%
transmute(budget_log10_z = scale(budget_log10),
rating_z = scale(rating))Nach langem Beratschschlagen mit Film-Experten gehen wir davon aus, dass es im Prinzip - “normalerweise” keinen Zusammenhang gibt. Allerdings hat man in Filmen schon Pferde kotzen gesehen, deswegen wollen wir unseren Golem zu verstehen geben, dass er sich auch auf Überraschungen einstellen soll … (Schließlich wurden wir, oder zumindest unser Mäzen gerade selber von den Daten überrascht.)
Kurz gesagt, wir definieren \(\beta\) so, auf Golem-Sprech:
\[\beta = \mathcal{N}(0,0.2)\]
Bei z-standardisierten Variablen in einem einfachen Regressionsmodell kann \(\beta\) wir \(\rho\) (Korrelation) interpretiert werden.
Oben formulierter Prior ist also der Meinung.
Es gibt eine schwache (\(r=0.2\)) Korrelation zwischen Verzehnfachung des Budgetsund der Filmbewertung.
Mit anderen Worten, ein Unterschied von einer SD im Log-Budget geht mit einer Veränduerng von (in 95% der Fälle) nicht mehr als ±0.4 Rating-Punkte einher.
Ein Nachteil der z-Standardisierung ist die schwierigere Interpretation; daher ist die Zentrierung in einigen Fällen die bessere Transformation. Bei der Zentrierung bleiben die Einheiten (wie Dollar) erhalten, was oft wünschenswert ist. Sind die Eiheiten aber wenig “sprechend” (z.B. Punkte auf einer Fragebogenskala), so verliert man wenig, wenn man die Roh-Einheiten mit der z-Transformation aufgibt.
Für die übrigen Priorwerte greifen wir auf wenig ambitionierte Standardwerte zurück.
transmute() ist wie mutate(), nur dass nur die transformierten Variablen im Datensatz behalten werden, alle übrigen Variablen werden entfernt.
tic()
post2 <- stan_glm(rating_z ~ budget_log10_z,
data = movies3,
prior = normal(0, .2), # beta
refresh = 0) # Nicht so viel Ausgabe
toc()## 1.747 sec elapsedEin “Gefühl” zur Exponentialverteilung zu entwickeln, ist nicht so einfach wie bei der Normalverteilung, da weniger gebräuchlich.
Es hilft, sich die Quantile q der exponentialverteilung anzuschauen, hier die “Standardverteilung” mit der Rate 1:
qexp(p = .95, rate = 1)## [1] 2.995732Die Ausgabe sagt uns, dass 95% der Beobachtungen einen Wert kleiner als ca. 3 haben.
Auf unsere Fallstudie bezogen: Wir erwarten, dass 95% der Werte höchstens 3 Einheiten der AV vom vorhergesagten Wert entfernt sind. Das ist ziemlich lax (schwach informativ), denn unsere AV ist z-skaliert.
Um Rechenzeit zu sparen, kann man das Modell auch speichern:
save(post1, file = "post2.rda")
load(file = "post2.rda")5.2.2 Priors
Hier ein Überblick der Priors, die unser Modell post2 verwendet:
prior_summary(post2)## Priors for model 'post2'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 2.1e-16, scale = 2.5)
##
## Coefficients
## ~ normal(location = 0, scale = 0.2)
##
## Auxiliary (sigma)
## ~ exponential(rate = 1)
## ------
## See help('prior_summary.stanreg') for more details5.2.3 Prior-Prädiktiv-Verteilung
5.2.3.1 Berechnen der Prior-Prädiktiv-Verteilung
tic()
post2_prior_pred <- stan_glm(rating_z ~ budget_log10_z,
data = movies3,
prior_PD = TRUE, # Prior-Prädiktiv-Verteilung
prior = normal(0, .2), # beta
refresh = 0) # Nicht so viel Ausgabe
toc()## 2.111 sec elapsed5.2.3.2 Visualisierung der Prior-Prädiktiv-Verteilung
Dazu gehen wir vor wie für Modell 1, post1. Zuerst wandeln wir das Objekt post2 in eine Tabelle mit den Stichproben aus der Post-Verteilung um:
post2_prior_pred_draws <-
post2_prior_pred %>%
as_tibble() %>%
rename(a = `(Intercept)`, # schönere, einfachere Namen
b = budget_log10_z) %>%
slice_sample(n = 100)Und dann visualisieren
movies3 %>%
ggplot() +
geom_point(aes(x = budget_log10_z, y = rating_z)) +
geom_abline(data = post2_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2)
Sieht doch schon viel besser aus. 🏆
5.2.4 Posteriori-Verteilung
5.2.4.1 Kurzfassung
Überblick über die Parameter:
print(post2)## stan_glm
## family: gaussian [identity]
## formula: rating_z ~ budget_log10_z
## observations: 5183
## predictors: 2
## ------
## Median MAD_SD
## (Intercept) 0.0 0.0
## budget_log10_z -0.1 0.0
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 1.0 0.0
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregNur die mittleren Schätzwerte für die Regression:
coef(post2)## (Intercept) budget_log10_z
## 0.0002337387 -0.07868737965.2.4.2 Langfassung
Ausführlicher:
summary(post2)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: rating_z ~ budget_log10_z
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 5183
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## budget_log10_z -0.1 0.0 -0.1 -0.1 -0.1
## sigma 1.0 0.0 1.0 1.0 1.0
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 0.0 0.0 0.0
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3674
## budget_log10_z 0.0 1.0 3871
## sigma 0.0 1.0 3395
## mean_PPD 0.0 1.0 3861
## log-posterior 0.0 1.0 1957
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).5.2.5 Visualisieren von post2
5.2.5.1 Regressionsgerade
plot0_z <- movies3 %>%
ggplot(aes(x = budget_log10_z, y = rating_z)) +
geom_point(alpha = .2)plot1_m2 <-
plot0_z +
geom_abline(intercept = coef(post2)[1],
slope = coef(post2)[2],
color = "blue")
plot1_m2
Erstellen wir eine Tabelle nur mit den Post-Samples:
col_names <- c("a", "b", "sigma")
draws_m2 <-
post2 %>%
as_tibble()
names(draws_m2) <- col_namesDie ersten paar Werte aus der Tabelle mit Post-Samples:
| a | b | sigma |
|---|---|---|
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
| 0.0 | −0.1 | 1.0 |
Und hier die Regressionsgerade mit dem “Unsicherheitsbereich” für die Regressionskoeffizienten.
plot0_z +
geom_abline(data = draws_m2,
aes(intercept = a,
slope = b),
color = "skyblue1",
alpha = .1) +
geom_abline(intercept = coef(post2)[1],
slope = coef(post2)[2],
color = "blue")
5.2.5.2 Verteilung von \(\beta\)
draws_m2 %>%
ggplot(aes(x = b)) +
geom_density()
5.2.5.3 Posterior-Intervalle
S. Hilfe hier
mcmc_areas(post2) +
labs(title = "Posteriori-Verteilung",
caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")
mcmc_intervals(post2,
pars = "budget_log10_z") 
mcmc_areas(post2,
pars = "budget_log10_z") +
labs(title = "Posteriori-Verteilung",
caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")
Das ist das Gleiche wie unser Dichte-Diagramme etwas weiter oben.
5.2.6 Quantile
summary(post2)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: rating_z ~ budget_log10_z
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 5183
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 0.0 0.0 0.0 0.0 0.0
## budget_log10_z -0.1 0.0 -0.1 -0.1 -0.1
## sigma 1.0 0.0 1.0 1.0 1.0
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 0.0 0.0 0.0 0.0 0.0
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3674
## budget_log10_z 0.0 1.0 3871
## sigma 0.0 1.0 3395
## mean_PPD 0.0 1.0 3861
## log-posterior 0.0 1.0 1957
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).Laut dem Modell (post2) liegt der Regressionskoeffizient mit 90% Wahrscheinlichkeit eng um -0.1 herum.
Genauer gesagt: \(90\%PI_b: (-0.13, -0.07)\):
draws_m2 %>%
summarise(b_90pi = quantile(b, probs = c(0.05, .95)))## # A tibble: 2 × 1
## b_90pi
## <dbl>
## 1 -0.102
## 2 -0.0549posterior_interval(post2)## 5% 95%
## (Intercept) -0.02234195 0.02282316
## budget_log10_z -0.10179989 -0.05486384
## sigma 0.98107701 1.013416645.2.7 Wahrscheinlichkeiten für Parameterwerte
5.2.7.1 Positiver Zusammenhang
\(p(b > 0|D)\)
mit “D”, den Daten des Modells.
draws_m2 %>%
count(b > 0)## # A tibble: 1 × 2
## `b > 0` n
## <lgl> <int>
## 1 FALSE 4000Die Wahrscheinlichkeit ist praktisch Null, dass der Zusammenhang positiv ist.
5.2.8 \(R^2\)
tic()
post2_r2 <-
bayes_R2(post2) %>%
as_tibble()
toc()## 1.79 sec elapsedpost2_r2 %>%
ggplot(aes(x=value)) +
geom_density()
post2_r2 %>%
summarise(r2_mean = mean(value),
r2_median = median(value))## # A tibble: 1 × 2
## r2_mean r2_median
## <dbl> <dbl>
## 1 0.00635 0.00619Der Anteil erklärter Varianz ist praktisch Null.
Der Golem blickt es nicht! Das Modell ist schlecht.
5.2.9 Verhältnis von Prior zu Posteriori
Ein interessanter Blick ist der Vergleich der Priori- zu Posteriori-Verteilungen: Um wieviel hat sich unsere Einschätzung (von Priori zu Posteriori) durch die Daten geändert?
Das kann man sich gut so visualisieren lassen:
posterior_vs_prior(post2)
Wir sehen, dass sich die zentrale Tendenz der Parameter (von Prior zu Posterior) kaum geändert hat: Die “Lage” ist fast gleich. Allerdings hat sich die Ungewissheit (des Golems) drastisch verringert. Ob wir seiner Meinung folgen, steht auf einem anderen Blatt: Seine Meinung bezieht sich ausschließlich auf seine kleine Welt. Er weiß nichts von der Welt außer den Informationen, die wir ihm gegeben haben. Wir dagegen vermuten vermutlich, dass noch andere Dinge eine Rolle spielen außer dem Budget, wenn es darum geht, die Beliebtheit eines Films vorherzusagen.
5.2.10 PPV
5.2.10.1 PPV berechnen
Simulieren wir den Erfolg neuer Filme; dabei betrachten wir das Budget von \(10^3\) bis \(10^8\) (6 Werte). Wir ziehen pro Budgetwert 1000 Stichproben aus der PPV.
neue_Filme <-
tibble(
budget_log10_z = 3:8)Warum 3 bis 8? Das sind genau die Werte für budget_log10, die wir in daten Daten haben.
Wie viel delogarithmiertem, also “richtigem” Budget entspricht das?
10^(3:8)## [1] 1e+03 1e+04 1e+05 1e+06 1e+07 1e+08- \(10^3 = 1,000\)
- \(10^4 = 10,000\)
- …
- \(10^8 = 100,000,000\)
ppv_m2 <-
posterior_predict(post2, neue_Filme, draws = 1e3) %>%
as_tibble()
dim(ppv_m2) # Zeilen, Spalten## [1] 1000 6Leider werden unsere Eingabewerte (3 bis 8) zerhauen, stattdessen wird 1 bis 6 zurückgegeben. Kümmern wir uns gleich darum.
Hier ein Blick in die Tabelle ppv_m2:
| 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|
| -0.7370431 | 0.1351947 | -0.4939387 | 0.4804310 | -2.46032163 | 0.05719559 |
| -0.0371196 | -1.3209433 | -0.5672978 | -1.7587757 | -2.71411266 | -0.95442101 |
| 1.5546158 | -0.3867539 | -1.4839460 | -2.1761649 | -0.52137867 | -0.92159574 |
| -0.5827967 | -0.4091292 | -0.4729163 | -0.9540337 | -0.05088593 | -0.27007721 |
| 0.1248318 | 0.2327066 | -0.1887572 | 0.9653534 | -0.67919977 | -0.49252981 |
| -1.7063626 | -1.4473596 | -0.4733521 | 0.4351542 | -1.77418030 | -1.93331043 |
Ändern wir die Spaltennamen von 1,2,…6, in 3,4,…,8 um.
Reine Zahlen werden als Spaltennamen nicht akzeptiert, deswegen wandeln wir noch die Zahl 3 in den Text "3" um, mit as.character():
names(ppv_m2) <- as.character(neue_Filme$budget_log10_z)Vom breiten ins lange Format überführen:
ppv_m2_long <-
ppv_m2 %>%
pivot_longer(everything(),
names_to="budget_log10_z",
values_to="rating")Ein paar Erklärungen zu pivot_wider():
ppv_m2_long %>%
ggplot(aes(x = budget_log10_z,
y = rating)) +
geom_boxplot() 
Tja, keine großen Unterschiede zwischen den Budget-Faktoren (also den Werte von budget_log10). Unser Modell prognostiziert die Daten ganz passabel. Nicht, dass das Ergebnis spektakulär wäre.
5.2.11 90%-Vorhersage-Intervalle
Mit der Funktion predictive_interval kann man sich obige Berechnung sparen, sondern bekommt sie genussfertig nach Hause.
Wir sehen hier die tatsächlichen Rating-Werte pro Budget-Wert, nicht nur \(\mu|b\).
post2_pred <-
predictive_interval(post2,
newdata=neue_Filme)
post2_pred## 5% 95%
## 1 -1.904588 1.415068
## 2 -1.961340 1.326470
## 3 -2.032938 1.282505
## 4 -2.110864 1.141463
## 5 -2.190812 1.140123
## 6 -2.250574 1.057674Wie man sieht, sind die Intervalle sehr groß: Das Modell ist sehr schlecht.
5.3 Modell 3
Schauen wir uns noch mal, was es eigentlich für Prädiktoren gibt bzw. zur Auswahl stehen:
movies %>%
select(-starts_with("r")) %>%
get_summary_stats(type = "robust") %>%
gt()| variable | n | median | iqr |
|---|---|---|---|
| Action | 58788 | 0.000e+00 | 0.000e+00 |
| Animation | 58788 | 0.000e+00 | 0.000e+00 |
| budget | 5215 | 3.000e+06 | 1.475e+07 |
| budget_log10 | 5215 | 6.477e+00 | 1.778e+00 |
| Comedy | 58788 | 0.000e+00 | 1.000e+00 |
| Documentary | 58788 | 0.000e+00 | 0.000e+00 |
| Drama | 58788 | 0.000e+00 | 1.000e+00 |
| length | 58788 | 9.000e+01 | 2.600e+01 |
| Short | 58788 | 0.000e+00 | 0.000e+00 |
| votes | 58788 | 3.000e+01 | 1.010e+02 |
| year | 58788 | 1.983e+03 | 3.900e+01 |
Dabei fällt auf, dass nur (relativ) gesehen sehr wenige Filme eine Info zum Budget haben. Unsere Modell mit budget haben also nur einen geringen Teil der Stichprobe verwendet.
Unser Mäzen schlägt vor, die Variable Action zu untersuchen. Aus seiner Sicht sind Achtionfilme zumeist Mist. Er steht auf Drama, wie er uns berichtet. Aha.
Nehmen wir die Meinung des Mäzens als Ersatz für eine fundierte Literaturrecherche und untersuchen seine Hypothese: Wenn ein Film vom Typ Action ist, dann wird er im Schnitt schlechter beurteilt.
movies4 <-
movies %>%
mutate(rating_z = scale(rating)) %>%
select(rating_z, Action) %>%
drop_na()nrow(movies4) / nrow(movies)## [1] 1Ah! Mit diesen Variablen haben wir kaum Beobachtungen eingebüßt.
post3 <-
stan_glm(rating_z ~ Action,
data = movies4,
refresh = 0)posterior_interval(post3)## 5% 95%
## (Intercept) 0.02872002 0.04271857
## Action -0.47288713 -0.42456944
## sigma 0.98804366 0.99737537Ah! Das sieht ganz gut aus.
Unser Mäzen hatte vielleicht einen ganz guten Riecher. Action scheint einen Einfluss auf das Rating zu haben (zumindest deskriptiv betrachtet).
6 Modellvergleich
modelle <-
tibble(
id = 1:3,
sigma_modell = c(
sigma(post1),
sigma(post2),
sigma(post3)
),
r2_modell = c(
median(bayes_R2(post1)),
median(bayes_R2(post2)),
median(bayes_R2(post3))
)
)
modelle %>%
gt() %>%
fmt_number(where(is.numeric), decimals = 2)| id | sigma_modell | r2_modell |
|---|---|---|
| 1.00 | 1.54 | 0.01 |
| 2.00 | 1.00 | 0.01 |
| 3.00 | 0.99 | 0.01 |
Tja. Das R-Quadrat ist jeweils gering.
7 Fazit
Die Forschungsfrage war, ob das Budget eines Films mit der Bewertung zusammenhängt.
Dazu wurden zwei einfache lineare Modelle berechnet, die sich in ihren Vorannahmen leicht unterschieden.
7.1 Schätzbereiche für \(\beta\)
Beide Modelle fanden konsistent einen schwachen, negativen linearen Zusammenhang \(\beta\) zwischen Budget und Bewertung: Filme mit mehr Budget wurden konsistent schlechter bewertet, laut den beiden Modellen.
Hier sind 90%-PI berichtet:
- Modell 1: [-0.13, -0.07]
- Modell 2: [-0.1, -0.05]
7.2 Medianer Effekt
- Modell 1: -0.1
- Modell 2: -0.08
7.3 Beantwortung der Forschungsfrage
Das Modell ist überzeugt, dass es einen leichten, negativen Zusammenhang gibt. Das Modell schließt aus, dass es keinen Zusammenhang gibt.
7.4 Limitationen
Die Grenzen dieser Analyse sind vielschichtig. Ein Hauptprobleme ist sicherlich, dass wir nur einen Prädiktor verwendet haben, um einen Eindruck von der “Relevanz” des Prädiktors zu bekommen. Dabei ist “Relevanz” wahrscheinlich bereits ein kausales Wort - und ohne ein Kausalmodell können wir über kausale Gegebenheiten keine Schlüsse ziehen. Es bleibt uns der, eher unbefriedigende deskriptive Schluss, dass es einen schwachen, negativen, linearen Zusammenhang zwischen den beiden Größen gibt - laut diesem Modell. Es darf nicht vergessen werden, dass sich der Zusammenhang ändern kann (und es in vielen Fällen auch tun wird), wenn wir das Modell verändern, z.B. weitere Prädiktoren dem Modell hinzufügen.
Ein komfortabler Ausweg ist oder wäre es, sich auf Vorhersagen zu begrenzen. Dafür sind Kausalmodelle nicht zwingend (aber auch nicht verkehrt).
Ein weiterer Kritikpunkt könnte die Wahl unserer Prior sein. Allerdings waren unsere Priors nur schwach informativ, so dass sie wenig Einfluss auf die Posteriori-Verteilung hatten.
8 SessionInfo
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.1.0 (2021-05-18)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Berlin
## date 2021-12-09
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.1.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
## backports 1.2.1 2020-12-09 [1] CRAN (R 4.1.0)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.1.0)
## bayesplot * 1.8.1 2021-06-14 [1] CRAN (R 4.1.0)
## bayestestR * 0.11.0 2021-09-03 [1] CRAN (R 4.1.0)
## boot 1.3-28 2021-05-03 [2] CRAN (R 4.1.0)
## broom 0.7.9 2021-07-27 [1] CRAN (R 4.1.0)
## bslib 0.3.1 2021-10-06 [1] CRAN (R 4.1.0)
## callr 3.7.0 2021-04-20 [1] CRAN (R 4.1.0)
## car 3.0-11 2021-06-27 [2] CRAN (R 4.1.0)
## carData 3.0-4 2020-05-22 [2] CRAN (R 4.1.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
## checkmate 2.0.0 2020-02-06 [2] CRAN (R 4.1.0)
## cli 3.1.0 2021-10-27 [1] CRAN (R 4.1.0)
## codetools 0.2-18 2020-11-04 [2] CRAN (R 4.1.0)
## colorspace 2.0-2 2021-06-24 [1] CRAN (R 4.1.0)
## colourpicker 1.1.0 2020-09-14 [2] CRAN (R 4.1.0)
## crayon 1.4.2 2021-10-29 [1] CRAN (R 4.1.0)
## crosstalk 1.1.1 2021-01-12 [2] CRAN (R 4.1.0)
## curl 4.3.2 2021-06-23 [1] CRAN (R 4.1.0)
## data.table 1.14.2 2021-09-27 [1] CRAN (R 4.1.0)
## datawizard 0.2.1 2021-10-04 [1] CRAN (R 4.1.0)
## DBI 1.1.1 2021-01-15 [1] CRAN (R 4.1.0)
## dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.1.0)
## digest 0.6.28 2021-09-23 [1] CRAN (R 4.1.0)
## dplyr * 1.0.7 2021-06-18 [1] CRAN (R 4.1.0)
## DT 0.19.1 2021-10-15 [1] Github (rstudio/DT@e223814)
## dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 4.1.0)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.1.0)
## fansi 0.5.0 2021-05-25 [1] CRAN (R 4.1.0)
## farver 2.1.0 2021-02-28 [1] CRAN (R 4.1.0)
## fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.1.0)
## forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.1.0)
## foreign 0.8-81 2020-12-22 [2] CRAN (R 4.1.0)
## fs 1.5.0 2020-07-31 [1] CRAN (R 4.1.0)
## generics 0.1.1 2021-10-25 [1] CRAN (R 4.1.0)
## ggplot2 * 3.3.5 2021-06-25 [1] CRAN (R 4.1.0)
## ggplot2movies * 0.0.1 2015-08-25 [1] CRAN (R 4.1.0)
## ggridges 0.5.3 2021-01-08 [1] CRAN (R 4.1.0)
## glue 1.5.1 2021-11-30 [1] CRAN (R 4.1.0)
## gridExtra 2.3 2017-09-09 [2] CRAN (R 4.1.0)
## gt * 0.3.1 2021-08-07 [1] CRAN (R 4.1.0)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.1.0)
## gtools 3.9.2 2021-06-06 [2] CRAN (R 4.1.0)
## haven 2.4.1 2021-04-23 [1] CRAN (R 4.1.0)
## highr 0.9 2021-04-16 [1] CRAN (R 4.1.0)
## hms 1.1.1 2021-09-26 [1] CRAN (R 4.1.0)
## htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.1.0)
## htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.1.0)
## httpuv 1.6.3 2021-09-09 [1] CRAN (R 4.1.0)
## httr 1.4.2 2020-07-20 [1] CRAN (R 4.1.0)
## igraph 1.2.6 2020-10-06 [1] CRAN (R 4.1.0)
## inline 0.3.19 2021-05-31 [1] CRAN (R 4.1.0)
## insight 0.14.5 2021-10-16 [1] CRAN (R 4.1.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
## jsonlite 1.7.2 2020-12-09 [1] CRAN (R 4.1.0)
## knitr 1.36 2021-09-29 [1] CRAN (R 4.1.0)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.1.0)
## lattice 0.20-44 2021-05-02 [2] CRAN (R 4.1.0)
## lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.1.0)
## lme4 1.1-27 2021-05-15 [1] CRAN (R 4.1.0)
## loo 2.4.1 2020-12-09 [1] CRAN (R 4.1.0)
## lubridate 1.7.10 2021-02-26 [1] CRAN (R 4.1.0)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.1.0)
## markdown 1.1 2019-08-07 [1] CRAN (R 4.1.0)
## MASS 7.3-54 2021-05-03 [2] CRAN (R 4.1.0)
## Matrix 1.3-4 2021-06-01 [2] CRAN (R 4.1.0)
## matrixStats 0.61.0 2021-09-17 [1] CRAN (R 4.1.0)
## mgcv 1.8-36 2021-06-01 [2] CRAN (R 4.1.0)
## mime 0.12 2021-09-28 [1] CRAN (R 4.1.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.0)
## minqa 1.2.4 2014-10-09 [1] CRAN (R 4.1.0)
## modelr 0.1.8 2020-05-19 [1] CRAN (R 4.1.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
## nlme 3.1-152 2021-02-04 [2] CRAN (R 4.1.0)
## nloptr 1.2.2.2 2020-07-02 [1] CRAN (R 4.1.0)
## openxlsx 4.2.3 2020-10-27 [1] CRAN (R 4.1.0)
## pillar 1.6.4 2021-10-18 [1] CRAN (R 4.1.0)
## pkgbuild 1.2.0 2020-12-15 [2] CRAN (R 4.1.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
## plyr 1.8.6 2020-03-03 [1] CRAN (R 4.1.0)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.1.0)
## processx 3.5.2 2021-04-30 [1] CRAN (R 4.1.0)
## promises 1.2.0.1 2021-02-11 [2] CRAN (R 4.1.0)
## ps 1.6.0 2021-02-28 [1] CRAN (R 4.1.0)
## purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.1.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
## Rcpp * 1.0.7 2021-07-07 [1] CRAN (R 4.1.0)
## RcppParallel 5.1.4 2021-05-04 [1] CRAN (R 4.1.0)
## readr * 2.0.0 2021-07-20 [1] CRAN (R 4.1.0)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 4.1.0)
## reprex 2.0.0 2021-04-02 [1] CRAN (R 4.1.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.1.0)
## rio 0.5.26 2021-03-01 [1] CRAN (R 4.1.0)
## rlang 0.4.12 2021-10-18 [1] CRAN (R 4.1.0)
## rmarkdown 2.11 2021-09-14 [1] CRAN (R 4.1.0)
## rsconnect 0.8.25 2021-11-19 [1] CRAN (R 4.1.0)
## rstan 2.21.2 2020-07-27 [1] CRAN (R 4.1.0)
## rstanarm * 2.21.1 2020-07-20 [1] CRAN (R 4.1.0)
## rstantools 2.1.1 2020-07-06 [1] CRAN (R 4.1.0)
## rstatix * 0.7.0 2021-02-13 [1] CRAN (R 4.1.0)
## rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.1.0)
## rvest 1.0.0 2021-03-09 [1] CRAN (R 4.1.0)
## sass 0.4.0 2021-05-12 [1] CRAN (R 4.1.0)
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.1.0)
## sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.1.0)
## shiny 1.7.1 2021-10-02 [1] CRAN (R 4.1.0)
## shinyjs 2.0.0 2020-09-09 [2] CRAN (R 4.1.0)
## shinystan 2.5.0 2018-05-01 [1] CRAN (R 4.1.0)
## shinythemes 1.2.0 2021-01-25 [1] CRAN (R 4.1.0)
## StanHeaders 2.21.0-7 2020-12-17 [1] CRAN (R 4.1.0)
## stringi 1.7.5 2021-10-04 [1] CRAN (R 4.1.0)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.1.0)
## survival 3.2-11 2021-04-26 [2] CRAN (R 4.1.0)
## threejs 0.3.3 2020-01-21 [1] CRAN (R 4.1.0)
## tibble * 3.1.6 2021-11-07 [1] CRAN (R 4.1.0)
## tictoc * 1.0.1 2021-04-19 [2] CRAN (R 4.1.0)
## tidyr * 1.1.4 2021-09-27 [1] CRAN (R 4.1.0)
## tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.1.0)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.1.0)
## tzdb 0.1.2 2021-07-20 [2] CRAN (R 4.1.0)
## utf8 1.2.2 2021-07-24 [1] CRAN (R 4.1.0)
## V8 3.4.2 2021-05-01 [1] CRAN (R 4.1.0)
## vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.1.0)
## withr 2.4.2 2021-04-18 [1] CRAN (R 4.1.0)
## xfun 0.28 2021-11-04 [1] CRAN (R 4.1.0)
## xml2 1.3.2 2020-04-23 [1] CRAN (R 4.1.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
## xts 0.12.1 2020-09-09 [1] CRAN (R 4.1.0)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.1.0)
## zip 2.2.0 2021-05-31 [2] CRAN (R 4.1.0)
## zoo 1.8-9 2021-03-09 [1] CRAN (R 4.1.0)
##
## [1] /Users/sebastiansaueruser/Library/R/x86_64/4.1/library
## [2] /Library/Frameworks/R.framework/Versions/4.1/Resources/library