Calculating an index of colocation for a large dataset R

Question

I have a complicated, multi-part question. My apologies if I do not make myself clear. I am also a fairly novice R user, so forgive me if this seems rudimentary. I want to calculate a index of colocation for whale dive data and prey distribution data. This entails:

1. Calculating a frequency distribution of whale depth of dive data BY DIVE into depth bins from prey (fish and zoop) data.
2. For each dive, calculate the center of gravity (CG) and inertia (I).
3. For each dive, calculate a global index of colocation (GIC) vs. each prey type.

I want to be able to write a function (or series of functions) such that I do not have to separate my data by dive and rerun the functions for each dive manually.

Example whale data, where number if the dive number (sometimes 40+ dives), dive is equal to the depth and classification is related to the type of dive it is. [IMG]http://i41.tinypic.com/33vc5rs.jpg[/IMG]

Depth bins come from a separate data set containing prey information:

I have the following codes that work for the dive data as a whole, but need to write a loop or include an apply function such that I can run this for the data for each dive which is contained in a single file. So, for a whale with 40 dives, I need 40 whale frequencies, 40 whale CGs, 40 whale Is, etc. The prey distributionss are the SAME for each dive! Ultimately, I'd like a table which contains a list of the delta GIC values.

#bin whale dive depths
dive.cut=cut(whale,c(0 ,depths), right=FALSE) 
dive.freq=table(dive.cut) 

# compute CG 
fish.CG=sum(depths*fish)/sum(fish)
whale.CG=sum(depths*whale.freq)/sum(whale.freq)
zoop.CG=sum(depths*zoop)/sum(zoop)

# compute Inertia 
fish.I=sum((depths-fish.CG)^2*fish)/sum(fish)
whale.I=sum((depths-whale.CG)^2*whale.freq)/sum(whale.freq)
zoop.I=sum((depths-zoop.CG)^2*zoop)/sum(zoop)

#compute GIC as per 
# compute delta CG
deltaCG.fish_whale=fish.CG-whale.CG
GIC.fish_whale= 1-((deltaCG.fish_whale)^2/((deltaCG.fish_whale)^2+fish.I+whale.I))
deltaCG.zoop_whale=zoop.CG-whale.CG
GIC.zoop_whale= 1-((deltaCG.zoop_whale)^2/((deltaCG.zoop_whale)^2+zoop.I+whale.I))

UPDATES I have pasted example data for both prey and whale dives.

Prey data

 depths        fish       zoop
1      5     0.00000    0.000000
2     10     0.00000    0.000000
3     15     0.00000    0.000000
4     20    21.24194    0.000000
5     25   149.51694   14.937945
6     30   170.43214    0.000000
7     35   296.93453    0.737109
8     40    16.61643    4.295556
9     45    92.68130   26.384844
10    50    50.68548   55.902301
11    55    37.47343  218.673781
12    60    32.74443  204.452678
13    65    20.62983  113.112452
14    70    13.75121   83.014457
15    75    16.15562   55.051358
16    80    22.65562   96.746271
17    85    42.99768  302.229135
18    90 16315.65099  783.868978
19    95 43006.20482 1713.133161
20   100 23476.24740 3440.034642
21   105 30513.66346 6667.914707
22   110 17411.64500 9398.790964
23   115 12127.70195 7580.233165
24   120  4526.63393 7205.768739
25   125  3328.89644 6567.175766
26   130  1864.21486 4567.446886
27   135  2202.07464 4295.772442
28   140  2719.29417 4419.903403
29   145  1710.75599 5102.689940
30   150  2033.69552 4496.121974
31   155  2796.81788 3269.193606
32   160   967.09406 2310.203528
33   165   437.30896  447.940140
34   170   193.15526   63.731336
35   175   143.88043   38.004799
36   180   406.31373   22.565211
37   185   786.30087   31.889927
38   190  1643.52542   36.580063
39   195  1665.69794   14.084152
40   200  1281.15790    0.000000
41   205   753.75309   35.343794
42   210   252.48867    0.000000

Whale data:

  Number Dive Class
1       1 95.1     F
2       1 95.9     F
3       1 95.1     F
4       1 95.9     F
5       1 96.8     F
6       1 97.2     F
7       1 96.8     F
8       2 95.5     N
9       2 94.2     N
10      3 94.7     F
11      3 94.2     F
12      3 94.2     F
13      3 95.9     F
14      3 95.9     F
15      4 93.8     F
16      4 97.7     F
17      4 99.4     F
18      4 94.7     F
19      4 92.5     F
20      4 98.1     F
21      5 97.2     N
22      5 98.5     N
23      5 95.5     N
24      5 97.2     N
25      5 98.5     N
26      5 96.4     N
27      5 94.7     N
28      5 95.5     N

Answer 1

Give this code a try. I tested it out on the data you posted. I used the depths from the prey data frame. Not sure if that's what you wanted to do. And, this time I guessed that you used the whale$Dive for your dive.freq. If not, you'll have to change that. (Note, this question was cross-posted to the r-help list, too.)

