In unserem Datensatz zur Fruchtbarkeit in der Schweiz gibt es weitere Variablen:

data(swiss)
names(swiss)

## [1] "Fertility"        "Agriculture"      "Examination"      "Education"       
## [5] "Catholic"         "Infant.Mortality"

Fertility, der Fertilitätsindex, ist unsere Zielvariable. Die restlichen fünf können wir als Kovariablen benutzen.

##   Agriculture     Examination      Education        Catholic      
##  Min.   : 1.20   Min.   : 3.00   Min.   : 1.00   Min.   :  2.150  
##  1st Qu.:35.90   1st Qu.:12.00   1st Qu.: 6.00   1st Qu.:  5.195  
##  Median :54.10   Median :16.00   Median : 8.00   Median : 15.140  
##  Mean   :50.66   Mean   :16.49   Mean   :10.98   Mean   : 41.144  
##  3rd Qu.:67.65   3rd Qu.:22.00   3rd Qu.:12.00   3rd Qu.: 93.125  
##  Max.   :89.70   Max.   :37.00   Max.   :53.00   Max.   :100.000  
##  Infant.Mortality
##  Min.   :10.80   
##  1st Qu.:18.15   
##  Median :20.00   
##  Mean   :19.94   
##  3rd Qu.:21.70   
##  Max.   :26.60

Die Kovariablen sind:

• Agriculture: Anteil der Arbeitsplätze (von Männern) in der Landwirtschaft
• Examination: Anteil Wehrpflichtiger, die bei der Musterung die höchste Note erhalten
• Education: Anteil an Wehrpflichtigen mit Bildung über die Volksschule hinaus
• Catholic: Anteil Katholiken
• Infant.Mortality: Zahl Lebendgeborener, die im ersten Lebensjahr starben

Wir erweiteren unser lineares Regressionsmodell für \(y\) mit mehreren Kovariablen \(x_1, \ldots,x_p\):

\[ y_i=\alpha+\beta_1 x_{1i}+\beta_2x_{2j}+\ldots\beta_px_{pi}+\epsilon_i, \quad i=1,\ldots,n\\ \]

Die Annahmen an den Fehler \(\epsilon_i\) sind identisch wie zuvor:

• \(\epsilon_i\sim\)N\((0,\sigma^2)\)
• Cov\((\epsilon_i,\epsilon_j)=0\)

Mit \(\beta_0:=\alpha\) und \(x_0=(1,1,\ldots,1)\) schreiben wir:

\[ y_i=\sum_{j=0}^p\beta_j x_{ji}+\epsilon_i, \quad i=1,\ldots,n\\ \]

Datendichte

Aus dem linearen Modell \(y_i=\sum_{j=0}^p\beta_j x_{ji}+\epsilon_i\) und der Annahme \(\epsilon_i\sim\)N\((0,\sigma^2)\) ergibt sich eine Verteilung für \(y\)

\[ y_i|\beta, \sigma^2 \sim \text{N}\left(\sum_{j=0}^p\beta_j x_{ji},\sigma^2\right) \]

Dabei ist \(\beta\) nun ein Vektor. Davon abgesehen hat sich im Vergleich zur einfachen Regression kaum etwas geändert, nur der Erwartungswert, der sich aus gewichteter Summe der Kovariablen angenommen wird (für \(p=1\) ist das Modell natürlich identisch zur einfachen Regression).

Die Datendichte aller \(y=(y_1,\ldots,y_n)\) gegeben den Parametern \(\theta=(\beta_0,\beta_1,\ldots,\beta_p,\sigma^2)\) ist damit:

\[ \begin{aligned} f(y|\theta) & = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{1}{2\sigma^2}\left(y_i-\sum_{j=0}^p\beta_j x_{ji}\right)^2\right)\\ & = \frac{1}{(2\pi\sigma^2)^{n/2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-\sum_{j=0}^p\beta_j x_{ji}\right)^2\right) \end{aligned} \]

Priori der Regressionsparameter

Unbekannte Parameter unseres Modells sind \(\theta=(\beta_0,\beta_1,\ldots,\beta_p,\sigma^2)\). Für diese müssen wir Prioris festlegen. Dafür übernehmen wir die Prioris aus der einfachen Regression:

