stan_glm_parameterzahl

Published

December 13, 2022

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, $β_{0}$
pro UV ein weiterer Parameter, $β_{1}, β_{2}, \dots$
für sigma ( $σ$ ) noch ein zusätzlicher Parameter

Zu beachten ist aber, dass bei einer nominalen Variablen mit zwei Stufen nur ein Regressionsgewicht ( $β_{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.

library(tidyverse)
library(rstanarm)

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 details

Wie 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 details

Zu 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 details

Grundsä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: