Exploring capture - recapture data (MUKA)
Serge MORAND
2024-07-28
1 Managing collected data
Three data types are monitored in the field:
- 1. data sensors: a
georeferenced list if sensors (trapping cages, camera traps,
recorders)
- 2. data animals: a list of individual animals trapped
in each session, with information on taxonomy, morphology, sensor, tag
…
- 3. data capture: a list of captures per date and per sensor
1.0.1 Libraries required
## library(leaflet)
## library(mapview)
## library(leafpop)
## library(knitr)
## library(swimplot)
## library(vegan)
## library(Rcapture)
## library(marked)
## library(ggplot2)
## library(ctmm)
## library(move)
## library(secr)
## library(anipaths)
## library(ggmap)
## library(gganimate)
## library(wildlifeDI)
## library(survival)
## library(coxme)
## library(survminer)
##
## We also called dplyr, tidyr, readxl, sp, sf, DT, ggpubr, kableExtra
1.1 Get locations of rodent traps (data.sensors) (excluding locations of camera traps)
trap_location <- readxl::read_excel("/Users/smorand/Documents/MUKA/data/MUKA_Rodents_Database/MUKA Grid11 Daily Trap Checking.xlsx", sheet="trapID_GPS") |>
dplyr::filter(!trapID %in% c("MUKACam37", "MUKACam39", "MUKACam31", "MUKACam41", "MUKACam28"))
1.1.1 Convert long lat coordinates into UTM coordinates (Zone = 47)
rodent_gps <- trap_location
sp::coordinates(rodent_gps) <- c("long", "lat")
sp::proj4string(rodent_gps) <- CRS("+proj=longlat +datum=WGS84")
res <- sp::spTransform(rodent_gps, CRS("+proj=utm +zone=47 ellps=WGS84"))
as(res, "SpatialPoints") # SpatialPoints object ## class : SpatialPoints
## features : 100
## extent : 515939.5, 516035.6, 1561551, 1561648 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=47 +ellps=WGS84 +units=m +no_defslatlon <- data.frame(as(res, "SpatialPoints"))
names(latlon)<- c("lon_utm","lat_utm")
latlon <- latlon |>
dplyr::mutate(lon_utm = as.integer(round(lon_utm,0))) |>
dplyr::mutate(lat_utm = as.integer(round(lat_utm,0)))
trap_location <- cbind(trap_location, latlon) |>
dplyr::select(trapID,long, lat, lon_utm, lat_utm)
1.1.2 Table of trap locations
DT::datatable(trap_location, class = 'cell-border stripe',
# colnames = c('trapID', 'Zone', 'Lat', 'Lon'),
caption = 'Table 1: "List of trap locations"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))
1.1.3 Prepare map using leaflet
opts_chunk$set(fig.width=10, fig.height=6)
map_trap <- leaflet() %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addMiniMap(tiles = providers$Esri.WorldStreetMap,position ="bottomright",
toggleDisplay = TRUE)%>%
addCircleMarkers(data=trap_location,
~long, ~lat,
label = ~as.character(trapID),
labelOptions = labelOptions(noHide = FALSE, textOnly = FALSE,
style=list(
'color'='red',
'font-family'= 'serif',
'font-style'= 'italic',
'box-shadow' = '3px 3px rgba(0,0,0,0.25)',
'font-size' = '12px',
'border-color' = 'rgba(0,0,0,0.5)'
)),
fillColor = "red",
stroke = FALSE, fillOpacity = 0.5,
clusterOptions = markerClusterOptions(removeOutsideVisibleBounds = F, maxClusterRadius = 10))%>%
addLayersControl(baseGroups = c("Google","WSM","OSM"),
options = layersControlOptions(collapsed = FALSE))
1.2 Get rodent data (data.animals)
rodent_data <- readxl::read_excel("/Users/smorand/Documents/MUKA/data/MUKA_Rodents_Database/MUKA Grid11 Daily Trap Checking.xlsx",
sheet = "Rodent_MUKA_database") |>
dplyr::select(id_indiv,trapID,taxonomyID_field,trappingDate,sex,maturity) |>
dplyr::mutate(taxonomyID_field = stringr::str_to_title(taxonomyID_field)) |>
dplyr::mutate(taxonomyID_field = gsub("_"," ",taxonomyID_field)) |>
dplyr::arrange(id_indiv)1.2.1 Table of rodent data
DT::datatable(rodent_data, class = 'cell-border stripe',
caption = 'Table 2: "List of rodent individuals"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))
1.2.2 Map of first invidual rodent captures (data.animals)
# join rodent_data with trap locations
rodent_location <- dplyr::left_join(rodent_data,trap_location,by="trapID") |>
tidyr::drop_na()
opts_chunk$set(fig.width=10, fig.height=6)
map_rodent <- leaflet() %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addMiniMap(tiles = providers$Esri.WorldStreetMap,position ="bottomright",
toggleDisplay = TRUE)%>%
addCircleMarkers(data=rodent_location,
~long, ~lat,
label = ~as.character(taxonomyID_field),
labelOptions = labelOptions(noHide = FALSE, textOnly = FALSE,
style=list(
'color'='red',
'font-family'= 'serif',
'font-style'= 'italic',
'box-shadow' = '3px 3px rgba(0,0,0,0.25)',
'font-size' = '12px',
'border-color' = 'rgba(0,0,0,0.5)'
)),
fillColor = "red",
group = "rodents",
stroke = FALSE, fillOpacity = 0.5,
clusterOptions = markerClusterOptions(removeOutsideVisibleBounds = F, maxClusterRadius = 10))%>%
addCircleMarkers(data=trap_location,
~long, ~lat,
label = ~as.character(trapID),
labelOptions = labelOptions(noHide = FALSE, textOnly = FALSE,
style=list(
'color'='yellow',
'font-family'= 'serif',
'font-style'= 'italic',
'box-shadow' = '3px 3px rgba(0,0,0,0.25)',
'font-size' = '12px',
'border-color' = 'rgba(0,0,0,0.5)'
)),
fillColor = "yellow",
group = "traps",
stroke = FALSE, fillOpacity = 0.5,
clusterOptions = markerClusterOptions(removeOutsideVisibleBounds = F, maxClusterRadius = 10))%>%
addLayersControl(baseGroups = c("Google","WSM","OSM"),
overlayGroups = c("rodents", "traps"),
options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup("traps")
1.4 Get the capture and recapture data (data.EpicCollect)
rodent_capture <- readxl::read_excel("/Users/smorand/Documents/MUKA/data/MUKA_Rodents_Database/MUKA Grid11 Daily Trap Checking.xlsx",
sheet = "trap_followup") |>
dplyr::select(trapID,gridID,session, trappingDate,id_indiv,release) |>
dplyr::filter(id_indiv !="na") |>
dplyr::mutate(trappingDate= as.POSIXlt(trappingDate, format="%Y-%m-%d")) |>
dplyr::mutate(trappingDateH = trappingDate + lubridate::hours(9) + lubridate::minutes(0) + lubridate::seconds(0)) |>
dplyr::mutate(occasion = 1) |>
dplyr::left_join(rodent_data |>
dplyr::select(id_indiv,taxonomyID_field), by=c("id_indiv")) |>
dplyr::left_join(trap_location |>
dplyr::select(trapID,long,lat), by = "trapID")
1.5 Save the 3 data sets (sensors, animals, followup) in a RData file
1.5.1 What is a RData file?
A RData format (with extension .RData) is a format designed for use with R for storing selected “objects” from a workspace in a form that can be loaded back by R
2 Visualizing individual captures
2.1 Preparing data
start_session <- min(rodent_capture$trappingDate)
individuals_capt <- rodent_capture |>
dplyr::group_by(id_indiv) |>
dplyr::summarise(start = as.numeric((min(trappingDate) - start_session), units = "days"),
end = as.numeric((max(trappingDate) - start_session), units = "days"),
start_date = min(trappingDate),
end_date = max(trappingDate),
n_days = end-start) |>
#dplyr::arrange(start_date) |>
dplyr::left_join(rodent_data |>
dplyr::select(id_indiv, taxonomyID_field, sex, maturity), by = "id_indiv") |>
dplyr::left_join(rodent_capture |>
dplyr::select(id_indiv, release) |>
dplyr::filter(release == "no"), by = "id_indiv")
2.1.1 Table of individual captures
DT::datatable(individuals_capt, class = 'cell-border stripe',
caption = 'Table 3: "List of rodent individuals"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))
2.1.2 Individual captures through time
individuals_capt <- as.data.frame(individuals_capt)
ggplot() +
swimmer_lines(df=individuals_capt,id='id_indiv',start="start", end='end',
name_col='maturity',size=4) +
swimmer_arrows(df_arrows = individuals_capt |>
dplyr::mutate(release = ifelse(is.na(release), "yes", release)) |>
dplyr::mutate(survival = dplyr::recode(release,
'yes' = 'alive',
'no' = 'dead')) |>
dplyr::filter(survival == "alive"),
id='id_indiv',arrow_start='end',
cont = 'survival',name_col='survival', type = "open",cex=1, angle = 30) +
swimmer_points(df = individuals_capt |>
dplyr::mutate(release = ifelse(is.na(release), "yes", release)) |>
dplyr::mutate(survival = dplyr::recode(release,
'yes' = 'alive',
'no' = 'dead')) |>
dplyr::filter(survival == "dead"),
time='end', name_col = "survival",
id='id_indiv',name_shape = 'survival',size=5) +
ggplot2::theme_bw() +
scale_y_continuous(name = paste("Time since starting session 5 (",start_session,") in days",
sep = " "),breaks = seq(0,35,by=5)) +
scale_x_discrete(name = "Individuals") +
coord_flip()
2.1.3 Number of individual captures according to sex
No differences in number of individual captures according to sex
ggpubr::ggviolin(individuals_capt |>
dplyr::select(n_days,sex) |>
dplyr::filter(sex != "na"),
x ="sex" , y = "n_days",
fill = "sex",
add = "boxplot", add.params = list(fill = "white")) +
ggpubr::stat_compare_means(comparisons = c("M","F"),
label = "p.signif")+ # Add significance levels
ggpubr::stat_compare_means(label.y = 30) # Add global the p-value 
2.1.4 Number of individual captures according to maturity
No differences in the number of individual captures according to maturity
ggpubr::ggviolin(individuals_capt,
x ="maturity" , y = "n_days",
fill = "maturity",
add = "boxplot", add.params = list(fill = "white")) +
ggpubr::stat_compare_means(comparisons = c("adult","juvenile"),
label = "p.signif")+ # Add significance levels
ggpubr::stat_compare_means(label.y = 30) # Add global the p-value 
3 Small mammal species richness
3.0.1 Estimation of species richness using the package “vegan”
## @Manual{,
## title = {vegan: Community Ecology Package},
## author = {Jari Oksanen and Gavin L. Simpson and F. Guillaume Blanchet and Roeland Kindt and Pierre Legendre and Peter R. Minchin and R.B. O'Hara and Peter Solymos and M. Henry H. Stevens and Eduard Szoecs and Helene Wagner and Matt Barbour and Michael Bedward and Ben Bolker and Daniel Borcard and Gustavo Carvalho and Michael Chirico and Miquel {De Caceres} and Sebastien Durand and Heloisa Beatriz Antoniazi Evangelista and Rich FitzJohn and Michael Friendly and Brendan Furneaux and Geoffrey Hannigan and Mark O. Hill and Leo Lahti and Dan McGlinn and Marie-Helene Ouellette and Eduardo {Ribeiro Cunha} and Tyler Smith and Adrian Stier and Cajo J.F. {Ter Braak} and James Weedon},
## year = {2024},
## note = {R package version 2.6-6.1},
## url = {https://CRAN.R-project.org/package=vegan},
## }
3.1 Estimation of species richness (all small mammals)
3.1.1 Preparing a table of number of captured animals by rodent species and by capture date
species_observed <- table(rodent_data$trappingDate,rodent_data$taxonomyID_field)
kableExtra::kable(species_observed)| Berylmys berdmorei | Leopoldamys neilli | Maxomys surifer | Rattus tanezumi | Tupaia belangeri | |
|---|---|---|---|---|---|
| 2024-06-26 | 0 | 0 | 2 | 0 | 0 |
| 2024-06-27 | 0 | 0 | 1 | 1 | 0 |
| 2024-06-28 | 0 | 0 | 1 | 1 | 0 |
| 2024-06-29 | 0 | 0 | 1 | 0 | 1 |
| 2024-06-30 | 0 | 0 | 2 | 0 | 0 |
| 2024-07-01 | 0 | 0 | 1 | 0 | 0 |
| 2024-07-02 | 0 | 0 | 4 | 0 | 0 |
| 2024-07-03 | 0 | 0 | 1 | 0 | 0 |
| 2024-07-04 | 1 | 0 | 4 | 0 | 0 |
| 2024-07-05 | 0 | 0 | 7 | 0 | 0 |
| 2024-07-06 | 0 | 0 | 0 | 2 | 0 |
| 2024-07-12 | 0 | 1 | 2 | 1 | 0 |
| 2024-07-17 | 0 | 0 | 1 | 0 | 1 |
| 2024-07-19 | 0 | 0 | 3 | 0 | 0 |
3.3 Species accumulation curve using ggplot2 and BiodiversityR
Accumulation <- BiodiversityR::accumresult(species_observed, method = "exact")
accum <- data.frame(Accumulation$sites,Accumulation$richness, Accumulation$sd)
names(accum) <- c("sites","richness","sd")
plotgg <- ggplot(accum, aes(x = sites, y = richness)) +
scale_x_continuous(expand=c(0, 1), sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_line(aes(colour="red"), linewidth=2) +
geom_point(data=subset(accum),
aes(colour="red"), size=5) +
geom_errorbar(aes(ymin=richness-sd, ymax=richness+sd,colour="red"), width=.2,
position=position_dodge(.9)) +
geom_ribbon(aes(ymin=richness-sd, ymax=richness+sd, x= sites,fill = "band", color ="pink"), alpha = 0.3)+
theme_classic() +
xlab("Day (sampling)") +
ylab("Rodent Species Richness") +
theme(legend.position = "none")
plotgg
3.4 Estimators of species richness: Observedn Chao, Jackknife 1, Jackknife 2, bootstrap.
The estimators gave 6 (bootstrap) to 7 species (Chao, Jackknife 1) to be present; Jackknife 2 usually sur-estimates the number of species
3.4.1 A list of diversity indices
## 1 "w" = (b+c)/(2*a+b+c)
## 2 "-1" = (b+c)/(2*a+b+c)
## 3 "c" = (b+c)/2
## 4 "wb" = b+c
## 5 "r" = 2*b*c/((a+b+c)^2-2*b*c)
## 6 "I" = log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) +
## (a+c)*log(a+c)) / (2*a+b+c)
## 7 "e" = exp(log(2*a+b+c) - 2*a*log(2)/(2*a+b+c) - ((a+b)*log(a+b) +
## (a+c)*log(a+c)) / (2*a+b+c))-1
## 8 "t" = (b+c)/(2*a+b+c)
## 9 "me" = (b+c)/(2*a+b+c)
## 10 "j" = a/(a+b+c)
## 11 "sor" = 2*a/(2*a+b+c)
## 12 "m" = (2*a+b+c)*(b+c)/(a+b+c)
## 13 "-2" = pmin.int(b,c)/(pmax.int(b,c)+a)
## 14 "co" = (a*c+a*b+2*b*c)/(2*(a+b)*(a+c))
## 15 "cc" = (b+c)/(a+b+c)
## 16 "g" = (b+c)/(a+b+c)
## 17 "-3" = pmin.int(b,c)/(a+b+c)
## 18 "l" = (b+c)/2
## 19 "19" = 2*(b*c+1)/(a+b+c)/(a+b+c-1)
## 20 "hk" = (b+c)/(2*a+b+c)
## 21 "rlb" = a/(a+c)
## 22 "sim" = pmin.int(b,c)/(pmin.int(b,c)+a)
## 23 "gl" = 2*abs(b-c)/(2*a+b+c)
## 24 "z" = (log(2)-log(2*a+b+c)+log(a+b+c))/log(2)
3.4.2 Diversity indices for MUKA first trapping session (Shannon and Simpson)
# Shannon
H <- diversity(species_observed)
# Simpson
simp <- diversity(species_observed, "simpson")
# Species richness (S)
S <- specnumber(species_observed)
results <- data.frame(cbind(S, H, simp))
DT::datatable(results, class = 'cell-border stripe',
colnames = c("Species richness",'Shannon (H)', 'Simpson'),
caption = 'Table 4: "Diversity indices by trapping date"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))
4 Abundance, density
4.0.1 Estimation of abundance using “Rcapture” and “marked”
## @Manual{,
## title = {Rcapture: Loglinear Models for Capture-Recapture Experiments},
## author = {Louis-Paul Rivest and Sophie Baillargeon},
## year = {2022},
## note = {R package version 1.4-4},
## url = {https://CRAN.R-project.org/package=Rcapture},
## }## @Article{,
## author = {J.L. Laake and D.S. Johnson and P.B. Conn},
## year = {2013},
## title = {marked: An R package for maximum-likelihood and MCMC analysis of capture-recapture data.},
## journal = {Methods in Ecology and Evolution},
## }
4.1 Abundance estimators and heterogeinity graph (select a species present in the grid; here Maxomys surifer)
The plot.descriptive function produces an exploratory heterogeneity graph. In the absence of heterogeneity, the relation(s) presented in the graph should be almost linear
species_selected <- "Maxomys surifer"
species_selected_recapture <- rodent_capture |>
dplyr::filter(taxonomyID_field == species_selected) |>
dplyr::filter(trappingDate > "2024-06-28" & trappingDate < "2024-07-10") |>
dplyr::select(id_indiv,trappingDate,occasion) |>
dplyr::distinct() |>
tidyr::spread(trappingDate,occasion, fill=0) |>
tibble::column_to_rownames(var = "id_indiv")
desc <- descriptive(species_selected_recapture)
plot(desc)
Loglinear Models for Closed Population Capture-Recapture
Experiments
closedp.0 works with fairly large data sets
##
## Number of captured units: 23
##
## Abundance estimations and model fits:
## abundance stderr deviance df AIC BIC infoFit
## M0 23.8 1.0 24.946 8 50.147 52.418 OK
## Mh Chao (LB) 29.9 5.1 1.772 6 30.973 35.515 OK
## Mh Poisson2 24.1 1.2 19.479 7 46.680 50.086 OK
## Mh Darroch 29.9 4.8 6.592 7 33.793 37.199 OK
## Mh Gamma3.5 50.4 22.0 4.680 7 31.881 35.288 OK4.2 Estimation of survival using marked
The Phi.(Intercept) term in a capture-recapture model represents the logit of the survival probability (Φ)
species_selected_recapture <- rodent_capture |>
dplyr::filter(taxonomyID_field == species_selected) |>
dplyr::filter(trappingDate > "2024-06-28" & trappingDate < "2024-07-10") |>
dplyr::select(id_indiv,trappingDate,occasion) |>
dplyr::distinct() |>
tidyr::spread(trappingDate,occasion, fill=0) |>
tibble::column_to_rownames(var = "id_indiv")
nc <- as.integer(ncol(species_selected_recapture))
species_selected_recapture_s <- species_selected_recapture |>
tidyr::unite(ch,1:nc,sep = "")
model= marked::crm(species_selected_recapture_s)## Number of evaluations: 100 -2lnl: 130.4849949##
## crm Model Summary
##
## Npar : 2
## -2lnL: 130.485
## AIC : 134.485
##
## Beta
## Estimate
## Phi.(Intercept) 2.0751606
## p.(Intercept) 0.4434585##
## crm Model Summary
##
## Npar : 2
## -2lnL: 130.485
## AIC : 134.485
##
## Beta
## Estimate se lcl ucl
## Phi.(Intercept) 2.0751606 0.5062353 1.0829394 3.067382
## p.(Intercept) 0.4434585 0.3121270 -0.1683104 1.055228
5 Autocorrelated Kernel Density Estimation
5.0.1 using ctmm
## @Manual{,
## title = {ctmm: Continuous-Time Movement Modeling},
## author = {Christen H. Fleming and Justin M. Calabrese},
## year = {2023},
## note = {R package version 1.2.0},
## url = {https://CRAN.R-project.org/package=ctmm},
## }
5.0.2 show individuals with more than n (> 4) captures
capture_n <- rodent_capture |>
dplyr::select(id_indiv,taxonomyID_field) |>
dplyr::group_by(id_indiv) |>
dplyr::summarise(n = dplyr::n()) |>
dplyr::filter(n > 4) |>
dplyr::arrange(desc(n))
DT::datatable(capture_n, class = 'cell-border stripe',
caption = 'Table 5: "Individuals with more than 4 captures"',
options = list(columnDefs =
list(list(className = 'dt-center',
targets = "_all"))))
5.0.3 Select an individual. Here RM0249 and convert to a telemetry object
id_selected <- "RM0268"
rodent_individual_selected <- rodent_capture |>
dplyr::filter(id_indiv == id_selected) |>
dplyr::select(trappingDateH,lat, long,id_indiv) |>
dplyr::rename(timestamp = trappingDateH)
# Convert the data frame to a telemetry object
rodent <- as.telemetry(rodent_individual_selected)
rodent
5.1 Estimation of individual density
M.IID <- ctmm.fit(rodent) # no autocorrelation timescales
GUESS <- ctmm.guess(rodent,interactive=FALSE) # automated model guess
M.OUF <- ctmm.fit(rodent,GUESS) # in general, use ctmm.select instead
KDE <- akde(rodent,M.IID) # KDE
AKDE <- akde(rodent,M.OUF) # AKDE
wAKDE <- akde(rodent,M.OUF,weights=TRUE) # weighted AKDE
5.1.1 Plot individual kernel density
# calculate one extent for all UDs
EXT <- extent(list(KDE,AKDE,wAKDE),level=0.95)
plot(rodent,UD=KDE,xlim=EXT$x,ylim=EXT$y)
exp_title <- paste(expression("IID KDE species:"),species_selected,sep = " ",";", "individual: ",id_selected)
title(exp_title)
plot(rodent,UD=AKDE,xlim=EXT$x,ylim=EXT$y)
exp_title <- paste(expression("OUF AKDE species:"),species_selected,sep = " ",";", "individual: ",id_selected)
title(exp_title)
plot(rodent,UD=wAKDE,xlim=EXT$x,ylim=EXT$y)
exp_title <- paste(expression("weighted OUF AKDE species:"),species_selected,sep = " ",";", "individual: ",id_selected)
title(exp_title)
5.2 Density using secr
## @Manual{,
## title = {secr: Spatially explicit capture-recapture models},
## author = {Murray Efford},
## year = {2024},
## note = {R package version 4.6.7},
## url = {https://CRAN.R-project.org/package=secr},
## }## https://globalsnowleopard.org/wp-content/uploads/2020/09/Data-analysis-with-SERC.pdf
5.2.2 Preparing data for secr
load("data_grid11.RData")
## select session
session_selected <- "5"
## select grid
grid_selected <- "G11"
## select species
species_selected <- "Maxomys surifer"
# rodent_selected <- "Rattus tanezumi"
## session starting day
date_initial <- rodent_capture$trappingDate[1]
data_cap <- rodent_capture |>
dplyr::filter(gridID == grid_selected) |>
dplyr::filter(session == session_selected)|>
dplyr::select(trapID,trappingDate,id_indiv) |>
dplyr::mutate(Session = paste("MUKA July 2025,
session = ", session_selected,",
grid = ", grid_selected,",
species = ", species_selected,
collapse=" ")) |>
dplyr::mutate(occasion = as.numeric(difftime(trappingDate,date_initial, units = "days")) +1) |>
dplyr::left_join(rodent_data |>
dplyr::select(-trapID, -trappingDate), by="id_indiv") |>
dplyr::filter(taxonomyID_field == species_selected)|>
dplyr::select(Session,id_indiv,occasion,trapID) |>
dplyr::left_join(trap_location |>
dplyr::select(trapID,lon_utm,lat_utm), by="trapID") |>
dplyr::rename(x = lon_utm, y = lat_utm) |>
dplyr::select(-trapID)
data_trap <- trap_location |>
dplyr::select(trapID,lon_utm, lat_utm) |>
dplyr::rename(x=lon_utm, y=lat_utm)
datatraps <- read.traps(data = data_trap)
data_cap <- data.frame(data_cap)
kableExtra::kable(head(data_cap))| Session | id_indiv | occasion | x | y |
|---|---|---|---|---|
| MUKA July 2025, |
session = 5 ,
grid = G11 ,
species = Maxomys surifer |RM0249 | 1| 515974| 1561615|
5.2.3 Compute
## Object class capthist
## Detector type multi
## Detector number 100
## Average spacing 8.062258 m
## x-range 515939 516036 m
## y-range 1561551 1561648 m
##
## Counts by occasion
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## n 2 1 1 1 5 5 8 8 10 13 12 0 0 6 9
## u 2 1 1 1 2 1 4 1 4 7 0 0 0 0 0
## f 8 5 3 2 1 0 0 2 2 4 0 2 0 0 1
## M(t+1) 2 3 4 5 7 8 12 13 17 24 24 24 24 24 24
## losses 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## detections 2 1 1 1 5 5 8 8 10 13 12 0 0 6 9
## detectors visited 2 1 1 1 5 5 8 8 10 13 12 0 0 6 9
## detectors used 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## n 9 9 7 0 0 5 6 7 8 5 0 0 0 7 5
## u 0 2 0 0 0 0 1 0 3 0 0 0 0 0 0
## f 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## M(t+1) 24 26 26 26 26 26 27 27 30 30 30 30 30 30 30
## losses 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## detections 9 9 7 0 0 5 6 7 8 5 0 0 0 7 5
## detectors visited 9 9 7 0 0 5 6 7 8 5 0 0 0 7 5
## detectors used 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 31 Total
## n 4 153
## u 0 30
## f 0 30
## M(t+1) 30 30
## losses 0 0
## detections 4 153
## detectors visited 4 153
## detectors used 100 3100## $`M(t+1)`
## 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
## M(t+1) 2 3 4 5 7 8 12 13 17 24 24 24 24 24 24 24 26 26 26 26 26 27 27 30 30 30
## 27 28 29 30 31 Total
## M(t+1) 30 30 30 30 30 30
5.3 Map
5.4 Estimation of density
initialsigma <- RPSV(CHxy, CC = FALSE)
fit <- secr.fit(CHxy, buffer = 4 * initialsigma, trace = FALSE)
fit <- secr.fit(CHxy, buffer= 100, trace = FALSE)
cat("Quick and biased estimate of sigma =", initialsigma, "m\n")## Quick and biased estimate of sigma = 24.53116 m##
## secr.fit(capthist = CHxy, buffer = 100, trace = FALSE)
## secr 4.6.7, 19:17:02 28 Jul 2024
##
## Detector type multi
## Detector number 100
## Average spacing 8.062258 m
## x-range 515939 516036 m
## y-range 1561551 1561648 m
##
## N animals : 30
## N detections : 153
## N occasions : 31
## Mask area : 7.496473 ha
##
## Model : D~1 g0~1 sigma~1
## Fixed (real) : none
## Detection fn : halfnormal
## Distribution : poisson
## N parameters : 3
## Log likelihood : -912.4522
## AIC : 1830.904
## AICc : 1831.827
##
## Beta parameters (coefficients)
## beta SE.beta lcl ucl
## D 2.178647 0.2190901 1.749238 2.608055
## g0 -3.820169 0.1566670 -4.127230 -3.513107
## sigma 3.238506 0.0879749 3.066078 3.410934
##
## Variance-covariance matrix of beta parameters
## D g0 sigma
## D 0.048000484 -0.007523862 -0.010263803
## g0 -0.007523862 0.024544540 0.001856563
## sigma -0.010263803 0.001856563 0.007739584
##
## Fitted (real) parameters evaluated at base levels of covariates
## link estimate SE.estimate lcl ucl
## D log 8.83434245 1.958977596 5.75021887 13.57263231
## g0 logit 0.02145375 0.003288986 0.01587152 0.02894159
## sigma log 25.49560510 2.247320335 21.45759034 30.29351707
5.6 Estimation of density (using several models)
# Estimation
secr.fit(CHxy, buffer = 100)
# density
fit0 <- list(CHxy, model = D ~1)
fitg0 <- list(CHxy, model = g0~b)
fitD <- list(CHxy, model = D ~ x + y)
fitD6 <- list(CHxy, model= D ~ s(x, y))
fitDxy2 <- list(CHxy, model = D ~ x + y + x2 + y2 + xy)
fits <- par.secr.fit (c('fit0','fitD','fitg0','fitDxy2'))
#
AIC(fits)
region.N(fits$fit.fitg0, alpha = 0.5,spacing = 40,
keep.region = TRUE)
temp <- region.N(fits$fit.fitg0,alpha = 0.5,spacing = 40,
keep.region = TRUE)
par.derived(fits, se.esa = FALSE)
par.region.N(fits)
#
surfaceDxy2 <- predictDsurface(fits$fit.fitDxy2)
5.7 Plot estimation density
trap_sf = sf::st_as_sf(trap_location, coords = c("lon_utm", "lat_utm"), crs = 24047, agr = "constant")
### plot
plot(surfaceDxy2, plottype = "shaded", poly = FALSE, breaks = seq(0,20,0.05),
title = "Density / ha", text.cex = 1)
plot(surfaceDxy2, plottype = "contour", poly = FALSE, breaks = seq(0,10,1),
add = TRUE)
plot(trap_sf, add=TRUE,pch = 1, col="red")
plot(CHxy,tracks = TRUE, cex=10,add=TRUE)
6 Movement
6.0.1 compute distance of recapture to initial capture (using geosphere)
rodent_capture <- rodent_capture |>
dplyr::mutate(distance = NA) # add a column for distance value
# a loop to compute the distance between each recapture to the initial capture for each individual
for (i in 1:nrow(rodent_capture)){
print(i)
log_i <- rodent_capture[i,] |>
dplyr::select(long,lat)
log_i[1,]
loc_old_i <-rodent_capture |>
dplyr::filter(id_indiv ==rodent_capture$id_indiv[i]) |>
dplyr::filter(trappingDate == min(trappingDate)) |>
dplyr::select(long,lat)
loc_old_i[1,]
dist_i <- geosphere::distGeo(log_i[1,],loc_old_i[1,])
print(dist_i)
rodent_capture$distance[i] <- as.integer(dist_i)
}
6.0.2 show table distance
kableExtra::kable(head(rodent_capture |>
dplyr::select(id_indiv,trapID, taxonomyID_field, distance)))| id_indiv | trapID | taxonomyID_field | distance |
|---|---|---|---|
| RM0249 | G11064 | Maxomys surifer | 0 |
| RM0248 | G11095 | Maxomys surifer | 0 |
| RM0253 | G11033 | Maxomys surifer | 0 |
| RM0254 | G11045 | Rattus tanezumi | 0 |
| RM0254 | G11033 | Rattus tanezumi | 19 |
| RM0257 | G11034 | Maxomys surifer | 0 |
6.1 plot travelled distances for all rodents
ggplot(rodent_capture,
aes(distance)) +
geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
theme_minimal()+
scale_fill_brewer() +
labs(title= "All individuals")+
theme(plot.title = element_text(size = 16, face = "bold.italic"),
axis.title=element_text(size=14,face="bold")) 
6.2 plot travelled distances for Maxomys surifer
species_selected <- "Maxomys surifer"
ggplot(rodent_capture_select_species <- rodent_capture |>
dplyr::filter(taxonomyID_field == species_selected),
aes(distance)) +
geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
theme_minimal()+
scale_fill_brewer() +
labs(title= paste0(species_selected))+
theme(plot.title = element_text(size = 16, face = "bold.italic"),
axis.title=element_text(size=14,face="bold")) 
6.3 plot distance for Rattus tanezumi
species_selected <- "Rattus tanezumi"
ggplot(rodent_capture_select_species <- rodent_capture |>
dplyr::filter(taxonomyID_field == species_selected),
aes(distance)) +
geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
theme_minimal()+
scale_fill_brewer() +
labs(title= paste0(species_selected))+
theme(plot.title = element_text(size = 16, face = "bold.italic"),
axis.title=element_text(size=14,face="bold")) 
6.4 individual movement (individual RM0268)
##
data_time <- rodent_capture |>
dplyr::mutate(time = paste(trappingDate,"01:01:01", sep = " ")) |>
dplyr::mutate(time = as.POSIXct(time,format="%Y-%m-%d%H:%M")) |>
dplyr::select(time,id_indiv,long,lat)
individual_selected <- "RM0268"
data_time_indiv <- data_time |>
dplyr::filter(id_indiv == individual_selected)|>
dplyr::select(long,lat,time)|>
dplyr::arrange(sort(time))|>
data.matrix()
line_rod <- sf::st_sf(
data.frame(id = individual_selected),
geometry = sf::st_sfc(sf::st_linestring(data_time_indiv)),
crs = 4326) |>
sf::st_zm(drop = T, what = "ZM")
individualMap <- mapview::mapview(line_rod)
grid11_sf <- sf::st_as_sf(trap_location, coords = c("long","lat"),
crs = 4326)
gridMap <- mapview::mapview(grid11_sf)
gridMap + individualMap
6.4.1 individual movement (individual RM0258)
##
data_time <- rodent_capture |>
dplyr::mutate(time = paste(trappingDate,"01:01:01", sep = " ")) |>
dplyr::mutate(time = as.POSIXct(time,format="%Y-%m-%d%H:%M")) |>
dplyr::select(time,id_indiv,long,lat)
individual_selected <- "RM0258"
data_time_indiv <- data_time |>
dplyr::filter(id_indiv == individual_selected)|>
dplyr::select(long,lat,time)|>
dplyr::arrange(sort(time))|>
data.matrix()
line_rod <- sf::st_sf(
data.frame(id = individual_selected),
geometry = sf::st_sfc(sf::st_linestring(data_time_indiv)),
crs = 4326) |>
sf::st_zm(drop = T, what = "ZM")
individualMap2 <- mapview::mapview(line_rod)
gridMap + individualMap2
6.4.2 individual movements of three individuals (“RM0249”,“RM0253”,“RM0258”, “RM0268”, “RM0257”,“RM0286”)
individual_selected <- c("RM0249","RM0253","RM0258", "RM0268", "RM0257","RM0286")
rodentmove <- rodent_capture |>
dplyr::filter(id_indiv %in% individual_selected)|>
dplyr::select(long, lat, trappingDateH,id_indiv) |>
as.data.frame()
p = ggplot()+
geom_path(data = data_time|>
dplyr::filter(id_indiv %in% individual_selected),
aes(x=long,y=lat,group=id_indiv,colour = id_indiv),
alpha = 0.3)+
geom_point(data = rodentmove,
aes(x=long ,y=lat,group=id_indiv,color = id_indiv),
alpha = 0.7, shape=21, size = 2) +
scale_size_continuous(range = c(0.1,10))+
theme_bw() +
theme(panel.grid = element_blank())
p
7 Interactions among individuals using move2 and wildlifeDI
## @Article{,
## title = {A critical examination of indices of dynamic interaction for wildlife telemetry data},
## author = {Jed Long and Trisalyn Nelson and Stephen Webb and Kenneth Gee},
## journal = {Journal of Animal Ecology},
## year = {2014},
## volume = {83},
## pages = {1216-1233},
## }
##
## @Article{,
## title = {Analyzing contacts and behavior from high frequency tracking data using the wildlifeDI R package},
## author = {Jed Long and Stephen Webb and Seth Harju and Kenneth Gee},
## journal = {Geographical Analysis},
## year = {2022},
## volume = {54},
## pages = {648-663},
## }## @Manual{,
## title = {move2: Processing and Analysing Animal Trajectories},
## author = {Bart Kranstauber and Kamran Safi and Anne K. Scharf},
## year = {2024},
## note = {R package version 0.3.0},
## url = {https://CRAN.R-project.org/package=move2},
## }
7.1 preparing a move2 file for two individuals of Maxomys surifer: RM0253 and RM0257
individuals_selected <- c("RM0249","RM0253","RM0258", "RM0268", "RM0257","RM0286")
rodent_capture_selected <- rodent_capture |>
dplyr::filter(id_indiv %in% individuals_selected) |>
dplyr::mutate(date = as.POSIXct(trappingDateH))|>
dplyr::select(date,id_indiv,trapID) |>
dplyr::left_join(trap_location |>
dplyr::select(trapID, lon_utm,lat_utm), by="trapID") |>
dplyr::select(- trapID) |>
move2::mt_as_move2(coords = c("lon_utm", "lat_utm"), time_column = "date",
track_id_column = "id_indiv") |>
sf::st_set_crs(24047)
checkTO(rodent_capture_selected)
7.2 plot the interaction overlap
idcol <- move2::mt_track_id_column(rodent_capture_selected)
mcphr <- rodent_capture_selected |>
dplyr::group_by_at(idcol) |>
dplyr::summarise() |>
sf::st_convex_hull()
plot(sf::st_geometry(mcphr),border=c("red","black"))
7.2.1 compute interaction indices
#Compute SI index
AIB <- mcphr |>
sf::st_intersection() |>
dplyr::filter(n.overlaps == 2)
sf::st_area(AIB)/sf::st_area(sf::st_union(mcphr))## Units: [1]
## [1] 0.0023477169 0.0487021769 0.0644816383 0.0378639808 0.0207455129
## [6] 0.0491196694 0.0024197958 0.0006530858 0.0216244255 0.0176755240
## [11] 0.0102437941 0.0000000000
8 Survival using survival,coxme and survminer
## @Manual{survival-package,
## title = {A Package for Survival Analysis in R},
## author = {Terry M Therneau},
## year = {2024},
## note = {R package version 3.5-8},
## url = {https://CRAN.R-project.org/package=survival},
## }
##
## @Book{survival-book,
## title = {Modeling Survival Data: Extending the {C}ox Model},
## author = {{Terry M. Therneau} and {Patricia M. Grambsch}},
## year = {2000},
## publisher = {Springer},
## address = {New York},
## isbn = {0-387-98784-3},
## }## @Manual{,
## title = {coxme: Mixed Effects Cox Models},
## author = {Terry M. Therneau},
## year = {2024},
## note = {R package version 2.2-20},
## url = {https://CRAN.R-project.org/package=coxme},
## }## @Manual{,
## title = {survminer: Drawing Survival Curves using 'ggplot2'},
## author = {Alboukadel Kassambara and Marcin Kosinski and Przemyslaw Biecek},
## year = {2021},
## note = {R package version 0.4.9},
## url = {https://CRAN.R-project.org/package=survminer},
## }
8.0.1 preparing data for Maxomys surifer (including sex and maturity)
species_selected <- "Maxomys surifer"
rodent_capture_date <- rodent_capture |>
dplyr::select(id_indiv,trappingDate) |>
dplyr::arrange(id_indiv,trappingDate)|>
dplyr::mutate(event = 1) |>
dplyr::group_by(id_indiv) |>
dplyr::mutate(date = as.Date(trappingDate))|>
dplyr::select(-trappingDate) |>
tidyr::complete(date = seq.Date(min(date),max(date), by="day"))|>
dplyr::mutate(event = tidyr::replace_na(event,0)) |>
dplyr::group_by(id_indiv) |>
dplyr::mutate(n_days = as.integer(as.Date(date)- as.Date(min(date)))+1) |>
dplyr::left_join(rodent_data |>
dplyr::select(id_indiv,taxonomyID_field,sex,maturity), by= "id_indiv") |>
dplyr::filter(taxonomyID_field == species_selected)
8.1 Fit a cox proportional hazards model with random effects to account for individual heterogeneity
km_fit <- survfit(Surv(n_days, event) ~ sex , data = rodent_capture_date |>
dplyr::filter(sex !="na"))
coxme_model <- coxme(Surv(n_days, event) ~ sex + (1 | id_indiv), data = rodent_capture_date)
# Print summary of model
summary(coxme_model)## Mixed effects coxme model
## Formula: Surv(n_days, event) ~ sex + (1 | id_indiv)
## Data: rodent_capture_date
##
## events, n = 154, 277
##
## Random effects:
## group variable sd variance
## 1 id_indiv Intercept 1.321401 1.746102
## Chisq df p AIC BIC
## Integrated loglik 63.3 3.0 1.16e-13 57.30 48.19
## Penalized loglik 139.1 21.6 0.00e+00 95.94 30.34
##
## Fixed effects:
## coef exp(coef) se(coef) z p
## sexM 0.03047 1.03094 0.67123 0.05 0.964
## sexna 0.63380 1.88475 0.69084 0.92 0.359
8.2 Survival according to sex
sfit <- survfit(Surv(n_days, event) ~ sex, data = rodent_capture_date|>
dplyr::filter(sex !="na"))
summary(sfit)
8.3 Plot survival according to sex
ggsurvplot(sfit, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
legend.labs=c("M", "F"), legend.title="Sex",
palette=c("dodgerblue2", "orchid2"),
title="Kaplan-Meier curve for Maxomys surifer survival in relation to sex",
risk.table.height=.15)
8.3.1 Survival according to maturity
km_fit <- survfit(Surv(n_days, event) ~ maturity , data = rodent_capture_date)
coxme_model <- coxme(Surv(n_days, event) ~ maturity + (1 | id_indiv), data = rodent_capture_date)
# Print summary of model
summary(coxme_model)## Mixed effects coxme model
## Formula: Surv(n_days, event) ~ maturity + (1 | id_indiv)
## Data: rodent_capture_date
##
## events, n = 154, 277
##
## Random effects:
## group variable sd variance
## 1 id_indiv Intercept 1.322404 1.748751
## Chisq df p AIC BIC
## Integrated loglik 63.23 2.00 1.865e-14 59.23 53.15
## Penalized loglik 139.05 21.37 0.000e+00 96.31 31.41
##
## Fixed effects:
## coef exp(coef) se(coef) z p
## maturityjuvenile 0.5018 1.6516 0.5734 0.88 0.382
8.4 Plot survival according to sex
ggsurvplot(sfit, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
legend.labs=c("adult", "juvenile"), legend.title="Maturity",
palette=c("dodgerblue2", "orchid2"),
title="Kaplan-Meier curve for Maxomys surifer survival in relation to maturity",
risk.table.height=.15)
8.5 Plot survival curve for a selected individual (RM0268)
select_indiv <- "RM0268"
selected_individual <- rodent_capture_date|>
dplyr::filter(id_indiv == select_indiv)
selected_individual$event_time <- ifelse(selected_individual$event,selected_individual$n_days, NA)
selected_individual$event_time[is.na(selected_individual$event_time)] <- max(selected_individual$n_days)
selected_individual_survival <- Surv(selected_individual$event_time, selected_individual$event)
title = paste(expression("Survival curve for individual ID: "),select_indiv)
plot(survfit(selected_individual_survival ~ 1), xlab = "Time", ylab = "Survival probability", main = title)