prey <- structure(list(depths = c(5L, 10L, 15L, 20L, 25L, 30L, 35L, 40L, 
    45L, 50L, 55L, 60L, 65L, 70L, 75L, 80L, 85L, 90L, 95L, 100L, 
    105L, 110L, 115L, 120L, 125L, 130L, 135L, 140L, 145L, 150L, 155L, 
    160L, 165L, 170L, 175L, 180L, 185L, 190L, 195L, 200L, 205L, 210L
    ), fish = c(0, 0, 0, 21.24194, 149.51694, 170.43214, 296.93453, 
    16.61643, 92.6813, 50.68548, 37.47343, 32.74443, 20.62983, 13.75121, 
    16.15562, 22.65562, 42.99768, 16315.65099, 43006.20482, 23476.2474, 
    30513.66346, 17411.645, 12127.70195, 4526.63393, 3328.89644, 
    1864.21486, 2202.07464, 2719.29417, 1710.75599, 2033.69552, 2796.81788, 
    967.09406, 437.30896, 193.15526, 143.88043, 406.31373, 786.30087, 
    1643.52542, 1665.69794, 1281.1579, 753.75309, 252.48867), zoop = c(0, 
    0, 0, 0, 14.937945, 0, 0.737109, 4.295556, 26.384844, 55.902301, 
    218.673781, 204.452678, 113.112452, 83.014457, 55.051358, 96.746271, 
    302.229135, 783.868978, 1713.133161, 3440.034642, 6667.914707, 
    9398.790964, 7580.233165, 7205.768739, 6567.175766, 4567.446886, 
    4295.772442, 4419.903403, 5102.68994, 4496.121974, 3269.193606, 
    2310.203528, 447.94014, 63.731336, 38.004799, 22.565211, 31.889927, 
    36.580063, 14.084152, 0, 35.343794, 0)), .Names = c("depths", 
    "fish", "zoop"), class = "data.frame", row.names = c("1", "2", 
    "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", 
    "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", 
    "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", 
    "37", "38", "39", "40", "41", "42"))

whale <- structure(list(Number = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
    3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 
    5L, 5L, 5L), Dive = c(95.1, 95.9, 95.1, 95.9, 96.8, 97.2, 96.8, 
    95.5, 94.2, 94.7, 94.2, 94.2, 95.9, 95.9, 93.8, 97.7, 99.4, 94.7, 
    92.5, 98.1, 97.2, 98.5, 95.5, 97.2, 98.5, 96.4, 94.7, 95.5), 
    Class = c("F", "F", "F", "F", "F", "F", "F", "N", "N", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "N", "N", 
    "N", "N", "N", "N", "N", "N")), .Names = c("Number", "Dive", 
    "Class"), class = "data.frame", row.names = c("1", "2", "3", 
    "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
    "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
    "27", "28"))

# split the data frame into a list with a different element for each dive
dives <- split(whale, whale$Dive)

# define a single function that does all of your computations
compute <- function(whale, depths, fish, zoop) {
    # you don't say what part of the whale data you are counting ... I'll assume it's the dive
    dive.freq <- table(cut(whale$Dive, c(0, depths)))
    #compute Center of Gravity
    fish.CG <- sum(depths*fish)/sum(fish) #calculate CG for fish distribution ONCE for each whale
    zoop.CG <- sum(depths*zoop)/sum(zoop) #calculate CG for zoop distribution ONCE for each whale
    whale.CG <- sum(depths*dive.freq/sum(dive.freq)) #calculate for EACH dive
    #compute Inertia
    fish.I <- sum((depths-fish.CG)^2*fish)/sum(fish) 
    zoop.I <- sum((depths-zoop.CG)^2*zoop)/sum(zoop)
    whale.I <- sum((depths-whale.CG)^2*dive.freq)/sum(dive.freq) #needs to be calculated for EACH dive
    # compute delta CG
    deltaCG.fish_whale <- fish.CG-whale.CG
    GIC.fish_whale <- 1-((deltaCG.fish_whale)^2/((deltaCG.fish_whale)^2+fish.I+whale.I))
    deltaCG.zoop_whale <- zoop.CG-whale.CG
    GIC.zoop_whale <- 1-((deltaCG.zoop_whale)^2/((deltaCG.zoop_whale)^2+zoop.I+whale.I))
    # then list off all the variables you want to keep as output from the function here
    c(fish.CG=fish.CG, whale.CG=whale.CG, zoop.CG=zoop.CG, fish.I=fish.I, whale.I=whale.I, zoop.I=zoop.I, 
        GIC.fish_whale=GIC.fish_whale, GIC.zoop_whale=GIC.zoop_whale)
    }

# apply the compute function to each element of the dives list
t(sapply(dives, function(dat) compute(whale=dat, depths=prey$depths, fish=prey$fish, zoop=prey$zoop)))