0

階層ベイジアン分析を実行しようとしていますが、R と WinBUGS コードに少し問題があります。バランスの取れたデータがなく、コーディングに苦労しています。トランセクトで iButton (温度記録デバイス) を使用して毎日収集された温度データがあり、これをリモート センシング データに関連付けるモデルを生成しようとしています。残念ながら、各トランセクトには異なる数の iButton があるため、トランセクト (j) でボタン (i) の 3D マトリックスを作成し、day(t) に繰り返し「サンプリング」することは私にとって問題です。

最終的に、私のモデルは次のようになります。

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

レベル 2 b0 = a00 + a01*y1 + . . . an*yn 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]
}

次に、WinBUGS または OpenBUGS でループできる 3D マトリックスを作成できます。

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
4

2 に答える 2

4

長さのベクトルまたは配列を作成し、サブインデックスを使用します。あなたの例を使用して：

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

うまくいくと思いますが、データがないのでテストしていません。

于 2011-12-17T04:23:54.967 に答える
1

list代わりにa を使用してみてください。

これにより、各インデックスがトランセクトに対応するリスト内の各項目の可変長が可能になります。

だから、このようなもの：

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])
      }
    }
  }
于 2011-12-17T04:25:48.207 に答える