尝试在 R 中创建并循环不平衡数据矩阵

Question

我正在尝试进行分层贝叶斯分析，但 R 和 WinBUGS 代码遇到了一些问题。我没有平衡的数据，并且正在努力编码。我每天使用 iButtons（温度记录设备）收集横断面的温度数据，并尝试生成一个将其与遥感数据相关的模型。不幸的是，每个横断面都有不同数量的 iButton，因此在横断面（j）中创建按钮（i）的 3D 矩阵，在日（t）重复“采样”对我来说是一个问题。

最终，我的模型将类似于：

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

2 级 b0 = a00 + a01*y1 + 。 。 。安*因 b1 = a10 + a11*y1 ...

3级（也许？） - 随机2级拦截

通常我会这样做： 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]
}

然后我可以有一个 3D 矩阵，我可以在 WinBUGS 或 OpenBUGS 中循环遍历这样：

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]
}}}

无论如何，不​​要担心上面代码的细节，它只是作为其他分析的示例放在一起的。我的主要问题是，当我没有每个横断面具有相同数量 iButton 的平衡设计时，如何进行此类分析？任何帮助将不胜感激。我显然对 R 和 WinBUGS 很陌生，并且以前没有太多的计算机编码经验。

谢谢！

哦，这是长（堆叠）格式的数据：

    > 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

Answer 1

创建长度向量或数组并使用子索引。
用你的例子：

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]
}}}

我认为它会起作用，但我没有测试过，因为没有数据。

Answer 2

您可以尝试使用 list 来代替吗？

这允许列表中每个项目的长度可变，其中每个索引对应于横断面。

所以像这样：

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])
      }
    }
  }