Bayes (Hierarcical Bayesian model)

1. JOHDANTO

1.1 BAYESMENETELMÄ

Bayes menetelmä perustuu bayesin sääntöön. Seuraavat kaavat löytyvät kirjasta Bayesian Data Analysis second edition sivuilta 7-8.

Bayesin sääntö on muotoa.

p(θ,y) = p(θ)p(y|θ)

Tässä p(θ) on priori jakauma ja p(y|θ) on datasta saatu jakauma. Ehdollistamalla θ:taa y:llä, joka on tunnettu, päästään muotoon.

p(θ|y) = p(θ)p(y|θ) / p(y)

Näin ehdollistamalla saadaan p(θ|y), joka on siis selkokielellä prioritodennäköisyys ehdollistettuna dataan. Tätä kutsutaan posterioritiheydeksi. Näillä kaavoilla saadaan siis tarkennettua ennakko-oletuksia datan tiheysfunktiosta (priorijakauma) käyttäen dataa itseään hyväksi. Lopputuloksena on siis päivitetty tiheysfunktio jota kutsutaan posterioritiheydeksi.

Muistisääntö:
Posteriorijakauma ~ priorijakauma x datasta estimoitu jakauma

1.2 HIERARKINEN BAYES

Hierarkisen bayesin idea on helpoin ymmärtää konkreettisen esimerkin kautta. Meidän tapauksessa havaitaan pitoisuuksia eri alueittain, joiden voidaan olettaa olevan alueittain samankaltaisia. Oletetaan, että pitoisuudet alueittain ovat otoksia alueiden yhteisestä pitoisuus jakaumasta. Näin ollen voimme käyttää havaittuja pitoisuuksia yli alueiden estimoitaessa tietyn alueen pitoisuuksia. Tämä tapahtuu ehdollistamalla havaittavissa olevia arvoja tiettyihin parametreihin. Näille parametreille oletetaan todennäköisyydet uusien parametrien avulla, joita kutsutaan hyperparametreiksi.

Kirjassa Bayesian Data Analysis second edition sivuilla 117-121 löytyy kattava selvennys ideasta ja sivulla 119 on kuva hierarkisesta mallista.

Käyttötarkoitus meillä on, että pyritään löytämään keinoja mallintaa tilastollisin keinoin pitoisuuksia kaloista. Tavoitteena saada aikaan malli, jonka avulla pienestäkin datasta saataisiin suhteellisen pitävä arvio koko alueen kalapopulaatiolle.

2. MENETELMÄT JA AINEISTOT

2.1 MENETELMÄT

Käytettävä menetelmä on hierarkkinen bayes. Teoria kirjasta Bayesian data analysis second edition. Työkaluna WinBUGS ja R.

http://www.mrc-bsu.cam.ac.uk/bugs/

Tästä osoitteesta voi hakea päivityksen ja lisenssin. Päivitysohjeet ja lisenssiohjeet kattavat, eli tässä ei pitäisi tulla ongelmia. Tekstin vain kopioi WinBUGS uuteen tiedostoon ja tools valikosta decode kuten ohjeissa mainitaan .

Minkä nimisiä ovat muuttujat mallissa, taulukko (muuttujan nimi ja parametrit)
Mallin graafinen kuva (vaikka sama kuva kuin WinBUGS:ssa) muuttujien vuorovaikutussuhteista

Taulukko1. Mallin muuttujien nimet ja parametrit

Kuva 1. Mallin graafinen kuva WinBUGS:ssa

2.2. AINEISTO

Datassa on silakan kokoja, ikiä sekä arvoja dioksiini-pitoisuuksille (pg/g WHOPCDD/F)

Ajatellaan pitoisuutta kalassa, hieman samaan tyyliin, kuin rotta esimerkissä on ajateltu painoa, eli malli on hyvin samankaltainen.

Muuttujat:

Rotan paino = Silakan pitoisuus
Rotan ikä = Silakan ikä
Monesko rotta = Mistä pyydetty

Tämän jälkeen WinBUGS:lla voi pyörittää rotta mallia silakka-datalla. Näin saadaan mielenkiintoiset muuttujat ja jakaumat esille.

Rotta esimerkissä siis estimoidaan jakaumat painolla jokaiselle rotalle erikseen. Meillä estimoidaan jakaumat jokaiselle pitoisuudelle alueittain.

3. TULOKSET

3.1 KUVAILEVIA TILASTOLLISIA TULOKSIA

Kuva 2. Silakan ikä pitoisuuden suhteen.

