ISF test
Morera-Pujol, Virginia
May 2021
Individual site fidelity test
Testing the IndEffectTest function developed for the Detecting recurrent sources of variability in animal tracking studies manuscript by Virginia Morera-Pujol et al. submitted to Ecological Applications in 2021
Description
In this document I’ll generate, using the package SIMRiv, simulated animal movement tracks to test the detection of Individual Site Fidelity by the function developed for the MS “Detecting recurrent sources of variability in animal tracking studies”, by Virginia Morera-Pujol et al. (2020). The details of how the tracks with ISF can be simulated can be found here
Here, we’ll run the function once in a dataset of 10 individuals with 10 trips each, with different origins, and once on a dataset of 10 individuals with 10 trips each with the same origin (central place forager scenario).
For animals with different origins
Generating the data
We use a value for resistance = 1, which means we’re simulating maximum ISF (simulated tracks are forced within the template area of no resistance because resistance outside the template area is 1).
'%!in%' <- function(x,y)!('%in%'(x,y))
set.seed(123)
# set movement parameters of the "species".
# We're going to configure a Correlated Random walker
CRW <- species(state.CRW(0.3))
test <- simulate(CRW, coords = c(10, 15), time = 100)
# create empty base raster for resistance surfaces
baseRaster <- raster(xmn = -50, xmx = 50, ymn = -50, ymx = 50, res = 0.1)
values(baseRaster) <- 1
baseGrid <- as(baseRaster, "SpatialPixels")
baseGrid@proj4string <- CRS(as.character(NA))
# generate 10 random origin points within a 100000x100000 square, one per "individual"
origin_points <- data.frame(x = runif(10, -25, 25),
y = runif(10, -25, 25))
# simulate 10 trips for each origin (i.e. 10 trips per individual)
template_tracks <- list()
template_UDs <- list()
resistance_surfaces <- list()
simulated_tracks <- list()
for (i in 1:10) {
# i = 2
# First we need to create a template trip to get the UD for the resistance surface
template <- simulate(CRW, coords = c(origin_points[i,1], origin_points[i, 2]),
time = 100)
template_tracks[i] <- list(template)
# plot(template, col = brewer.pal(10,"Set3")[i], type = "l")
# points(template, col = brewer.pal(10, "Set3")[i], pch = 20)
spTemplate <- SpatialPoints(coords = template[,1:2], proj4string = baseGrid@proj4string)
KDE <- kernelUD(spTemplate, h = "href", grid = baseGrid, extent = 2)
UD50 <- getverticeshr(KDE, percent = 50)
template_UDs[i] <- list(UD50)
# plot(KDE)
# plot(UD50, add = T)
# plot(spTemplate, pch = 20, cex = 0.5, add = T)
# Second, we create the resistance surface
res_surface <- resistanceFromShape(UD50, baseRaster = baseRaster, field = 0,
background = 1)
resistance_surfaces[i] <- list(res_surface)
# plot(res_surface)
# Third, we simulate 10 tracks from each surface (i.e. from each individual)
sim_tracks <- simulate(list(CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW),
start.resistance = 0,
time = 100, resist = res_surface)
spSim <- SpatialPointsDataFrame(coords = cbind(c(sim_tracks[,1], sim_tracks[,4],
sim_tracks[,7], sim_tracks[,10],
sim_tracks[,13],sim_tracks[,16],
sim_tracks[,19], sim_tracks[,22],
sim_tracks[,25], sim_tracks[,28]),
c(sim_tracks[,2], sim_tracks[,5],
sim_tracks[,8], sim_tracks[,11],
sim_tracks[,14], sim_tracks[,17],
sim_tracks[,20], sim_tracks[,23],
sim_tracks[,26], sim_tracks[,29])),
data = data.frame(TrackID = paste(i,
rep(1:10, each = 100),
sep = "_"),
BirdID = rep(i, 1000)))
spSim@proj4string <- baseGrid@proj4string
simulated_tracks[i] <- list(spSim)
spSim$BirdID <- NULL
KUD_sim <- kernelUD(spSim, grid = baseGrid)
}Now the list simulated_tracks has 10 elements, each of them containing the 10 tracks belonging to each “individual”. We need to reformat that in a way that we can proceed with the Individual Effect Test
list.df <- lapply(seq_along(simulated_tracks),
function(i){simulated_tracks[i][[1]]})
# apply rbind to the list
simulated_dataset <- do.call("rbind",
list.df)
# assign null projection as X/Y values don't mean anything in this case
simulated_dataset@proj4string <- CRS(as.character(NA))Visualising the data
# first create a version of the data without BirdID to calculate trip overlaps
sds <- as.data.frame(simulated_dataset)
# convert it in order to be able to plot it
sds %>%
st_as_sf(coords = c("coords.x1", "coords.x2")) %>%
tidyr::separate(TrackID, into = c("BirdID", "TrackID"), sep = "_") -> sdf
col1 <- ggplot() +
geom_sf(data = sdf, aes(col = BirdID), alpha = 0.3, size = 0.7) +
xlab("Longitude") + ylab("Latitude") +
# scale_shape_manual(values = 1:10, guide = "none") +
theme_bw() +
theme(legend.position = "right", legend.box = "horizontal") +
guides(col = guide_legend(override.aes = list(size = 3, alpha = 1))) +
# scale_shape(guide = "none") +
NULL
print(col1)
# ggsave(filename = "HighISF.png", width = 12*1.5, height = 10*1.5, dpi = 800, units = "cm")Testing the function
We now run a simplified version of the function IndEffectTest without all the code re the spatial projection, crossing datelines, etc. and as plain code outside the function for easy tweaking.
Calculate akde for all tracks and pairwise overlaps
# To calculate the overlap we:
# 1 - Turn simulated tracks into "telemetry" objects"
# 2 - fit a guesstimate movement model to each track
# 3 - using the parameters of the guesstimate, fit a series of competing movement models
# 4 - Calculate the debiased BA between top-selected models for each track
# (uncertainty from the model fits is propagated into the overlap estimate under the approximation that the BA is a chi-square random variable)
# 1
sds %>%
dplyr::group_by(TrackID) %>%
dplyr::mutate(Date_Time = seq(ymd_hms('2019-11-30 12:10:00'), ymd_hms('2020-01-30 12:10:00'),
length.out = 100)) %>%
mutate(Date_Time = format(Date_Time, format = "%Y/%m/%d %H:%M:%S", tz = "UTC")) %>%
dplyr::select(individual.local.identifier = TrackID,
Longitude = coords.x1, Latitude = coords.x2,
datetime = Date_Time) %>%
ungroup() -> sds
sds <- as.data.frame(sds)
sdst <- sds %>%
filter(str_detect(individual.local.identifier, "^7"))
sd_tel <- as.telemetry(object = sds, timeformat = "%Y/%m/%d %H:%M:%S")
FIT.list <- list()
file <- file("FIT.list_nonColonial.txt", open = "wt")
sink(file, type = "message") # sink text messages to .txt to avoid them appearing on the markdown
before <- Sys.time()
for(i in 1:length(sd_tel)){
print(i)
one <- sd_tel[[i]]
# 2
GUESS <- ctmm.guess(one, interactive=F)
# 3
FITS <- ctmm.select(one, GUESS, trace=TRUE, verbose=TRUE, cores=-1)
print(paste("The selected model for", names(sd_tel)[i], "is", names(FITS)[1], sep = " "))
FIT.list[[i]] <- FITS[[1]]
}## [1] 1
## [1] "The selected model for 1_1 is OUF isotropic"
## [1] 2
## [1] "The selected model for 1_10 is OUF isotropic"
## [1] 3
## [1] "The selected model for 1_2 is OUF isotropic"
## [1] 4
## [1] "The selected model for 1_3 is OUF isotropic"
## [1] 5
## [1] "The selected model for 1_4 is OUF isotropic"
## [1] 6
## [1] "The selected model for 1_5 is OUF isotropic"
## [1] 7
## [1] "The selected model for 1_6 is OUF isotropic"
## [1] 8
## [1] "The selected model for 1_7 is OUF isotropic"
## [1] 9
## [1] "The selected model for 1_8 is OUF isotropic"
## [1] 10
## [1] "The selected model for 1_9 is OUF isotropic"
## [1] 11
## [1] "The selected model for 10_1 is OUF isotropic"
## [1] 12
## [1] "The selected model for 10_10 is OUF isotropic"
## [1] 13
## [1] "The selected model for 10_2 is OUF isotropic"
## [1] 14
## [1] "The selected model for 10_3 is OUF isotropic"
## [1] 15
## [1] "The selected model for 10_4 is OUF isotropic"
## [1] 16
## [1] "The selected model for 10_5 is OUF isotropic"
## [1] 17
## [1] "The selected model for 10_6 is OUF isotropic"
## [1] 18
## [1] "The selected model for 10_7 is OUF isotropic"
## [1] 19
## [1] "The selected model for 10_8 is OU isotropic"
## [1] 20
## [1] "The selected model for 10_9 is OUF isotropic"
## [1] 21
## [1] "The selected model for 2_1 is OUF isotropic"
## [1] 22
## [1] "The selected model for 2_10 is OUF isotropic"
## [1] 23
## [1] "The selected model for 2_2 is OUF isotropic"
## [1] 24
## [1] "The selected model for 2_3 is OUF isotropic"
## [1] 25
## [1] "The selected model for 2_4 is OUF isotropic"
## [1] 26
## [1] "The selected model for 2_5 is OUF isotropic"
## [1] 27
## [1] "The selected model for 2_6 is OUF isotropic"
## [1] 28
## [1] "The selected model for 2_7 is OUF isotropic"
## [1] 29
## [1] "The selected model for 2_8 is OUF isotropic"
## [1] 30
## [1] "The selected model for 2_9 is OUF isotropic"
## [1] 31
## [1] "The selected model for 3_1 is OUF isotropic"
## [1] 32
## [1] "The selected model for 3_10 is OUF isotropic"
## [1] 33
## [1] "The selected model for 3_2 is OUF isotropic"
## [1] 34
## [1] "The selected model for 3_3 is OUF isotropic"
## [1] 35
## [1] "The selected model for 3_4 is OUF isotropic"
## [1] 36
## [1] "The selected model for 3_5 is OUF isotropic"
## [1] 37
## [1] "The selected model for 3_6 is OUF isotropic"
## [1] 38
## [1] "The selected model for 3_7 is OUF isotropic"
## [1] 39
## [1] "The selected model for 3_8 is OUF isotropic"
## [1] 40
## [1] "The selected model for 3_9 is OUF isotropic"
## [1] 41
## [1] "The selected model for 4_1 is OUF isotropic"
## [1] 42
## [1] "The selected model for 4_10 is OUF anisotropic"
## [1] 43
## [1] "The selected model for 4_2 is OUF isotropic"
## [1] 44
## [1] "The selected model for 4_3 is OUF isotropic"
## [1] 45
## [1] "The selected model for 4_4 is OUF anisotropic"
## [1] 46
## [1] "The selected model for 4_5 is OU anisotropic"
## [1] 47
## [1] "The selected model for 4_6 is OUF isotropic"
## [1] 48
## [1] "The selected model for 4_7 is OUF anisotropic"
## [1] 49
## [1] "The selected model for 4_8 is OUF isotropic"
## [1] 50
## [1] "The selected model for 4_9 is OUF isotropic"
## [1] 51
## [1] "The selected model for 5_1 is OUF isotropic"
## [1] 52
## [1] "The selected model for 5_10 is OUF isotropic"
## [1] 53
## [1] "The selected model for 5_2 is OUF isotropic"
## [1] 54
## [1] "The selected model for 5_3 is OUF isotropic"
## [1] 55
## [1] "The selected model for 5_4 is OUF isotropic"
## [1] 56
## [1] "The selected model for 5_5 is OUF isotropic"
## [1] 57
## [1] "The selected model for 5_6 is OUF isotropic"
## [1] 58
## [1] "The selected model for 5_7 is OUF isotropic"
## [1] 59
## [1] "The selected model for 5_8 is OUF isotropic"
## [1] 60
## [1] "The selected model for 5_9 is OUF anisotropic"
## [1] 61
## [1] "The selected model for 6_1 is OUF isotropic"
## [1] 62
## [1] "The selected model for 6_10 is OUF isotropic"
## [1] 63
## [1] "The selected model for 6_2 is OUF isotropic"
## [1] 64
## [1] "The selected model for 6_3 is OUF isotropic"
## [1] 65
## [1] "The selected model for 6_4 is OUF isotropic"
## [1] 66
## [1] "The selected model for 6_5 is OUF isotropic"
## [1] 67
## [1] "The selected model for 6_6 is OUF anisotropic"
## [1] 68
## [1] "The selected model for 6_7 is OUF isotropic"
## [1] 69
## [1] "The selected model for 6_8 is OUF isotropic"
## [1] 70
## [1] "The selected model for 6_9 is OUF isotropic"
## [1] 71
## [1] "The selected model for 7_1 is OUF isotropic"
## [1] 72
## [1] "The selected model for 7_10 is OUF isotropic"
## [1] 73
## [1] "The selected model for 7_2 is OUF isotropic"
## [1] 74
## [1] "The selected model for 7_3 is OUF isotropic"
## [1] 75
## [1] "The selected model for 7_4 is OUF isotropic"
## [1] 76
## [1] "The selected model for 7_5 is OUF isotropic"
## [1] 77
## [1] "The selected model for 7_6 is OUF anisotropic"
## [1] 78
## [1] "The selected model for 7_7 is OUF isotropic"
## [1] 79
## [1] "The selected model for 7_8 is OUF isotropic"
## [1] 80
## [1] "The selected model for 7_9 is OUF isotropic"
## [1] 81
## [1] "The selected model for 8_1 is OUF isotropic"
## [1] 82
## [1] "The selected model for 8_10 is OUF isotropic"
## [1] 83
## [1] "The selected model for 8_2 is OUF isotropic"
## [1] 84
## [1] "The selected model for 8_3 is OUF isotropic"
## [1] 85
## [1] "The selected model for 8_4 is OUF isotropic"
## [1] 86
## [1] "The selected model for 8_5 is OUF isotropic"
## [1] 87
## [1] "The selected model for 8_6 is OUF isotropic"
## [1] 88
## [1] "The selected model for 8_7 is OUF isotropic"
## [1] 89
## [1] "The selected model for 8_8 is OUF isotropic"
## [1] 90
## [1] "The selected model for 8_9 is OUF isotropic"
## [1] 91
## [1] "The selected model for 9_1 is OUF isotropic"
## [1] 92
## [1] "The selected model for 9_10 is OUF isotropic"
## [1] 93
## [1] "The selected model for 9_2 is OUF isotropic"
## [1] 94
## [1] "The selected model for 9_3 is OUF isotropic"
## [1] 95
## [1] "The selected model for 9_4 is OUF isotropic"
## [1] 96
## [1] "The selected model for 9_5 is OUF isotropic"
## [1] 97
## [1] "The selected model for 9_6 is OUF isotropic"
## [1] 98
## [1] "The selected model for 9_7 is OUF isotropic"
## [1] 99
## [1] "The selected model for 9_8 is OUF isotropic"
## [1] 100
## [1] "The selected model for 9_9 is OUF isotropic"Sys.time() - before## Time difference of 1.110633 hourssink()
# 4
ov <- ctmm::overlap(FIT.list, level = 0.95)Obtain median, lower and upper bounds of CI of the overlap and name the matrices with the “BirdID”
The array ov contains three “layers”, with the median of the overlap distribution in the middle one, and the upper and lower bounds of the CI selected in de function in the first and last. We extract the three matrices from there
# separate the layers of the median and CI
X05 <- ov[,,2]
Xlo <- ov[,,1]
Xhi <- ov[,,3]
# remove lower triangle of each (repeated values since it's a symetric matrix)
X05[lower.tri(X05, diag = T)] <- NA
Xlo[lower.tri(Xlo, diag = T)] <- NA
Xhi[lower.tri(Xhi, diag = T)] <- NA
# assign value of BirdID that we rescue from the simulated_dataset df to rows and columns
gid <- simulated_dataset@data[!duplicated(simulated_dataset@data[, "TrackID"]), ][["BirdID"]]
rownames(X05) <- colnames(X05) <- gid
rownames(Xlo) <- colnames(Xlo) <- gid
rownames(Xhi) <- colnames(Xhi) <- gidSeparate within (WI) and between (BW) group overlaps
The function works by testing how many of the Between Individual overlap CIs contain the median of all the within individual overlaps. It provides three values, the proportion of between individual overlap CIs that contain the median of the WI overlaps, the proportion of BW overlap CIs in which the median of the WI overlaps falls outside of the BW and is higher (above) and the proportion in which is lower (below).
# median
WI05 <- NULL
BW05 <- NULL
# BW lower and higher CI bounds
BWlo <- NULL
BWhi <- NULL
for (i in seq_along(rownames(X05))) {
# i = 2
# medians
x1 <- X05[i,]
x2 <- x1[which(names(x1) == rownames(X05)[i])]
x3 <- x1[which(names(x1) != rownames(X05)[i])]
WI05 <- c(WI05, x2)
BW05 <- c(BW05, x3)
# BW intervals
y1 <- Xlo[i,]
y2 <- y1[which(names(y1) != rownames(Xlo)[i])]
BWlo <- c(BWlo, y2)
z1 <- Xhi[i,]
z2 <- z1[which(names(z1) != rownames(Xhi)[i])]
BWhi <- c(BWhi, z2)
}
# remove NAs
WI05 <- WI05[!is.na(WI05)]
BW05 <- BW05[!is.na(BW05)]
BWlo <- BWlo[!is.na(BWlo)]
BWhi <- BWhi[!is.na(BWhi)]
# organize values in a dataframe for plotting
Overlaps <- data.frame(
Overlap = c(WI05, BW05),
Type = c(rep("Within", length(WI05)), rep("Between", length(BW05)))
)Visualising results
First of all we visualise the different densities of the WI and BW medians (without taking into account CI of individual overlaps)
(plot <- ggplot(data = Overlaps) +
geom_density(aes(x = Overlap, col = Type, fill = Type), alpha = 0.5) +
ggtitle("Density of within individual vs. between individual median overlaps") +
theme_bw())
Calculating the WBR (Within-Between ratio)
Now we calculate the median of the WI overlaps, to use as reference values (see if it is contained, above or below the BW CIs)
WI_median <- quantile(WI05, 0.5)
BWint <- data.frame(Low = BWlo,
High = BWhi)Finally we produce a one row dataset containing the proportion of Between individual overlaps in which the Within Individual median is below, inside or above the CI.
BWint <- BWint %>%
mutate(Contains = "NA") %>%
mutate(Contains = if_else(High > WI_median, "NotAbove", "Above (WBR)"))
proportions <- as.data.frame(prop.table(table(BWint$Contains)))
names(proportions) <- c("WI_Median_Position", "Percentage")
proportions$Percentage <- round(proportions$Percentage*100, 2)
(proportions_wide <- spread(proportions, WI_Median_Position, Percentage))## Above (WBR) NotAbove
## 1 83.09 16.91The Within-Between ratio is 76.8% (i.e., in 76.8% of the cases, the median of the Within Individual overlaps falls above the 95% CI of the between individual overlap distribution). This would indicate that, when run on an artificial dataset, purposely simulating a high degree of individual site fidelity, our function is able to detect individual site fidelity in the simulated dataset.
For animals with the same origin but different outward angles
Generating the data
In case of foraging trips of central place foragers, all trips will have the same origin. We create a dataset that simulates those conditions too, to see if the indEffectTest function is still able to detect individual site fidelity.
Again, we simulate the dataset with resistance = 1, which will simulate high ISF conditions
set.seed(2106)
CRW <- species(state.CRW(0.2))
test <- simulate(CRW, coords = c(0, 0), time = 100)
# create empty base raster for resistance surfaces
baseRaster <- raster(xmn = -50, xmx = 50, ymn = -50, ymx = 50, res = 0.1)
values(baseRaster) <- 1
baseGrid <- as(baseRaster, "SpatialPixels")
# create a vector of random angles
angles <- runif(10, 0, 360)
# simulate 10 trips for each origin (i.e. 10 trips per individual)
template_tracks <- list()
template_UDs <- list()
resistance_surfaces <- list()
simulated_tracks <- list()
for (i in 1:10) {
# i = 1
# First we need to create a template trip to get the UD for the resistance surface
template <- simulate(CRW, coords = c(0, 0), angles = angles[i], time = 100)
template_tracks[i] <- list(template)
spTemplate <- SpatialPoints(coords = template[,1:2], proj4string = baseGrid@proj4string)
KDE <- kernelUD(spTemplate, h = "href", grid = baseGrid)
UD50 <- getverticeshr(KDE, percent = 50)
template_UDs[i] <- list(UD50)
# plot(KDE)
# plot(UD50, add = T)
# plot(spTemplate, pch = 20, cex = 0.5, add = T)
# Second, we create the resistance surface
res_surface <- resistanceFromShape(UD50, baseRaster = baseRaster, field = 0,
background = 1)
resistance_surfaces[i] <- list(res_surface)
# plot(res_surface)
# Third, we simulate 10 tracks from each surface (i.e. from each individual)
sim_tracks <- simulate(list(CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW, CRW),
start.resistance = 0,
time = 100, resist = res_surface)
spSim <- SpatialPointsDataFrame(coords = cbind(c(sim_tracks[,1], sim_tracks[,4],
sim_tracks[,7], sim_tracks[,10],
sim_tracks[,13], sim_tracks[,16],
sim_tracks[,19], sim_tracks[,22],
sim_tracks[,25],sim_tracks[,28]),
c(sim_tracks[,2], sim_tracks[,5],
sim_tracks[,8], sim_tracks[,11],
sim_tracks[,14],sim_tracks[,17],
sim_tracks[,20], sim_tracks[,23],
sim_tracks[,26], sim_tracks[,29])),
data = data.frame(TrackID = paste(i, rep(1:10,
each = 100),
sep = "_"),
BirdID = rep(i, 1000)))
spSim@proj4string <- baseGrid@proj4string
simulated_tracks[i] <- list(spSim)
spSim$BirdID <- NULL
KUD_sim <- kernelUD(spSim, grid = baseGrid)
# print(mean(kerneloverlap(spSim, method = "BA", percent = 50, conditional = T)))
# plot(spSim, add = T, col = spSim$TrackID, pch = 20, cex = 0.5)
}We adapt data to indEffectTest same as before
list.df <- lapply(seq_along(simulated_tracks),
function(i){simulated_tracks[i][[1]]})
# apply rbind to the list
simulated_dataset <- do.call("rbind",
list.df)
# assign null projection as X/Y values don't mean anything in this case
simulated_dataset@proj4string <- CRS(as.character(NA))Visualising the data
# first create a version of the data without BirdID to calculate trip overlaps
sds <- as.data.frame(simulated_dataset)
# convert it in order to be able to plot it
sds %>%
st_as_sf(coords = c("coords.x1", "coords.x2")) %>%
tidyr::separate(TrackID, into = c("BirdID", "TrackID"), sep = "_") -> sdf
col2 <- ggplot() +
geom_sf(data = sdf, aes(col = BirdID), alpha = 0.3, size = 0.7) +
xlab("Longitude") + ylab("Latitude") +
# scale_shape_manual(values = 1:10) +
theme_bw() +
theme(legend.position = "right", legend.box = "horizontal") +
guides(col = guide_legend(override.aes = list(size = 3, alpha = 1))) +
NULL
print(col2)
#ggsave(filename = "HighISF_colonial.png", width = 12*1.5, height = 10*1.5, dpi = 800, units = "cm")Testing the function
Calculate akde for all tracks and pairwise overlaps
# To calculate the overlap we:
# 1 - Turn simulated tracks into "telemetry" objects"
# 2 - fit a guesstimate movement model to each track
# 3 - using the parameters of the guesstimate, fit a series of competing movement models
# 4 - Calculate the debiased BA between top-selected models for each track
# (uncertainty from the model fits is propagated into the overlap estimate under the approximation that the BA is a chi-square random variable)
# 1
sds %>%
dplyr::group_by(TrackID) %>%
dplyr::mutate(Date_Time = seq(ymd_hms('2019-11-30 12:10:00'), ymd_hms('2020-01-30 12:10:00'),
length.out = 100)) %>%
mutate(Date_Time = format(Date_Time, format = "%Y/%m/%d %H:%M:%S", tz = "UTC")) %>%
dplyr::select(individual.local.identifier = TrackID,
Longitude = coords.x1, Latitude = coords.x2,
datetime = Date_Time) %>%
ungroup() -> sds
sds <- as.data.frame(sds)
sd_tel <- as.telemetry(object = sds, timeformat = "%Y/%m/%d %H:%M:%S")
FIT.list <- list()
file <- file("FIT.list_Colonial.txt", open = "wt")
sink(file, type = "message") # sink text messages to .txt to avoid them appearing on the markdown
before <- Sys.time()
for(i in 1:length(sd_tel)){
print(i)
one <- sd_tel[[i]]
# 2
GUESS <- ctmm.guess(one, interactive=F)
# 3
FITS <- ctmm.select(one, GUESS, trace=TRUE, verbose=TRUE, cores=-1)
print(paste("The selected model for", names(sd_tel)[i], "is", names(FITS)[1], sep = " "))
FIT.list[[i]] <- FITS[[1]]
}
Sys.time() - before
sink()
# 4
ov <- ctmm::overlap(FIT.list, level = 0.95)Obtain median, lower and upper bounds of CI of the overlap and name the matrices with the “BirdID”
The array ov contains three “layers”, with the median of the overlap distribution in the middle one, and the upper and lower bounds of the CI selected in de function in the first and last. We extract the three matrices from there
# separate the layers of the median and CI
X05 <- ov[,,2]
Xlo <- ov[,,1]
Xhi <- ov[,,3]
# remove lower triangle of each (repeated values since it's a symetric matrix)
X05[lower.tri(X05, diag = T)] <- NA
Xlo[lower.tri(Xlo, diag = T)] <- NA
Xhi[lower.tri(Xhi, diag = T)] <- NA
# assign value of BirdID that we rescue from the simulated_dataset df to rows and columns
gid <- simulated_dataset@data[!duplicated(simulated_dataset@data[, "TrackID"]), ][["BirdID"]]
rownames(X05) <- colnames(X05) <- gid
rownames(Xlo) <- colnames(Xlo) <- gid
rownames(Xhi) <- colnames(Xhi) <- gidSeparate within (WI) and between (BW) group overlaps
The function works by testing how many of the Between Individual overlap CIs contain the median of all the within individual overlaps. It provides three values, the proportion of between individual overlap CIs that contain the median of the WI overlaps, the proportion of BW overlap CIs in which the median of the WI overlaps falls outside of the BW and is higher (above) and the proportion in which is lower (below).
# mediana
WI05 <- NULL
BW05 <- NULL
# BW lower and higher CI bounds
BWlo <- NULL
BWhi <- NULL
for (i in seq_along(rownames(X05))) {
# i = 2
# medians
x1 <- X05[i,]
x2 <- x1[which(names(x1) == rownames(X05)[i])]
x3 <- x1[which(names(x1) != rownames(X05)[i])]
WI05 <- c(WI05, x2)
BW05 <- c(BW05, x3)
# BW intervals
y1 <- Xlo[i,]
y2 <- y1[which(names(y1) != rownames(Xlo)[i])]
BWlo <- c(BWlo, y2)
z1 <- Xhi[i,]
z2 <- z1[which(names(z1) != rownames(Xhi)[i])]
BWhi <- c(BWhi, z2)
}
# remove NAs
WI05 <- WI05[!is.na(WI05)]
BW05 <- BW05[!is.na(BW05)]
BWlo <- BWlo[!is.na(BWlo)]
BWhi <- BWhi[!is.na(BWhi)]
# organize values in a dataframe for plotting
Overlaps <- data.frame(
Overlap = c(WI05, BW05),
Type = c(rep("Within", length(WI05)), rep("Between", length(BW05)))
)Visualising results
First of all we visualise the different densities of the WI and BW medians (without taking into account CI of individual overlaps)
(plot <- ggplot(data = Overlaps) +
geom_density(aes(x = Overlap, col = Type, fill = Type), alpha = 0.5) +
ggtitle("Density of within individual vs. between individual median overlaps") +
theme_bw())
Calculating the WBR (Within-Between ratio)
Now we calculate the median of the WI overlaps, to use as reference values (see if it is contained, above or below the BW CIs)
WI_median <- quantile(WI05, 0.5)
BWint <- data.frame(Low = BWlo,
High = BWhi)Finally we produce a one row dataset containing the proportion of Between individual overlaps in which the Within Individual median is below, inside or above the CI.
BWint <- BWint %>%
mutate(Contains = "NA") %>%
mutate(Contains = if_else(High > WI_median, "NotAbove", "Above (WBR)"))
proportions <- as.data.frame(prop.table(table(BWint$Contains)))
names(proportions) <- c("WI_Median_Position", "Percentage")
proportions$Percentage <- round(proportions$Percentage*100, 2)
(proportions_wide <- spread(proportions, WI_Median_Position, Percentage))## Above (WBR) NotAbove
## 1 62.69 37.31The Within-Between ratio is 62% (i.e., in 62% of the cases, the median of the Within Individual overlaps falls above the 95% CI of the between individual overlap distribution). This would indicate that, even when all trips start in the same origin (central place foragers scenario), our function is able to detect individual site fidelity in the simulated dataset.