2011-09-21 101 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