ABindoff / TwilightFree

R package for TwilightFree method of geolocation from (noisy) light data
8 stars 1 forks source link

Error in if (s == 0) x else x/sum(x) :missing value where TRUE/FALSE needed #3

Open VirginiaMorera opened 6 years ago

VirginiaMorera commented 6 years ago

Hi,

I am testing the method with the data obtained from a Migrate Tech geolocator. The data looks like this:

head(d.lig) 
# A tibble: 6 x 3
   Date                Light Temp
    <dttm>              <dbl> <lgl>
1 2015-07-05 17:20:11 655   NA
2 2015-07-05 17:25:11 655   NA    
3 2015-07-05 17:30:11  85.2 NA    
4 2015-07-05 17:35:11  89.7 NA    
5 2015-07-05 17:40:11  80.6 NA    
6 2015-07-05 17:45:11  80.6 NA

The results from the calibrate function look like this: imagen "max light in night window: 537.118 assuming a solar zenith angle of: 97"

when I use a single day, but give a warning when I use several days:

day <- seq(as.POSIXct("2015-07-06 00:00:00", "GMT"), as.POSIXct("2015-07-10 00:00:00", "GMT"), "day") imagen [1] "max light in night window: 537.118 assuming a solar zenith angle of: 97" Warning messages: 1: In>=.default(df$Date, as.POSIXct(day, tz = "GMT")) : longer object length is not a multiple of shorter object length 2: In<.default(df$Date, as.POSIXct(day + 24 * 60 * 60, tz = "GMT")) : longer object length is not a multiple of shorter object length

When I run the model with just 1 day of calibration (to avoid the warning) with these parameters:

model <- TwilightFree(d.lig,
                      alpha=c(1, 1/10),
                      beta=c(1, 1/4),
                      zenith = zen, threshold = thresh,
                      deployed.at = c(-15.788, 27.844), 
                      retrieved.at = c(-15.788, 27.844))

using the function

fit <- SGAT::essie(model,grid,epsilon1=1.0E-4, epsilon2 = 1.0E-4)

where the grid was made with

grid <- makeGrid(c(-62, 80), c(-47, 70), cell.size = 1, pacific = F)

Apparently the model starts to run, as I get the first two Fwd 366 Bwd 365

but it stops and gives the error message

Error in if (s == 0) x else x/sum(x) :missing value where TRUE/FALSE needed

Any idea what I might be doing wrong?

Thanks

ABindoff commented 6 years ago

Thanks for opening an issue (and for the detail), it is a great way to disseminate information to the community of TwilightFree users.

The warning generated by calibrate is nothing to worry about (you've told it to look at a sequence of dates, it's expecting a single date, it will just calculate a single threshold and throw a warning).

The most likely culprit for the terminating error would be the log scale (lux?) of the Migrate Tech tags.

The first attempt to resolve the issue would be to convert the Migrate Tech light data to a linear scale, I suspect there is an issue with the estimated threshold. Try,

d.lig$Light <- exp(d.lig$Light)

Then calibrate again. It's desirable to add a small margin of error to the threshold, but if the data are on a log scale then a "small" margin won't even make a dent.

The algorithm failed on the second day of the backwards path, so it's possible that the data from the last day were not at the retrieved.at location. This could be because you have trimmed the data or because the tag failed before it was retrieved. The simplest way to try and resolve this is to remove the retrieved.at parameter and let the algorithm find it from the data alone.

I've not tested the method on Migrate Tech data, but my intuition tells me that the log scale is wreaking havoc on estimates of the threshold while the animal is at sea.

VirginiaMorera commented 6 years ago

Hello. So, I tried your suggestions, and here is the outcome:

First, I don't think the light data is in a log scale. It can be represented in a log scale when you use the dedicated software (Intiproc), but it is not stored like that in the tag. The difference with the light recordings from BAS or Biotrack loggers is that (at least in the mode we use) the light is not clipped, so instead of having values between 0 and 64 it can have much higher values. In this particular case, range(d.lig$Light) gave [1] 1.136 74418.605. So, what actually would make the light measurements more similar to those of a Biotrack tag would be log(d.lig$Light) I guess. But, I've tried a couple of trips both with and without transforming and the outcome is very similar, so I am happy to report havoc is not being wreaked at all.

