Revision as of 07:49, 31 July 2008

This page describes utilization of Hierarcical Bayesin method. The following page briefly introduces theory behind the method, and illustrates the use the method by examples.

Introduction

The Bayesian method

The Bayes method is based on the bayesian rule. The following equations can be found from the reference 1, pages 7-8. The Bayesian rule is:

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

,where p(θ) is prior distribution and on p(y|θ) is distribution acquired from the data. By contidionalizing θ with y (which is known), we get:

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

By conditionalizing this way, we obtain p(θ|y), which is prior probability conditionalised by the data. This is also called posterior probability. With these equations, we can adjust presuppositions about probability density function of the data using data itself. As a end result we obtain upgraded probability density function, which is called posterior density.   

Posterior distribution ~ prior distribution x distribution estimated by the data

Hierarcical Bayesian method

We create an example of the method in this article. We study dioxin concentrations in different sea area which we assume to be share a similar form of distribution. In other words, we assume that the concentraions in different areas are samples of the total concentration distribution. Thus, we can use observations of concentrations in all areas to estimate concentrations in one particular area of interest. This means conditionalizing observated values to certain parameters. We assume probabilites to these parameters with new parameters, which are called hyperparameters.

There is an extensive explanation of the idea and the method in reference 1, pages 117-121, and you can find a picture of the method in page 119. In our example, we try to create a model which gives us a relatively accurate concentration estimate for all different areas.

Material and methods

Methods

The method used herein is called hierarcical Bayes method. Theory can be read from reference 1. We used WinBUGS and R as a modelling tool. 

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

You can get the program and updates from the address above. They also include quite extensive manuals. You need to copy the text to WinBUGS to a new file, and select decode from tools menu and follow the manual. 

Table 1. Names of the variables and parameters in the model. 

Picture 1. Picture of the model in Win BUGS

Material

The data we used in our example consists of dioxin concentration of herring with ages, sizes and location

We create a model similar to example model in WinBUGS. Variables are analogous as demonstrated below:
Variables:

Weight of rat = Concentration of herring
Age of rat = Age of herring
Number of rat (specimen #) = Location

We ran the model with Win BUGS with herring data, and examined the distribution of variables of interest.In rat example the model estimates distributions for weight for each rat. In our model we estimate concentration distributions for each location.

Results

Descriptive statistics

Kuva 2. Age against concentration in herring

Image:Ikä_vs_poitoisuus_silakka.jpg

Kuva 3.Weight against concentration in herring

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ä.

Results of WinBUGS calculations

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 3. WinBUGSin tuloslistaus silakkadatalle kolmen vuoden ikä-intervallilla.

Syntaxes and iterations 

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

Conculsions

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

Rat example

This model can be found from WinBUGS help page. Our example is based directly to this model.

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 "noninformativepriors, 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.