Kuva 3. Silakan paino pituuden suhteen.

Lineaarisen regressiomallin tuloksia, joissa pitoisuutta selitetään iällä ja painolla.

Call:
lm(formula = pitoisuus ~ ikä + paino, data = data)

Residuals:
Min 1Q Median 3Q Max
-7.9395 -1.3808 -0.1211 1.2573 6.4274

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.68568 0.78077 -2.159 0.0363 *
ikä 1.66347 0.14978 11.106 2.39e-14 ***
paino -0.03660 0.01205 -3.038 0.0040 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.575 on 44 degrees of freedom
Multiple R-squared: 0.7495, Adjusted R-squared: 0.7381
F-statistic: 65.82 on 2 and 44 DF, p-value: 5.952e-14

Huomataan siis että ikä ja paino merkitseviä selittäjiä.

3.2 WinBUGS ohjelmistolla laskettuja tuloksia.

Seuraavissa taulukoissa 2 ja 3 esitetään hierarkinen bayes malli silakka-datalle. Huomaamisen arvoista on, että Y[1,1] on ensimmäisen alueen 1 vuotiaan silakan pitoisuus jne. Ensimmäinen indeksi on siis paikka indeksi ja toinen on ikä indeksi. Osa arvoista pienempiä kuin nolla, mutta Y:lle voi asettaa alarajaksi nolla, jolloin tästä ongelmasta päästään. Liitteet kulkukuva ja tiheydet liittyvät tähän ja aukeavat bugs:ssa. Taulukossa 3 iät on jaettu neljään ryhmään kolmen vuoden jaksoissa. 1-3, 4-6 jne. Sekä näitä vastaavat tulokset. Huomataan ettei miinuksia ole keskiarvoissa, joskin luottamusväleistä niitä löytyy. Taulukon 2 estimaatit ovat vain puuttuville arvoille. Mediaanit ovat melko hyviä arvioita, mutta luottamusvälit leviävät kohtuuttoman isoiksi, johtuen liian suuresta määrästä puuttuvia havaintoja.

Taulukko 2. WinBUGSin tuloslistaus silakkadatalle yhden vuoden ikä-intervallilla.

Taulukko 3. WinBUGSin tuloslistaus silakkadatalle kolmen vuoden ikä-intervallilla.

3.2. Koodit ja niiden ajaminen

Liitteenä tulee koodit model1 , noodle_1, ja data ja data_1 joista data:ssa on 1-4 ikäluokat ja data_1 löytyy kaikille iän arvoille vastaava, sekä inits jossa alkuarvoja, mikäli niitä ei tahdota simuloida.

Kaikki koodit aukeavat WinBUGS ohjelmalla. Ja ovat vapaasti päivitettävissä. Ajot on tehty mallilla model1 ja datalla data ja data_1 ja todettu toimiviksi alkuarvoilla inits. Vastaavanlainen malli on esimerkeissä joka löytyy helpin alta examples vol1 Rats.

rajan asettaminen tapahtuu laittamalla Y:n koodin eli jakauman ja arvojen perään symbolin I(0,)

Ajaminen:

Mikäli ongelmia niin soita/tule kysymään. Tarvitset kolme ikkunaa auki. Modelin alta model specification, Inferencen alta samples ja modelin alta update, niin pääsee tekemään samat ajot kuin mitä on tässä tehty.

Koodit löytyy: N:\Huippuyksikko\Tutkimus\R55_PMvsDioxin\MeHgvsOmega-3\mallit\Bayes\Aiemmat ajot ja pähkäilyt

5. JOHTOPÄÄTÖKSET

Bugs ohjelmistolla saatuja tuloksia voi verrata silakka-dataan, ja sen pohjalta tehdä päätelmiä.

ESIMERKKI 1.

Hierarkinen rotta esimerkki, joka löytyy winbugsin helpistä. Tähän malliin perustuu suoraan meidän käyttämät mallit.

Rats: a normal hierarchical model

This example is taken from section 6 of Gelfand et al (1990), and concerns 30 young rats whose weights were measured weekly for five weeks. Part of the data is shown below, where Yij is the weight of the ith rat measured at age xj.

Weights Yij of rat i on day xj
xj = 8 15 22 29 36
__________________________________
Rat 1 151 199 246 283 320
Rat 2 145 199 249 293 354
.......
Rat 30 153 200 244 286 324

A plot of the 30 growth curves suggests some evidence of downward curvature.

The model is essentially a random effects linear growth curve

Yij ~ Normal(i + i(xj - xbar), c)

i ~ Normal(c, )