• \(\beta_j\sim\)N\((m_j,v_j^2)\) für \(j=0,\ldots,p\).

a priori seien die \(\beta_j\) alle voneinander unabhängig.

Bemerkung: Man kann hier auch Abhängigkeiten in die Priori aufnehmen, wenn entsprechendes Vorwissen vorliegt.

Priori der Varianz

Analog benutzen wir für die Varianz \(\sigma^2\) die (semi-)konjugierte Priori:

\[ \sigma^2 \sim IG(a,b), \]

Gemeinsame Priori

Zusätzlich gehen wir a priori davon aus, das \(\sigma^2\) von \(\alpha\) und \(\beta\) stochastisch unabhängig ist. Damit erhalten wir die gemeinsame Priori-Dichte aller Parameter:

\[ p(\beta, \sigma^2)=\prod_{j=0}^pp(\beta_j)\cdot p(\sigma^2) \]

Posteriori

Leiten wir nun die Posteriori der Parameter gegeben der Daten her. Es gilt mit Daten \(y=(y_1,\ldots,y_n)\) und Parametern \(\theta=(\beta_0,\beta_1,\ldots,\beta_p,\sigma^2)\):

\[\begin{eqnarray*} p(\theta|y)&\propto&f(y|\theta)p(\theta)\\ &=& \prod_{i=1}^n f(y_i|\theta)\prod_{j=0}^pp(\beta_j)p(\sigma^2)\\ &\propto & \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-\sum_{j=0}^p\beta_j x_{ji}\right)^2\right)\cdot\\ &&\cdot v_j^{-1} \exp \left(-\frac{1}{2v_j^2} (\beta_j-m_j)^2\right)\\ &&\cdot \sigma^{-2(a+1)}\exp\left(b\sigma^{-2}\right) \end{eqnarray*}\]

Vollständig bedingte Posterioris

Aus der Posterior leiten wir uns nun wieder die vollständig bedingten Posterioris (full conditionals) her. Zur Erinnerung:

\[ f(\theta_i|\theta_{-i})=\frac{f(\theta)}{f(\theta_{-i})}\propto f(\theta) \]

Full conditional von \(\beta_j\)

Wir benutzten wieder die Notation \(\beta_{-j}=(\beta_0,\ldots,\beta_{j-1},\beta_{j+1},\ldots,\beta_p)\) und suchen die Terme aus der Posteriori, die von \(\beta_j\) abhängen:

\[\begin{eqnarray*} p(\beta_j|\beta_{-j},\sigma^2;y) &\propto& \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-\sum_{k=0}^p\beta_k x_{ki}\right)^2\right) \prod_{j=0}^p \exp \left(-\frac{1}{2v_j^2} (\beta_j-m_j)^2\right)\\ & \propto& \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n \left(\beta_jx_{ij}-\left(y_i-\sum_{k\not=j}\beta_k x_{ki}\right)\right)^2\right) \exp \left(-\frac{1}{2v_j^2} (\beta_j-m_j)^2\right) \end{eqnarray*}\]

Wir haben in der Exponentialfunktiong wieder die Summe von zwei quadratischen Formen. Das heißt, wir bekommen den Kern einer Normalverteilung. Wie bei der einfachen Regression kann man weiter machen und erhält:

• \(\beta_j|\beta_{-j},\sigma^2;y\) ist normalverteilt mit

• Inverser Varianz \(\tilde{v}_{j}^2=v_j^{-2}+\sigma^{-2}\sum_{i=1}^nx_{ji}^2\) und

• Erwartungwert \(\tilde{m}_j=\tilde{v}_j^{-2}\left(v_{\beta}^{-2}m_\beta +\sigma^{-2}\sum_{i=1}^n\left(x_{ji}\left(y_i-\sum_{k\not=j}\beta_kx_{ki}\right)\right)\right)\).

Das Ergebnis für die Full conditional unterscheidet sich kaum von der full conditional in der einfachen Regression. Bis auf die Notation wurde \(\alpha\) ersetzt durch \(\sum_{k\not=j}\beta_kx_{ki}\), also dem Teil des Modells, der durch die anderen Kovariablen erklärt wird.

