0

I am trying to conduct an hierarchical bayesian analysis but am having a little trouble with R and WinBUGS code. I don't have balanced data and am struggling with the coding. I have temperature data collected daily with iButtons (temperature recording devices) in transects and am trying to generate a model that relates this to remote sensing data. Unfortunately, each transect has a different number of iButtons so creating a 3D matrix of button(i), in transect(j), repeatedly "sampled" on day(t) is a problem for me.

Ultimately, my model will be something like:

Level 1 Temp[ijk] ~ N(theta[ijk], tau) theta[ijk] = b0 + b1*x1 + . . . + bn*xn

Level 2 b0 = a00 + a01*y1 + . . . an*yn b1 = a10 + a11*y1 ...

Level 3 (maybe?) - random level 2 intercepts

Normally I would do something like this: Wide <- reshape(Data1, idvar = c("iButton","block"), timevar = "julian", direction = "wide")

J <- length(unique(Data$block))
I <- length(unique(Data$iButton))
Ti <- length(unique(Data$julian))

Temp <- array(NA, dim = c(I, Ti, J))

for(t in 1:Ti) {
sel.rows <- Wide$block == t
Temp[,,t] <- as.matrix(Wide)[sel.rows, 3:Ti]
}

Then I could have a 3D matrix that I could loop through in WinBUGS or OpenBUGS as such:

for(i in 1:J) {          # Loop over transects/blocks
  for(j in 1:I) {        # Loop over buttons
    for(t in 1:Ti) {     # Loop over days
    Temp[i,j,t] ~ dnorm(theta[i,j,t])    
    theta[i,j,t] <- alpha.lam[i] + blam1*radiation[i,j] + blam2*cwd[i,j] + blam3*swd[i,j]
}}}

Anyway, don't worry about the details of the code above, it's just thrown together as an example from other analyses. My main question is how to do this type of analysis when I don't have a balanced design with equal numbers of iButtons per transect? Any help would be greatly appreciated. I'm clearly new to R and WinBUGS and don't have much previous computer coding experience.

Thanks!

oh and here is what the data look like in long (stacked) format:

    > Data[1:15, 1:4]
   iButton julian block       aveT
1        1      1     1 -4.5000000
2        1      2     1 -5.7500000
3        1      3     1 -3.5833333
4        1      4     1 -4.6666667
5        1      5     1 -2.5833333
6        1      6     1 -3.0833333
7        1      7     1 -1.5833333
8        1      8     1 -8.3333333
9        1      9     1 -5.0000000
10       1     10     1 -2.4166667
11       1     11     1 -1.7500000
12       1     12     1 -3.2500000
13       1     13     1 -3.4166667
14       1     14     1 -2.0833333
15       1     15     1 -1.7500000
djhocking
  • 1,072
  • 3
  • 16
  • 28
  • no time to write a whole answer now but... If you've got a matrix theta and you want dnorm at each point try...Temp <- dnorm(theta). And then, when you look at Temp and the miracle that occurred, consider thinking completely differently about your problem. – John Dec 17 '11 at 00:59
  • I don't have a matrix theta. I would use WinBUGS as an MCMC machine to estimate theta and the beta terms. Basically my problem is I have Temp data at i, j, t but the length of i differs for each j so I can't figure out how to code Temp to be in the right form to loop through it. The only other thing I can think of that my current, poor coding skills would allow would be to have an [i,t] matrix and then use block (j) as a covariate using dummy variables. If the design was balanced I would code it as above but that won't work in this case. – djhocking Dec 17 '11 at 03:20

2 Answers2

4

Create a vector or array of lengths and use subindexing. Using your example:

J <- length(unique(Data$block))
I <- tapply(Data$iButton, Data$block, function(x) length(unique(x))
Ti <- tapply(Data$julian, list(Data$iButton, Data$block), function(x) length(unique(x))


for(i in 1:J) {          # Loop over transects/blocks
  for(j in 1:I[i]) {        # Loop over buttons
    for(t in 1:Ti[i, j]) {     # Loop over days
    Temp[i,j,t] ~ dnorm(theta[i,j,t])    
    theta[i,j,t] <- alpha.lam[i] + blam1*radiation[i,j] + blam2*cwd[i,j] + blam3*swd[i,j]
}}}

I think it would work, but I haven't tested since there no data.

Luciano Selzer
  • 9,806
  • 3
  • 42
  • 40
  • @Iselzer - Thanks, I think this will work. I have never used subindexing before (like I said, I'm new to this). Tomorrow I will have to think more about how to get my temperature data (Temp) into the right format. One challenge is that WinBUGS cannot handle any NAs in the data so my matrix of Temp[i,j,t] will have to have a different size 2D matrix for each level of block J. – djhocking Dec 17 '11 at 07:44
  • @Iselzer - I will when I ca but it won't let me yet because my reputation is only 11 and it has to be 15 to upvote. – djhocking Dec 17 '11 at 15:26
  • @djhocking: But if this is an effective answer, then surely you could "checkmark" it since you asked the question. – IRTFM Dec 18 '11 at 14:58
1

Can you try using a list instead?

This allows a variable length for each item in the list where each index would correspond to the transect.

So something like this:

theta <- list()

for(i in unique(Data$block)) {
  ibuttons <- unique(Data$iButton[Data$block==i])
  days <- unique(Data$julian[Data$block==i])
  theta[[i]] <- matrix(NA, length(ibuttons), length(days)) # Empty matrix with NA's
    for(j in 1:length(ibuttons)) {
      for(t in 1:length(days)) {
        theta[[i]][j,t] <- fn(i, ibuttons[j], days[t])
      }
    }
  }
oeo4b
  • 1,574
  • 1
  • 10
  • 5
  • Thanks, I am an inexperienced programmer, I had never thought to create a list of matrices to cycle through. This would definitely work if I did my MCMC in R. I don't know if the BUGS language handles indexing this way. I use R2WinBUGS to pass my data to WinBUGS to take advantage of it's simple, transparent language and use it's Gibbs sampler. If WinBUGS can't handle this for some reason, I may just try doing my MCMC sampling in R. There are plenty of packages and writing my own sampler wouldn't be too bad for a simple model. – djhocking Dec 17 '11 at 19:33
  • I only have a limited amount of experience with WinBUGS, but dealing with `list` objects in R definitely comes in handy when you're dealing with multiple matrices of varying dimensions. – oeo4b Dec 18 '11 at 07:10