Regarding the failure to run the model, I only managed to make it work after removing both deployed.at and recovered.at from the specifications, even though I am sure of the coordinates as I have selected a complete track to test this. My guess is that the problem has to do with the land mask, since the loggers where deployed and recovered obviously on land... Unfortunately, without those two parameters, the model fails to correctly estimate the locations of the first and last positions. In the image, the yellow dot is the deployment-recovery location, the black one is the first position and the red one the last (as estimated by the model). I guess an easy workaround would be to specify as deployment and recovery the closest at-sea location from the actual place, if you think that might be the problem. first_last_location

Lastly, I have tried a couple of trips that I have already analysed with the regular method of manually adjusting transitions, for comparison. Although the tracks look pretty alike, to me it looks a bit weird that it estimates the stationary positions during winter as a kind of "path" with very short steps, but with a directionality, instead of random movements in all directions. I don't know if that might be related to the movement hyperparameters I give to the model, and what changes could improve that. imagen imagen These two are Migrate tags, so it seems to be working pretty well!

I suspect the positions in the Mediterranean of the second track, which are definitely not real, have to do with the animal spending days incubating inside the burrow (around june-july in the plot below), but I don't know if such a strong effect can be fixed changing the shading hyperparameters without affecting the rest of the track...

imagen

Anyway, great job and a great improvement over the traditional process of manually adjusting transitions!

ABindoff commented 6 years ago

Apologies for the confusion, I did mean try a log transformation (sometime the words don't match the pictures in my head!) I had a look at the specs for these tags and they are definitely logarithmic. This means small differences in the light (say, the difference between cloudy days during calibration and clear days during deployment) will be represented by large differences in the data.

I can see from the thresholded light trace (last figure) that you have light where it should be dark around late Nov, and some more around July. If a log transformation doesn't give you a better threshold, then choose a larger threshold (until that errant above-threshold light disappears). This will definitely cause severe errors (not just on the affected days).

I suspect you're right about the deployed.at and retrieved.at locations falling inside the land mask. I've never had this problem myself, but it makes sense and it's easy to fix by changing the locations a little to put them off the coast. I had always intended to allow more flexible priors on land/sea masks (for a number of reasons).

Stationary positions

Some good news here - there is a very good solution. But first I'll explain the problem. Because the algorithm essentially has very poor information about location on those days what you'll get is the mathematically plausible but ecologically implausible solution - a solution is found using whatever light is available and whatever locations are estimated either side (hence the drift).

The solution is straightforward and greatly improves estimates if you have good reason to believe the animal is staying in the one place (seems obvious from your light trace, and may be supported by known information about the animal), and you know where that place is.

Instead of using deployed.at and retrieved.at, used fixd as per the example given here, https://github.com/ABindoff/TwilightFree/blob/master/tutorials/Using_SST_and_sightings.Rmd This will force the estimated positions to stay put at the nesting location. You just have to go through your light data and figure out when those days were.

If you're less confident that the animal was nesting at a particular location, you can set the movement hyperparameter to c(1, 1/2) or c(1, 3/5) or something like that. What this will do is make the scale parameter of the assumed Gamma distribution high. This means you tell the model you have more reason to believe that the animal stayed within a cell over consecutive days, because more of the probability distribution is at 0. It seems fairly clear that your animals are nesting at a known nesting site on particular days, though, so supplying that information through the fixd parameter is a much more reasonable solution.

ABindoff commented 6 years ago

Having a think about my previous reply and the issue of the deployed.at and retrieved.at locations possibly being on land, try the following:

grid <- makeGrid(c(-62, 80), c(-47, 70), cell.size = 1, pacific = F) grid[!grid] <- 1E-10

This will allow travel over land, but make it much less likely than travel over sea.