2011-09-21 204 views
0

建議使用adehabitat來計算交集量後,我偶然發現了一個輕微的(希望是簡單的)問題。在這個庫中,我使用了kerneloverlap命令,因爲我需要計算交集的體積。我想知道你是否可以幫助我解決一些編程問題。我需要修改腳本以使其「批量」處理友好。我知道足夠的R讓自己陷入麻煩並失去理智,因爲我知道某些事情應該是可能的,但無法弄清楚如何讓它工作。修改代碼以批量處理r

的命令是非常簡單:

kerneloverlap(loc[,c("X","Y")], loc$year, lev = 90, grid=30, meth="VI", conditional=TRUE) 

它從數據文件LOC取x,y座標,由年,並且計算相交的體積爲30的網格單元尺寸在利用分配90.

輸入文件(請參閱下面摘錄)是援助,X,Y,年和季節。對於這個例子,只有1個賽季(記住我有3個賽季)。對於這個例子,我想比較每個交叉點之間每年的一個季節。所以測試數據有2年1季和2個人。我希望能夠說的是「2003年至2004年間產犢季節的動物1交叉口的體積爲0.8,這表明交疊和保真度高的地方」。

我還想再比較一下季節。使得它在夏季路口動物1和越冬季節,2003年的成交量爲0.04指示重疊的低水平,並沒有忠誠的位置」

有些事情要記住:並非所有的人都存在每年或每個季節都活着,因此可能需要某種類型的水滴

這是我迄今爲止的R腳本(它不起作用)請注意,輸出沒有很好地連接在一起,似乎無法得到一個編譯的文件,它喜歡它告訴我什麼年份,​​個人或季節,它是比較東西與

IDNames= levels(loc$anid) 
Year = unique(loc$year) 
for (i in 1:(length(IDNames))){ 
vi90 = kerneloverlap(loc[,c("X","Y")], loc$year, lev = 90, grid=30, meth="VI", conditional=TRUE) 
    } 
colnames(vi)= c(paste(IDNames[i],Year[n], sep =""),paste(IDNames[i], Year[n], sep ="")) 
} 
write.csv(vi,"VolInter_indiv.csv") 


    structure(list(anid = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L 
), .Label = c("c_002", "c_104"), class = "factor"), X = c(276646.0514, 
276485.0397, 278102.4193, 278045.4716, 278993.8807, 274834.5677, 
278516.0218, 296741.8328, 299080.2451, 291874.5068, 168540.0024, 
168360.8211, 169538.2299, 164538.2592, 157321.7524, 148090.3478, 
140575.2442, 133369.7162, 134375.0805, 138763.5342, 232347.5137, 
231989.4609, 231793.1066, 234923.4012, 233374.4531, 232256.4667, 
233660.3445, 239317.3128, 246354.664, 145161.8922, 144148.7895, 
145154.7652, 145399.3515, 144581.4836, 143646.7295, 145055.3165, 
144613.1393, 145037.3035, 144701.2676), Y = c(2217588.648, 2216616.387, 
2219879.777, 2220818.804, 2216908.127, 2220423.322, 2216589.91, 
2234167.287, 2239351.696, 2232338.072, 2273737.333, 2273954.782, 
2269418.423, 2271308.607, 2264694.484, 2263710.512, 2254030.274, 
2253352.426, 2248644.946, 2262359.026, 2231404.821, 2229583.89, 
2231700.485, 2231598.882, 2237122.967, 2233302.185, 2240092.997, 
2237702.817, 2249213.958, 2261841.308, 2263064.156, 2262236.452, 
2264147.03, 2263214.877, 2263336.363, 2261417.946, 2256289.995, 
2256694.953, 2253352.576), year = c(2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2004L, 2004L, 
2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 
2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 2004L, 2004L), season = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L), .Label = "calving", class = "factor")), .Names = c("anid", 
"X", "Y", "year", "season"), class = "data.frame", row.names = c(NA, 
-39L)) 

回答

6

好吧,我會咬。

您的代碼有一些錯別字(我希望),使它無法運行。讓我們把它扔出去並重新開始。函數kerneloverlap爲其第二個參數中指定的每對項目返回一個重疊值矩陣。在你的第一個例子中,你要比較幾年。

首先讓我們來想象我們會用數據做只是一種動物,並編寫輸出我們想要的值只是簡單的情況下的函數:

kernMod <- function(x){ 
    #x is the data for a single animal 
    rs <- kerneloverlap(x[,c("X","Y")], 
         x$year,lev = 90, 
         grid = 30, 
         meth = "VI", 
         conditional = TRUE) 
    #Assumes we're only comparing two years 
    out <- data.frame(year = paste(colnames(rs),collapse="-"), val = rs[2,1]) 
    out 
} 

現在我們可以應用此分別每隻動物:

kernMod(subset(loc,anid == 'c_002')) 
     year val 
1 2003-2004 0 
> kernMod(subset(loc,anid == 'c_104')) 
     year  val 
1 2003-2004 0.06033966 

,或者我們可以使用ddplyplyr包將其應用到每個動物依次爲:

ddply(loc,.(anid),.fun = kernMod) 

    anid  year  val 
1 c_002 2003-2004 0.00000000 
2 c_104 2003-2004 0.06033966 

要包括多發季節,您只需將它添加到變量列表中ddply分裂了(未經測試):

ddply(loc,.(anid,season),.fun = kernMod) 

要一年內四季之間的比較,您將需要修改kernMod傳遞x$season作爲第二個參數,然後調用類似(未經測試):

ddply(loc,.(anid,year),.fun = kernMod) 

如果完全數據中有多年,kernMod將需要一些更多的修改,因爲kerneloverlap返回正x n矩陣,其中n是數據中的年數。也許是這樣的(未經測試)

kernMod <- function(x){ 
    #x is the data for a single animal 
    rs <- kerneloverlap(x[,c("X","Y")], 
         x$year,lev = 90, 
         grid = 30, 
         meth = "VI", 
         conditional = TRUE) 
    rs[lower.tri(rs,diag = TRUE)] <- NA 
    rs <- melt(rs) 
    rs <- subset(rs, !is.na(value)) 
    out <- data.frame(year = paste(rs$X1,rs$X2,collapse="-"), val = rs$value) 
    out 
} 

這種做法應該處理僅計算您的數據值「失蹤」的動物。

好的。我很樂意爲此獲得第3-4名作者,但我願意承認。 ;)

+0

謝謝!太棒了。我很感激你花時間來解釋這些步驟。它幫助我邏輯推理其他腳本,因爲我發現自己需要爲各種計算執行相同的步驟。再次感謝你!雖然你可能一直在開玩笑,但我完全承認你! – Kerry

+0

我整個上午都在測試這個腳本,並將這個批量處理的結果與一次只有一個人的數據進行比較,結果我得到了完全不同的結果。 – Kerry