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
