Dear
all,
I'm new to Bayesian stats and WinBUGS, and I'd really appreciate
your advice
regarding a hierarchical Bayes model, and bayesian jargon/usage
in general.
I'm trying to build a hierarchical model for a meta-analysis
study.
Basically, the model has two components: a lower-level model that I
think
would be a data model, and a higher-level model that I think is a
process
model (below I explain each in depth). I want to propagate the
uncertainty
from the original samples through to the posterior of the
process model. I'm
particularly intersted in the effect of Type (see below),
as well as in
variance conponents. I wrote WinBUGS codes both for the data
and process
models, and both codes run nicely on their own (with simulated
data).
However, I've been unable to combine both models into a single
hierarchical
Bayes model; here's where a hand would be very much
appreciated.
So, the idea is the following: I want to model the response
variable H as a
function of some fixed and random effects using a nested,
mixed model. This
would be the process model. However, H is a sort of random
difference,
inferred from data available in the literature. Specifically, H
is
calculated using the following data model:
H = abs(mean1 - mean2)
/ SDp
where abs is absolute value, mean1 and mean2 are mean trait values
of
samples 1 and 2 with sample size N1 and N2, and SDp is the pooled
standard
deviation of both samples. Means, SDs, and Ns are extracted from
the
literature, and the posterior distribution of H is calculated
propagating
the uncertainty from means and SDs (see code below), which are a
function of
sample size. This is done for each i entry in the
database.
Now, I need to model the (random) response variable H as a
function of some
fixed and random factors, for example using the following
process model:
H(i) = Intercept + Type + System(type) + e
where
Intercept is the overall mean, Type is a fixed factor with two levels,
System(Type) is a random factor nested in Type, and e is the residual
error.
Both, the data model to obtain the posterior of H for each entry,
and the
process model to obtain the posterior for the effects of Type and
System,
are working nicely when run independently from each other. But, how
can I
put them into a single hierarchical model????
I was checking
Clark and Gelfand's book, and one option seems to be using
'empirical
Bayes'. I think this would involve the following steps: 1) run
the data
model and record H-hat as well as the SD of H. 2) Assume that the
posterior
of H is normally distributed. And 3), fit the process model to the
H-hat
data and use SD of H as an additional source of error at the right
side of
the equation. However, the assumption in step 2) can be problematic,
because
we're interested in the absolute value of differences, and
distributions can
therefore be biased to the right. A proper hierarchical
Bayes analysis
should be much more robust.
Please, help me see the light...
Below
are the two (still) independent codes with some simulated data to run
them
independently.
Any clue/critique will be
appreciated.
Cheers,
Cristian
-----------------------------------------------------------------------------
Cristián
Correa Guzmán, PhD Candidate
Redpath Museum, Dept of Biology, McGill
University
859 Sherbrooke Street West
Montreal, Quebec H3A
2K6
http://biology.mcgill.ca/grad/cristian/
----------------------------------------------------------------------------
Data
model:
Model{
for (i in 1:5) { # for
five cases compiled from the literature (in
#reality there're hundreds)
M1[i] ~
dnorm(X1[i], tau1[i]) # expected mean of sample 1
M2[i] ~ dnorm(X2[i],
tau2[i])
SD1[i] ~
dnorm(S1[i], tauSD1[i]) I(0, ) # expected SD of sample 1; >= 0.
SD2[i] ~
dnorm(S2[i], tauSD2[i]) I(0, )
SE1[i] <-
SD1[i]/sqrt(N1[i]) #SE
of observed mean 1 (SE of X1)
SE2[i] <-
SD2[i]/sqrt(N2[i])
tau1[i] <- 1
/ (SE1[i] * SE1[i]) #Presicion for
expected mean 1.
tau2[i] <- 1/ (SE2[i] * SE2[i])
tauSD1[i] <- 1
/ pow((S1[i] * sqrt(2 / N1[i])),2) #Presicion for expected SD1
tauSD2[i]
<- 1 / pow((S2[i] * sqrt(2 / N2[i])), 2) # SD of SD formula
var1[i]
<- pow(SD1[i], 2) # Expected variance 1
var2[i] <- pow(SD2[i],
2)
SS1[i] <- var1[i] * (N1[i]-1) # Expected sum of squares
SS2[i]
<- var2[i] * (N2[i]-1)
SDP[i] <- sqrt((SS1[i] + SS2[i]) / (N1[i] +
N2[i] - 2)) # Expected SD pooled
H[i] <- abs(M1[i] - M2[i])/SDP[i] #
Expected H (posterior)
}}
Data
#here a few
data entries (rows) to try out the code,
#the real
dataset has hundreds of rows nested in systems and types (see next
model)
N1[] N2[] X1[] X2[] S1[] S2[]
5 5 25 27 3 3
10 10 25 27 3
3
20 20 25 27 3 3
50 50 25 27 3 3
100 100 25 27 3
3
END
------------------------------------------------------------------------------------------
Process
model:
#MODEL y = a + cont + Sys(cont) + e,
#With X fixed and
Sys(cont) random and nested in cont.
# Nesting was achieved through a double
loop.
model{
a ~ dnorm(0, 1.0E-6) # prior for intercept (global
mean)
type[1] ~ dunif(0,50) # prior for type (symetric offset from
a)
type[2] <- type[1]*-1
SDSys ~ dunif(0, 100) #prior, among System
variation
tauSys <- 1/(SDSys*SDSys)
RMSE ~ dunif(0, 100) # prior,
within System variation. Root mean square error
tau <-
1/(RMSE*RMSE)
for (c in 1:2){
for (s in 1:6){
# For each System nested in
Type
# Nesting is achieved by double
loop.
mSys[c, s] ~ dnorm(0, tauSys) # randdom effect of system nested in
Type
}}
for (i in 1:144){ #For each record (H,
from the data model)
H[i] ~ dnorm(mH[i], tau)
mH[i] <- a +
type[X[i]] + mSys[X[i],System[i]]
}}
Initials list(a = 10, type =
c(5, NA), SDSys = 3) # initil values
Data
# Here some simulated data
where H is provided as an observed value.
# The idea is
to provide means, SDs, and Ns, however,
# and let H be
a random variable based on the data model.
X[] System[] H[]
1 1
21.40526982
1 1 22.96069847
1 1 22.08790818
1 1 21.42811912
1 1
20.4203889
1 1 20.87540562
1 1 22.69079362
1 1 20.63009231
1 1
21.50064958
1 1 20.94080642
1 1 19.53246408
1 1 19.67526267
1 2
23.70950204
1 2 22.72182131
1 2 22.86289987
1 2 22.79907877
1 2
22.58214768
1 2 21.70751405
1 2 23.46384169
1 2 21.95452943
1 2
21.4559775
1 2 21.90409722
1 2 22.58854097
1 2 22.33987684
1 3
20.02665508
1 3 22.01765585
1 3 21.87157499
1 3 21.28572564
1 3
19.52806341
1 3 20.30449887
1 3 19.90207673
1 3 22.19582165
1 3
20.25635854
1 3 19.9580788
1 3 20.49807313
1 3 21.77696825
1 4
21.54902098
1 4 22.16761695
1 4 23.37282613
1 4 22.08414268
1 4
20.39198022
1 4 20.00225098
1 4 21.62097735
1 4 21.39075705
1 4
21.62931963
1 4 20.94289489
1 4 20.89058214
1 4 21.17044182
1 5
23.93100653
1 5 21.91187502
1 5 22.04685022
1 5 20.43306248
1 5
21.6970992
1 5 22.36825873
1 5 21.9484764
1 5 21.24235915
1 5
23.83174111
1 5 20.52429564
1 5 24.12862272
1 5 22.12816677
1 6
22.3585235
1 6 23.68226229
1 6 23.91913943
1 6 23.41560378
1 6
23.38594401
1 6 23.18804779
1 6 22.45934697
1 6 23.43975144
1 6
23.20264347
1 6 22.78201585
1 6 22.63947057
1 6 23.15890862
2 1
18.90443667
2 1 16.52527693
2 1 16.39718828
2 1 18.83546672
2 1
16.95609958
2 1 16.55602424
2 1 15.93888648
2 1 16.24977523
2 1
16.74221991
2 1 16.65873049
2 1 17.05805293
2 1 16.44065973
2 2
16.98697446
2 2 15.60740868
2 2 16.0380606
2 2 16.32339374
2 2
15.85330468
2 2 16.0702883
2 2 16.23681495
2 2 16.08152729
2 2
17.52479701
2 2 16.68068865
2 2 15.78972641
2 2 17.1039102
2 3
18.83569314
2 3 16.97038274
2 3 18.78421037
2 3 16.57855431
2 3
15.47891875
2 3 16.67819926
2 3 16.42623461
2 3 19.02827806
2 3
17.14768946
2 3 15.4561619
2 3 16.04244614
2 3 17.15557746
2 4
18.46772785
2 4 17.15541072
2 4 17.83362105
2 4 17.75359394
2 4
18.51484984
2 4 19.39443054
2 4 19.47705076
2 4 18.09696832
2 4
20.75061711
2 4 18.9564764
2 4 19.84962105
2 4 19.4744911
2 5
17.22561034
2 5 20.16334588
2 5 18.39367953
2 5 18.42276135
2 5
19.17101536
2 5 18.5779392
2 5 16.60984269
2 5 17.52040153
2 5
16.8098906
2 5 20.01838996
2 5 18.04699373
2 5 17.82742983
2 6
21.12148807
2 6 20.74469858
2 6 19.84250989
2 6 19.23402202
2 6
20.44772861
2 6 21.12731162
2 6 20.37193643
2 6 20.70035173
2 6
18.38081773
2 6 18.62533173
2 6 18.67498957
2 6
18.96553952
END
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list