i ~ Normal(c, )

where xbar = 22, and  represents the precision (1/variance) of a normal distribution. We note the absence of a parameter representing correlation between i and i unlike in Gelfand et al 1990. However, see the Birats example in Volume 2 which does explicitly model the covariance between i and i. For now, we standardise the xj's around their mean to reduce dependence between i and i in their likelihood: in fact for the full balanced data, complete independence is achieved. (Note that, in general, prior independence does not force the posterior distributions to be independent).

c , , c , , c are given independent "noninformative priors, with two alternatives considered for  and : prior 1 is uniform on the scale of the standard deviations  = 1/sqrt() and  = 1/sqrt(), and prior 2 is a gamma(0.001, 0.001) on the precisions  and . Interest particularly focuses on the intercept at zero time (birth), denoted 0 = c - c xbar.

Graphical model for rats example (using prior 1):

BUGS language for rats example:

Note the use of a very flat but conjugate prior for the population effects: a locally uniform prior could also have been used.

Data list(x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, N = 30, T = 5,
Y = structure(
.Data = c(151, 199, 246, 283, 320,
145, 199, 249, 293, 354,
147, 214, 263, 312, 328,
155, 200, 237, 272, 297,
135, 188, 230, 280, 323,
159, 210, 252, 298, 331,
141, 189, 231, 275, 305,
159, 201, 248, 297, 338,
177, 236, 285, 350, 376,
134, 182, 220, 260, 296,
160, 208, 261, 313, 352,
143, 188, 220, 273, 314,
154, 200, 244, 289, 325,
171, 221, 270, 326, 358,
163, 216, 242, 281, 312,
160, 207, 248, 288, 324,
142, 187, 234, 280, 316,
156, 203, 243, 283, 317,
157, 212, 259, 307, 336,
152, 203, 246, 286, 321,
154, 205, 253, 298, 334,
139, 190, 225, 267, 302,
146, 191, 229, 272, 302,
157, 211, 250, 285, 323,
132, 185, 237, 286, 331,
160, 207, 257, 303, 345,
169, 216, 261, 295, 333,
157, 205, 248, 289, 316,
137, 180, 219, 258, 291,
153, 200, 244, 286, 324),
.Dim = c(30,5)))

(Note: the response data (Y) for the rats example can also be found in the file ratsy.odc in rectangular format. The covariate data (X) can be found in S-Plus format in file ratsx.odc. To load data from each of these files, focus the window containing the open data file before clicking on "Data" from the "Model" menu.)

Inits1 list(alpha = c(250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,
250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250),
beta = c(6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6),

alpha.c = 150, beta.c = 10,
tau.c = 1, sigma.alpha = 1, sigma.beta = 1)

Inits2 list(alpha = c(250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,
250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250),
beta = c(6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6),
alpha.c = 150, beta.c = 10,
tau.c = 1, tau.alpha = 1, tau.beta = 1)

Results

A 1000 update burn in followed by a further 10000 updates gave the parameter estimates:

node mean sd MC error 2.5% median 97.5% start sample
alpha0 106.6 3.65 0.04151 99.43 106.5 113.9 1001 10000
beta.c 6.185 0.1102 0.001294 5.967 6.185 6.404 1001 10000
sigma 6.074 0.4673 0.007724 5.247 6.044 7.068 1001 10000

Using Prior 2 (not recommended) has negligible impact

These results may be compared with Figure 5 of Gelfand et al 1990 --- we note that the mean gradient of independent fitted straight lines is 6.19.

Gelfand et al 1990 also consider the problem of missing data, and delete the last observation of cases 6-10, the last two from 11-20, the last 3 from 21-25 and the last 4 from 26-30. The appropriate data file is obtained by simply replacing data values by NA (see below). The model specification is unchanged, since the distinction between observed and unobserved quantities is made in the data file and not the model specification.

click on one of the arrows to open the data for the missing value analysis

Gelfand et al 1990 focus on the parameter estimates and the predictions for the final 4 observations on rat 26. These predictions are obtained automatically in BUGS by monitoring the relevant Y[] nodes. The following estimates were obtained:

We note that our estimate 6.58 of bc is substantially greater than that shown in Figure 6 of Gelfand et al 1990. However, plotting the growth curves indicates some curvature with steeper gradients at the beginning: the mean of the estimated gradients of the reduced data is 6.66, compared to 6.19 for the full data. Hence we are inclined to believe our analysis. The observed weights for rat 26 were 207, 257, 303 and 345, compared to our predictions of 204, 250, 295 and 341.