2.2 Estimating clock drift
For correcting clock drift, we need to first provide some reasonably good guesses of the true drift values.
find_clock_drift() can take care of this.
This function takes a brute force approach and uses cross correlation to see how closely sync-tag detections and emissions line up over a time series of possible clock drift values.
This works well for getting values on the scale of, say, 0.1 second resolution.
For each clock_time, we’ll check possible drift values within a 1 hour-period
int_drft <- find_clock_drift(
sync_model, clock_resolution = 0.1, clock_window = 60*60*1, quiet = T)
print(int_drft)## Receiver pair: 62074→62368
## Drift (s) xcorr
## 2021-05-08 18:00:00 34.7 0.92
## 2021-05-08 19:00:00 34.8 0.92
## 2021-05-08 20:00:00 34.9 0.67
## 2021-05-08 21:00:00 34.9 0.82
## 2021-05-08 22:00:00 35.0 0.83
## 2021-05-08 23:00:00 35.1 0.88
## 2021-05-09 00:00:00 35.2 0.67
## 2021-05-09 01:00:00 35.3 0.78
## 2021-05-09 02:00:00 35.4 0.88
## 2021-05-09 03:00:00 35.4 0.87
## 2021-05-09 04:00:00 35.5 0.87
## 2021-05-09 05:00:00 35.6 0.60
## 2021-05-09 06:00:00 35.6 0.53
##
## Receiver pair: 62074→62581
## Drift (s) xcorr
## 2021-05-08 18:00:00 30.2 0.93
## 2021-05-08 19:00:00 30.3 0.92
## 2021-05-08 20:00:00 30.4 0.67
## 2021-05-08 21:00:00 30.4 0.82
## 2021-05-08 22:00:00 30.5 0.83
## 2021-05-08 23:00:00 30.6 0.88
## 2021-05-09 00:00:00 30.7 0.67
## 2021-05-09 01:00:00 30.7 0.73
## 2021-05-09 02:00:00 30.8 0.88
## 2021-05-09 03:00:00 30.9 0.87
## 2021-05-09 04:00:00 30.9 0.78
## 2021-05-09 05:00:00 31.0 0.67
## 2021-05-09 06:00:00 31.0 0.53
##
## Receiver pair: 62074→62041
## Drift (s) xcorr
## 2021-05-08 18:00:00 0.5 0.83
## 2021-05-08 19:00:00 0.5 0.92
## 2021-05-08 20:00:00 0.5 0.67
## 2021-05-08 21:00:00 0.5 0.82
## 2021-05-08 22:00:00 0.5 0.88
## 2021-05-08 23:00:00 0.5 0.92
## 2021-05-09 00:00:00 0.5 0.88
## 2021-05-09 01:00:00 0.5 0.78
## 2021-05-09 02:00:00 0.5 0.83
## 2021-05-09 03:00:00 0.6 0.87
## 2021-05-09 04:00:00 0.6 0.78
## 2021-05-09 05:00:00 0.6 0.82
## 2021-05-09 06:00:00 0.6 0.57
##
## Receiver pair: 62074→62042
## Drift (s) xcorr
## 2021-05-08 18:00:00 33.3 0.93
## 2021-05-08 19:00:00 33.4 0.92
## 2021-05-08 20:00:00 33.5 0.67
## 2021-05-08 21:00:00 33.5 0.82
## 2021-05-08 22:00:00 33.6 0.88
## 2021-05-08 23:00:00 33.7 0.78
## 2021-05-09 00:00:00 33.8 0.88
## 2021-05-09 01:00:00 33.9 0.83
## 2021-05-09 02:00:00 34.0 0.88
## 2021-05-09 03:00:00 34.0 0.88
## 2021-05-09 04:00:00 34.1 0.87
## 2021-05-09 05:00:00 34.2 0.73
## 2021-05-09 06:00:00 34.2 0.65
##
## Receiver pair: 62368→62584
## Drift (s) xcorr
## 2021-05-08 18:00:00 4.4 0.92
## 2021-05-08 19:00:00 4.4 0.88
## 2021-05-08 20:00:00 4.4 0.88
## 2021-05-08 21:00:00 4.4 0.96
## 2021-05-08 22:00:00 4.5 0.83
## 2021-05-08 23:00:00 4.4 0.83
## 2021-05-09 00:00:00 4.5 0.96
## 2021-05-09 01:00:00 4.5 0.83
## 2021-05-09 02:00:00 4.5 0.91
## 2021-05-09 03:00:00 4.5 0.88
## 2021-05-09 04:00:00 4.5 0.84
## 2021-05-09 05:00:00 4.5 0.92
## 2021-05-09 06:00:00 4.5 0.92To run the cross correlation, find_clock_drift first creates a discrete time series of drift values with a resolution given by clock_resolution.
For a receiver pair, there is one time series for detections and another for emissions.
clock_window determines how long these time series are (here, the series each span 1 hour).
Next, detections & emissions are marked as ones while all other values in the series are left at zero.
These series are then convolved with a Gaussian kernel (the standard deviation can be specified with the argument sd).
The convolution causes each detection & emission to now be represented by a Gaussian shaped curve in each time series.
Finally, the pair of time series are cross correlated with each other.
For each possible drift value, a xcorr value is returned, giving a standardized cross correlation score.
Drift values with higher xcorr indicate the pair of time series match up well when that particular clock drift correction is applied to it.
For each receiver pair and clock time, the clock drift value which gave the best xcorr score is returned.
For datasets with more extreme clock drift, its important that sd is set sufficiently high.
Setting a higher sd value causes sync-tag emission/detection values that are almost aligned to give higher xcorr scores.
Lastly, make sure to set clock_window large enough so that you have quite a lot of sync-tag detections within that period.
Keep in mind that sometimes periods of high background noise can cause ‘dark zones’ in your detections data.
Our values look good, so we’ll go ahead and set them in our model.
# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- int_drft
# Plot sync-tag diagnostics
plot(sync_model, layout = c(4,4))
As you can see in the plot, our clocks are roughly corrected within a 0.1 second resolution. For fine-scale positioning, we need to refine this further.
fit_clock_drift will automatically refine these clock drift values.
For each receiver pair, the clock drift is treated a a Gaussian random walk along some linear trend line.
The drift values are fitted as random effects in the R package TMB and the detection error is treated as a mixture distribution (see ?ddetect).
This function can provide extremely high resolution drift values.
However, for the model to fit, it requires quite good starting values.
We suggest using find_clock_drift to get drift values with at least a 0.1 second resolution before applying fit_clock_drift.
## --- Kaltoa receiver drift model
## Detection sigma (ms) Std. Err. Drift (ms h⁻¹) Std. Err.
## 62074→62368 0.036527 NaN 82.284 NaN
## 62074→62581 0.065194 NaN 64.283 NaN
## 62074→62041 0.048325 NaN 10.607 NaN
## 62074→62042 0.16552 NaN 77.527 NaN
## 62368→62584 1.5826 NaN 9.9174 NaN
## Drift sigma (ms) Std. Err.
## 62074→62368 177.64 NaN
## 62074→62581 148.21 NaN
## 62074→62041 0.033131 NaN
## 62074→62042 137.3 NaN
## 62368→62584 3.0461e-06 NaN
##
## Drift values
## 62074→62368 62074→62581 62074→62041 62074→62042 62368→62584
## 2021-05-08 18:00:00 34.583 30.106 0.37325 33.210 4.3139
## 2021-05-08 19:00:00 34.668 30.180 0.38388 33.290 4.3238
## 2021-05-08 20:00:00 34.753 30.254 0.39454 33.370 4.3338
## 2021-05-08 21:00:00 34.838 30.328 0.40516 33.450 4.3437
## 2021-05-08 22:00:00 34.922 30.402 0.41577 33.531 4.3536
## 2021-05-08 23:00:00 35.007 30.476 0.42638 33.611 4.3635
## 2021-05-09 00:00:00 35.092 30.550 0.43701 33.691 4.3734
## 2021-05-09 01:00:00 35.177 30.624 0.44761 33.772 4.3833
## 2021-05-09 02:00:00 35.262 30.698 0.45821 33.852 4.3933
## 2021-05-09 03:00:00 35.347 30.772 0.46882 33.932 4.4032
## 2021-05-09 04:00:00 35.431 30.846 0.47941 34.013 4.4131
## 2021-05-09 05:00:00 35.516 30.919 0.48997 34.093 4.4230
## 2021-05-09 06:00:00 35.601 30.993 0.50054 34.173 4.4329# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- fit_drft
# Plot sync-tag diagnostics
plot(sync_model, type = 'error', layout = c(4,4))
Now the sync-tag diagnostic plots show that the transmission latency is more-less constant over time. Note how some receiver pairs have a bias above or below their expected transmission latencies. These can be correct for in the next section.