2011-12-17 74 views
0

我試圖進行分層貝葉斯分析,但在R和WinBUGS代碼中遇到了一些問題。我沒有平衡的數據,並且在編碼上掙扎。我每天都用iButtons(溫度記錄設備)在橫斷面上收集溫度數據,並試圖生成一個與遙感數據相關的模型。不幸的是,每個橫斷面的iButton數量都不相同,因此在橫斷面(j)中創建按鈕(i)的3D矩陣在第(t)天重複「採樣」對我來說是一個問題。試圖創建並通過R中的非平衡數據矩陣循環

最終,我的模式將是這樣的:

級別1 溫度[IJK]〜N(THETA [IJK],tau蛋白) THETA [IJK] = B0 + B1 * X1 +。 。 。 + bn * xn

2級 b0 = a00 + a01 * y1 +。 。 。 (?也許)有* YN B1 = A10 + A11 * Y1 ...

等級3 - 隨機2級攔截

通常我會做這樣的事情: 寬< - 重塑(數據1,idvar = C( 「iButton的」, 「塊」),timevar = 「儒略」,方向= 「寬」)

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

無論如何,不​​要擔心上面的代碼的細節,它只是作爲其他分析的一個例子。我的主要問題是如何進行這種類型的分析,當我沒有平均的設計和相同數量的iButtons橫斷面?任何幫助將不勝感激。我明顯是R和WinBUGS的新手,並沒有很多以前的計算機編碼經驗。

謝謝!

OH和這裏是數據看起來像在長(堆疊)格式的內容:

> 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 
+0

沒時間寫信全部答案,但是...如果你已經有了一個矩陣theta和你想在每一點嘗試dnorm ...溫度< - dnorm (THETA)。然後,當你看到Temp和發生的奇蹟時,考慮完全不同的問題。 – John 2011-12-17 00:59:55

+0

我沒有矩陣theta。我將使用WinBUGS作爲MCMC機器來估計theta和beta項。基本上我的問題是,我有臨時數據在i,j,t但是我的長度不同,每個j,所以我不知道如何編寫Temp以正確的形式來循環它。我能想到的另外一件事是,我目前的編程能力不足將會產生[i,t]矩陣,然後使用塊(j)作爲使用虛擬變量的協變量。如果設計是平衡的,我會像上面那樣編碼,但在這種情況下不起作用。 – djhocking 2011-12-17 03:20:13

回答

1

你可以嘗試使用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]) 
     } 
    } 
    } 
4

創建長度的矢量或陣列和使用subindexing。使用 你的例子:

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

我認爲這是可行的,但我因爲沒有數據還沒有測試。

+0

@Iselzer - 謝謝,我認爲這會奏效。我從來沒有使用過subindexing(就像我說的,我是新來的)。明天我將不得不更多地考慮如何將我的溫度數據(Temp)轉換爲正確的格式。一個挑戰是WinBUGS無法處理數據中的任何NA,所以我的Temp [i,j,t]矩陣必須爲塊J的每個級別具有不同大小的2D矩陣。 – djhocking 2011-12-17 07:44:07