library(tidyverse)
library(rstanarm)stan_glm_parameterzahl
bayes
regression
parameters
Exercise
Berechnet man eine Posteriori-Verteilung mit stan_glm(), so kann man entweder die schwach informativen Prioriwerte der Standardeinstellung verwenden, oder selber Prioriwerte definieren.
Betrachten Sie dazu dieses Modell:
stan_glm(price ~ cut, data = diamonds,
prior = normal(location = c(100, 100, 100, 100),
scale = c(100, 100, 100, 100)),
prior_aux = exponential(1),
prior_intercept = normal(3000, 500))Wie viele Parameter gibt es in diesem Modell?
Hinweise:
- Geben Sie nur eine (ganze) Zahl ein.
Solution
Grundsätzlich hat ein Regressionsmodell die folgenden Parameter:
- einen Parameter für den Intercept, \(\beta_0\)
- pro UV ein weiterer Parameter, \(\beta_1, \beta_2, \ldots\)
- für sigma (\(\sigma\)) noch ein zusätzlicher Parameter
Zu beachten ist aber, dass bei einer nominalen Variablen mit zwei Stufen nur ein Regressionsgewicht (\(\beta_1\)) berechnet wird. Allgemein gilt bei nominalen also, dass bei \(k\) Stufen nur \(k-1\) Regressionsgewichte berechnet werden.
Im vorliegenden Fall hat die Variable cut 5 Stufen, also werden 4 Regressiongewichte berechnet.
Fazit: In Summe werden also 6 Parameter berechnet.
Berechnet man das Modell, so kann man sich auch Infos über die Prioris ausgeben lassen:
m1 <- stan_glm(price ~ cut, data = diamonds,
prior = normal(location = c(100, 100, 100, 100),
scale = c(100, 100, 100, 100)),
prior_intercept = normal(3000, 500),
prior_aux = exponential(1),
refresh = 0)
prior_summary(m1)Priors for model 'm1'
------
Intercept (after predictors centered)
~ normal(location = 3000, scale = 500)
Coefficients
~ normal(location = [100,100,100,...], scale = [100,100,100,...])
Auxiliary (sigma)
~ exponential(rate = 1)
------
See help('prior_summary.stanreg') for more detailsWie man sieht, wird für die Streuung im Standard eine Exponentialverteilung verwendet von stan_glm(). Gibt man also nicht an - wie im Beispiel m1 oben, so wird stan_glm() für die Streuung, d.h. prior_aux eine Exponentialverteilung verwenden. Zu beachten ist, dass stan_glm() ein automatische Skalierung vornimmt.
S. hier für weitere Erläuterung.
Möchte man den Prior für die Streuung direkt ansprechen, so kann man das so formulieren:
m2 <- stan_glm(price ~ cut, data = diamonds,
prior = normal(location = c(100, 100, 100, 100),
scale = c(100, 100, 100, 100)),
prior_intercept = normal(3000, 500),
prior_aux = exponential(1),
refresh = 0)
prior_summary(m1)Priors for model 'm1'
------
Intercept (after predictors centered)
~ normal(location = 3000, scale = 500)
Coefficients
~ normal(location = [100,100,100,...], scale = [100,100,100,...])
Auxiliary (sigma)
~ exponential(rate = 1)
------
See help('prior_summary.stanreg') for more detailsZu beachten ist beim selber definieren der Prioris, dass dann keine Auto-Skalierung von stan_glm() vorgenommen wird, es sei denn, man weist es explizit an:
m3 <- stan_glm(price ~ cut, data = diamonds,
prior = normal(location = c(100, 100, 100, 100),
scale = c(100, 100, 100, 100),
autoscale = TRUE),
prior_intercept = normal(3000, 500, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chain = 1, # nur 1 mal Stichproben ziehen, um Zeit zu sparen
refresh = 0)
prior_summary(m3)Priors for model 'm3'
------
Intercept (after predictors centered)
Specified prior:
~ normal(location = 3000, scale = 500)
Adjusted prior:
~ normal(location = 3000, scale = 2e+06)
Coefficients
Specified prior:
~ normal(location = [100,100,100,...], scale = [100,100,100,...])
Adjusted prior:
~ normal(location = [100,100,100,...], scale = [1129833.17, 868199.02, 936606.47,...])
Auxiliary (sigma)
Specified prior:
~ exponential(rate = 1)
Adjusted prior:
~ exponential(rate = 0.00025)
------
See help('prior_summary.stanreg') for more detailsGrundsätzlich ist es nützlich für die numerische Stabilität, dass die Zahlen (hier die Parameterwerte) etwa die gleiche Größenordnung haben, am besten um die 0-1 herum. Daher bietet sich oft eine z-Standardisierung an.
Unabhängig von der der Art der Parameter ist die Anzahl immer gleich.
Die Anzahl der geschätzten Parameter werden im Modell-Summary unter Estimates gezeigt:
summary(m1)
Model Info:
function: stan_glm
family: gaussian [identity]
formula: price ~ cut
algorithm: sampling
sample: 4000 (posterior sample size)
priors: see help('prior_summary')
observations: 53940
predictors: 5
Estimates:
mean sd 10% 50% 90%
(Intercept) 4027.5 21.7 4000.1 4027.7 4055.3
cut.L -249.1 51.5 -315.9 -249.4 -182.0
cut.Q -283.4 45.5 -340.8 -283.8 -223.7
cut.C -581.1 43.7 -637.0 -582.1 -523.9
cut^4 -265.7 37.6 -313.3 -265.6 -217.4
sigma 3830.2 11.3 3815.7 3830.1 3844.6
Fit Diagnostics:
mean sd 10% 50% 90%
mean_PPD 3931.8 23.1 3902.7 3932.1 3961.4
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.4 1.0 2750
cut.L 1.0 1.0 2604
cut.Q 1.0 1.0 2246
cut.C 0.9 1.0 2536
cut^4 0.7 1.0 3024
sigma 0.1 1.0 6932
mean_PPD 0.4 1.0 3604
log-posterior 0.0 1.0 1701
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).Das sind:
- 1 Intercept (Achsenabschnitt) -
prior_intercept - 4 Gruppen (zusätzlich zur Referenzgruppe, die mit dem Achsenabschnitt dargestellt ist) -
prior_normal - 1 Sigma (Ungewissheit “innerhalb des Modells”) -
prior_aux
Categories:
~