3.3 Estimating clock drift: Mixed-model
fit_clock_drift will estimate clock drift using a mixed model implemented in TMB.
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 and the detection error is treated as a mixture distribution (see ?ddetect).
For more details on the model, see Campbell et al. (in prep).
This function can provide extremely high resolution drift values.
However, for the model to fit, it requires quite good starting values.
We suggest using permute_clock_drift to get drift values with at least a 0.1 second resolution before applying fit_clock_drift.
fit_drft_1 <- fit_clock_drift(
sync_model, quiet = T,
sigma_det = 0.1)# Use clock_resolution as initial detection error value.
print(fit_drft_1)## --- Kaltoa receiver drift model
## Detection σ (ms) | Drift rate (ms h⁻¹) | Drift σ (ms h⁻¹) | TOA offset (ms) | TDOA offset σ (ms)
## Value, SE | Value, SE | Value, SE | Value, SE | Value, SE
## 3→1 0.10000, 0.0044737 | -6.0597, 10.7637 | 37.286, 7.6110 | 0.86029, 2.6137e-01 | 1.8955, 0.67323
## 3→2 0.33702, 0.0108150 | 5.2948, 12.5508 | 43.477, 8.8750 | -1.19639, 9.4900e-02 | 1.2771, 0.45356
## 3→4 0.26161, 0.0095992 | 9.8684, 12.5214 | 43.375, 8.8540 | 1.35560, 2.2941e-01 | 3.1968, 1.13045
## 3→5 0.28389, 0.0126186 | 7.2387, 13.5897 | 47.076, 9.6095 | -104.98715, 8.1944e+05 | 2.1594, 0.88195
## 4→6 0.89992, 0.0337254 | 1.6084, 7.9666 | 27.597, 5.6343 | -3.16099, 6.1703e-02 | 4.4181, 1.56872
##
## Drift values
## 3→1 3→2 3→4 3→5 4→6
## 2021-05-08 18:00:00 0.37241 33.210 34.583 30.108 4.3173
## 2021-05-08 19:00:00 0.38302 33.291 34.667 30.181 4.3273
## 2021-05-08 20:00:00 0.39369 33.372 34.751 30.255 4.3371
## 2021-05-08 21:00:00 0.40433 33.451 34.836 30.329 4.3467
## 2021-05-08 22:00:00 0.41496 33.532 34.921 30.403 4.3571
## 2021-05-08 23:00:00 0.42556 33.612 35.006 30.477 4.3666
## 2021-05-09 00:00:00 0.43615 33.693 35.091 30.551 4.3771
## 2021-05-09 01:00:00 0.44676 33.774 35.176 30.624 4.3861
## 2021-05-09 02:00:00 0.45733 33.853 35.260 30.699 4.3962
## 2021-05-09 03:00:00 0.46797 33.933 35.345 30.773 4.4059
## 2021-05-09 04:00:00 0.47858 34.014 35.430 30.846 4.4157
## 2021-05-09 05:00:00 0.48907 34.095 35.515 30.920 4.4260
## 2021-05-09 06:00:00 0.49969 34.174 35.601 30.995 4.4367
sigma_det is the initial value for the detection error standard deviation (for direct transmissions), in seconds.
It’s recommended to set this to the clock_resolution you used to create the initial drift values.
Setting this too small can cause the model fail.
Fitting these drift values can be a messy process.
There are many parameters being estimated simultaneously by the model, so we often run into problems such as clock drift being mistaken as large detection error values.
To work around this you can opt to fit the model while keeping the inital parameter values fixed.
fix_sigma_det = T tells our model to treat our provided detection error parameter sigma_det = 0.01 as the truth rather than try to estimate it.
Below, we’ll refit the clock drift model a few times using progressively smaller drift error values that are set as fixed values.
sync_model_tmp <- sync_model
fit_drft_tmp <- fit_clock_drift(
sync_model_tmp, quiet = T,
sigma_det = 0.1,
fix_sigma_det = T)# Do not fit the detection error
print(fit_drft_tmp)## --- Kaltoa receiver drift model
## Detection σ (ms) | Drift rate (ms h⁻¹) | Drift σ (ms h⁻¹) | TOA offset (ms) | TDOA offset σ (ms)
## Value, SE | Value, SE | Value, SE | Value, SE | Value, SE
## 3→1 100, 0 | -6.8306, 10.0036 | 34.585, 10.9434 | -29.547, 4.4841 | 0.0026706, 6.3765
## 3→2 100, 0 | 6.8892, 10.8072 | 37.369, 13.2410 | -32.885, 4.9385 | 0.0023298, 6.4655
## 3→4 100, 0 | 7.5229, 11.4339 | 37.787, 13.9303 | -30.634, 5.2949 | 0.0026376, 6.6525
## 3→5 100, 0 | 5.2222, 11.3071 | 36.840, 17.9434 | -34.853, 6.2695 | 0.0022849, 6.4820
## 4→6 100, 0 | 1.8194, 7.7189 | 27.131, 8.9827 | -32.870, 4.3027 | 0.0031275, 8.8260
##
## Drift values
## 3→1 3→2 3→4 3→5 4→6
## 2021-05-08 18:00:00 0.37919 33.203 34.603 30.130 4.3232
## 2021-05-08 19:00:00 0.38298 33.308 34.650 30.163 4.3290
## 2021-05-08 20:00:00 0.37513 33.354 34.749 30.252 4.3385
## 2021-05-08 21:00:00 0.42149 33.450 34.835 30.328 4.3486
## 2021-05-08 22:00:00 0.41674 33.533 34.921 30.405 4.3583
## 2021-05-08 23:00:00 0.42626 33.613 35.006 30.496 4.3670
## 2021-05-09 00:00:00 0.43695 33.710 35.108 30.532 4.3743
## 2021-05-09 01:00:00 0.44764 33.754 35.157 30.626 4.3663
## 2021-05-09 02:00:00 0.45792 33.854 35.260 30.719 4.4168
## 2021-05-09 03:00:00 0.46763 33.936 35.346 30.754 4.4103
## 2021-05-09 04:00:00 0.45934 34.034 35.433 30.845 4.4177
## 2021-05-09 05:00:00 0.50787 34.077 35.534 30.941 4.4258
## 2021-05-09 06:00:00 0.49320 34.184 35.570 30.965 4.4313# Store fitted drift values
clock_drift(sync_model_tmp) <- fit_drft_tmp
# Refit with a smaller detection error
fit_drft_tmp <- fit_clock_drift(
sync_model_tmp, quiet = T,
sigma_det = 0.01,
fix_sigma_det = T)
print(fit_drft_tmp)## --- Kaltoa receiver drift model
## Detection σ (ms) | Drift rate (ms h⁻¹) | Drift σ (ms h⁻¹) | TOA offset (ms) | TDOA offset σ (ms)
## Value, SE | Value, SE | Value, SE | Value, SE | Value, SE
## 3→1 10, 0 | 0.93122, 5.1364 | 17.735, 3.9530 | 0.31031, 0.67380 | 1.13183592, 1.09533
## 3→2 10, 0 | -1.23499, 6.4926 | 22.446, 4.8565 | -1.69568, 0.68332 | 0.00081611, 0.90157
## 3→4 10, 0 | 4.05460, 6.2017 | 21.426, 4.6826 | 0.78371, 0.69283 | 2.98480165, 1.33102
## 3→5 10, 0 | 4.14128, 7.7305 | 26.739, 5.6461 | -1.32099, 0.69064 | 1.14739693, 1.07107
## 4→6 10, 0 | 0.84612, 3.5597 | 12.238, 2.9203 | -3.72635, 0.65620 | 3.04548050, 1.91303
##
## Drift values
## 3→1 3→2 3→4 3→5 4→6
## 2021-05-08 18:00:00 0.37258 33.210 34.583 30.108 4.3181
## 2021-05-08 19:00:00 0.38369 33.292 34.666 30.181 4.3274
## 2021-05-08 20:00:00 0.39297 33.371 34.752 30.255 4.3375
## 2021-05-08 21:00:00 0.40561 33.452 34.836 30.329 4.3474
## 2021-05-08 22:00:00 0.41485 33.533 34.921 30.403 4.3576
## 2021-05-08 23:00:00 0.42600 33.612 35.005 30.478 4.3669
## 2021-05-09 00:00:00 0.43647 33.694 35.092 30.551 4.3784
## 2021-05-09 01:00:00 0.44708 33.773 35.175 30.625 4.3839
## 2021-05-09 02:00:00 0.45738 33.854 35.261 30.699 4.3994
## 2021-05-09 03:00:00 0.46894 33.934 35.345 30.772 4.4052
## 2021-05-09 04:00:00 0.47747 34.015 35.429 30.846 4.4161
## 2021-05-09 05:00:00 0.49127 34.094 35.516 30.921 4.4264
## 2021-05-09 06:00:00 0.49779 34.176 35.598 30.993 4.4363# Store fitted drift values
clock_drift(sync_model_tmp) <- fit_drft_tmp
# Refit with a smaller detection error
fit_drft_tmp <- fit_clock_drift(
sync_model_tmp, quiet = T,
sigma_det = 0.001,
fix_sigma_det = T)
print(fit_drft_tmp)## --- Kaltoa receiver drift model
## Detection σ (ms) | Drift rate (ms h⁻¹) | Drift σ (ms h⁻¹) | TOA offset (ms) | TDOA offset σ (ms)
## Value, SE | Value, SE | Value, SE | Value, SE | Value, SE
## 3→1 1, 0 | 0.14762, 0.47908 | 1.65349, 0.38606 | 0.86394, 0.066412 | 1.8774, 0.66822
## 3→2 1, 0 | -0.13275, 0.33397 | 1.14846, 0.29472 | -1.13409, 0.068088 | 1.3346, 0.47965
## 3→4 1, 0 | 0.21968, 0.37926 | 1.30454, 0.32731 | 1.35383, 0.068259 | 3.2232, 1.14205
## 3→5 1, 0 | 0.10335, 0.23256 | 0.79243, 0.23528 | -0.77206, 0.068083 | 1.7831, 0.63446
## 4→6 1, 0 | 0.15260, 0.57211 | 1.97216, 0.45521 | -14.25550, 2.446931 | 7.4018, 3.52446
##
## Drift values
## 3→1 3→2 3→4 3→5 4→6
## 2021-05-08 18:00:00 0.37237 33.211 34.582 30.107 4.3283
## 2021-05-08 19:00:00 0.38302 33.291 34.666 30.181 4.3382
## 2021-05-08 20:00:00 0.39361 33.371 34.752 30.255 4.3481
## 2021-05-08 21:00:00 0.40440 33.452 34.836 30.329 4.3577
## 2021-05-08 22:00:00 0.41493 33.532 34.921 30.403 4.3683
## 2021-05-08 23:00:00 0.42561 33.612 35.006 30.477 4.3777
## 2021-05-09 00:00:00 0.43619 33.693 35.091 30.551 4.3885
## 2021-05-09 01:00:00 0.44678 33.773 35.175 30.625 4.3970
## 2021-05-09 02:00:00 0.45725 33.853 35.260 30.699 4.4077
## 2021-05-09 03:00:00 0.46805 33.934 35.345 30.772 4.4166
## 2021-05-09 04:00:00 0.47836 34.014 35.430 30.846 4.4269
## 2021-05-09 05:00:00 0.48931 34.094 35.515 30.920 4.4372
## 2021-05-09 06:00:00 0.49934 34.175 35.599 30.993 4.4483
The above method of iteratively fitting smaller fixed values should be seen as a last resort. It’s a messy way to get rough drift values when the model otherwise won’t fit.
Every telemetry array is unique, so some trial and error will have to be used to find which fitting approach works best for the dataset on hand.
Finally, lets store our fitted drift values in our KaltoaClockSync model.
# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- fit_drft_1
# 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 corrected for in the next section.