Full conditional von \(\sigma^2\)

Die Full Conditional von \(\sigma^2\) ist ebenfalls praktisch identisch mit der in der einfachen Regression:

• \(\tilde{a}=a+n/2\)
• \(\tilde{b}=b+0.5\sum_{i=1}^n\left(y_i-\sum_{j=0}^p\beta_j x_{ji}\right)^2\)

• Wir sehen, dass sich die Full conditionals immer wieder wiederholen.
• Wir könnten für die Implementierung immer wieder Code-Teil wiederverwenden.

Hier benutzen wir rJAGS.

Modellierung

In der Modellierungssprache von JAGS werden die verschiedenen Verteilungen beschrieben:

• Die Datendichte für jede Beobachtung \(i=1,\ldots,n\), hier Fertility\(\sim\)N\((\mu_i,\tau)\), wobei \(\tau\) die Präzision (inverse Varianz) ist
• Die Prioris
  • \(\beta_j\sim N(0,10^6)\)
  • \(\tau\sim Ga(0.001,0.001)\)

Dazu kommen deterministische Umrechnungen, wie

• \(\mu_i=\sum_j\beta_jx_{ji}\)
• Varianz \(\sigma^2=1/tau\)

linmodel <- "model{
    for (i in 1:n) {
        Fertility[i] ~ dnorm(mu[i], tau)
        mu[i] <- beta[1] + beta[2] * Agriculture[i] + beta[3] * Examination[i] + beta[4] * Education[i]  + beta[5] * Catholic[i] + beta[6] * Infant.Mortality[i]
    }
    ## Prioris
    for(j in 1:p){
    beta[j] ~ dnorm(0, 0.001)
    }
    tau ~ dgamma(0.001, 0.001)
    sigma2 <- 1/tau
}"

In rjags muss das Modell als Text (oder als Datei) übergeben werden. Daher alles in Anführungszeichen.

Die Daten werden als Liste übergeben. Neben Ziel- und Kovariable sind dies hier auch \(n\) (Anzahl Beobachtungen) und \(p\) (Anzahl Kovariablen plus Intercept).

daten<-as.list(swiss)
daten$p=dim(swiss)[2]
daten$n=dim(swiss)[1]

Startwerte kann man übergeben, diese kann JAGS aber auch selbst erzeugen.

startwerte <-
    list(beta = c(60, 0.2, 0, 0, 0, 0), tau = 1/144)

Komplilieren

JAGS kompililiert aus dem Modell und den Daten den passenden Algorithmus

model <- jags.model(file = textConnection(linmodel), data = daten, inits = startwerte)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 47
##    Unobserved stochastic nodes: 7
##    Total graph size: 513
## 
## Initializing model

Durchführen

Wir führen 10000 Iterationen als burn-in durch.

update(model, 10000, progress.bar="none"); 

Danach ziehen wir 20000-mal aus den Parametern

samp <- jags.samples(model, 
        variable.names=c("beta","sigma2"), 
        n.iter=20000, progress.bar="none")

Schätzer

summary(samp$beta,mean)
summary(samp$beta,median)
summary(samp$beta,quantile,0.025)
summary(samp$beta,quantile,0.975)

Interpretation der Ergebnisse

• Die “Zahl der Lebendgeborenen, die im ersten Lebensjahr starben” (Infant.Mortality) hat einen positiven Einfluss auf die Fruchtbarkeit.
• Der “Anteil der Katholiken” (Catholic) hat einen leicht positiven Einfluss auf die Fruchtbarkeit.
• “Anteil an Wehrpflichtigen mit Bildung über die Volksschule hinaus” (Education) hat einen negativen Einfluss auf die Fruchtbarkeit.
• Der Erwartungswert des Einflusses von Examination (“Anteil Wehrpflichtiger, die bei der Musterung die höchste Note erhalten haben”) ist negativ. Das 95%-Quantil des Einflusses enthält dagegen die Null, dementsprechend können wir hier keinen Effekt auf die Fruchtbarkeit feststellen.