## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>"
)

## ----setup--------------------------------------------------------------------
library(soReta)

head(camOp_soReta[, 1:6])   # only the first 6 of 120 days, or the table gets too wide
head(recordTable_soReta)
head(recordTableIndividuals_soReta)

## -----------------------------------------------------------------------------
# 1 row = 1 site x 1 day (no RAI column here: with n_days_active always
# equal to 1, RAI would just be N x 100 -- no extra information over N)
ds_day <- build_site_day(recordTable_soReta, camOp_soReta, threshold_min = 30)
head(ds_day[ds_day$N_sp > 0, ])

## -----------------------------------------------------------------------------
wolf_burst <- recordTable_soReta[
  recordTable_soReta$Station == "S_02" &
    recordTable_soReta$Species == "wolf" &
    format(recordTable_soReta$DateTimeOriginal, "%Y-%m-%d") == "2026-02-14",
]
wolf_burst
wolf_burst_camOp <- camOp_soReta["S_02", "2026-02-14", drop = FALSE]

build_site_day(wolf_burst, wolf_burst_camOp, threshold_min = 30, independence_method = "chain")$wolf_N
build_site_day(wolf_burst, wolf_burst_camOp, threshold_min = 30, independence_method = "window")$wolf_N

## -----------------------------------------------------------------------------
interrupted_burst <- recordTable_soReta[
  recordTable_soReta$Station == "S_04" &
    format(recordTable_soReta$DateTimeOriginal, "%Y-%m-%d") == "2026-03-10",
]
interrupted_burst

interrupted_burst_camOp <- camOp_soReta["S_04", "2026-03-10", drop = FALSE]

build_site_day(interrupted_burst, interrupted_burst_camOp, threshold_min = 30, require_uninterrupted = FALSE)$wolf_N
build_site_day(interrupted_burst, interrupted_burst_camOp, threshold_min = 30, require_uninterrupted = TRUE)$wolf_N

## -----------------------------------------------------------------------------
# 1 row = 1 site x 1 fixed N-day block
ds_week <- build_site_block(recordTable_soReta, camOp_soReta, block_days = 7, threshold_min = 30)
head(ds_week)

# 1 row = 1 site x 1 calendar month
ds_month <- build_site_month(recordTable_soReta, camOp_soReta, threshold_min = 30)
head(ds_month)

# 1 row = 1 site x 1 named period, recurring every year -- our example
# data only spans January to April 2026, so a two-period split fits
# better here than a full four-season year
ds_period <- build_site_period(
  recordTable_soReta, camOp_soReta,
  period_names  = c("early", "late"),
  period_starts = c("01/01/2026", "01/03/2026"),
  threshold_min = 30
)
head(ds_period)

## ----eval = FALSE-------------------------------------------------------------
# # example: joining a per-day covariate you already have, e.g. lunar fraction
# ds_month |> dplyr::left_join(my_lunar_fraction, by = c("mid_day" = "Date"))

## ----eval = FALSE-------------------------------------------------------------
# ds_month |> dplyr::left_join(my_station_covariates, by = "Station")

## ----eval = FALSE-------------------------------------------------------------
# recordTable_soReta |> dplyr::filter(!Species %in% c("...")) |> build_site_month(camOp_soReta, threshold_min = 30)

## -----------------------------------------------------------------------------
# 1 row = 1 site
ds_site_tot <- build_site_total(recordTable_soReta, camOp_soReta, threshold_min = 30)
head(ds_site_tot)

# 1 row = 1 day
ds_day_tot <- build_day_total(recordTable_soReta, camOp_soReta, threshold_min = 30)
head(ds_day_tot)

## -----------------------------------------------------------------------------
gr_size_ev <- build_group_size_events(recordTable_soReta, camOp_soReta, countCol = "N_individuals",
                         independence_method = "window", threshold_min = 30)
head(gr_size_ev)


## -----------------------------------------------------------------------------
# most site x day x species combinations are genuinely empty
# at this grain, so filter for the informative rows
ds_group_day <- build_group_size_day(recordTable_soReta, camOp_soReta,
                                       countCol = "N_individuals",
                                       threshold_min = 30)
ds_group_day[!is.na(ds_group_day$red_deer_mean_group_size), ] |>
  head() |>
  print(width = 90)

## -----------------------------------------------------------------------------
ds_group_total <- build_group_size_total(recordTable_soReta, camOp_soReta,
                                           countCol = "N_individuals",
                                           threshold_min = 30)
print(ds_group_total, width = 100)

## -----------------------------------------------------------------------------
occ_day <- build_occupancy_day(recordTable_soReta, camOp_soReta)
head(occ_day[["wolf"]])

occ_week <- build_occupancy_block(recordTable_soReta, camOp_soReta, block_days = 7, min_days = 4)
head(occ_week[["wolf"]])

## -----------------------------------------------------------------------------
rad <- extract_radians(recordTable_soReta, threshold_min = 30)
names(rad)
# overlap::densityPlot(rad[["wolf"]], xcenter = "midnight")

## -----------------------------------------------------------------------------
recordTable_soReta$bimonth <- paste0("bim", ceiling(lubridate::month(recordTable_soReta$DateTimeOriginal) / 2))

# one level: split by station
rad_by_station <- extract_radians(recordTable_soReta, threshold_min = 30, group_col = "Station")
rad_by_station[["wolf"]][["S_01"]]

# two levels: station, then two-month period within station
rad_by_station_bimonth <- extract_radians(recordTable_soReta, threshold_min = 30,
                                            group_col = c("Station", "bimonth"))
names(rad_by_station_bimonth[["wolf"]][["S_01"]])   # check which periods actually
                                                       # exist for this species/station
                                                       # before indexing further
rad_by_station_bimonth[["wolf"]][["S_01"]][["bim1"]]

## -----------------------------------------------------------------------------
recordTable_soReta$Cluster <- ifelse(
  recordTable_soReta$Station %in% c("S_01", "S_02"),
  "S_01_S_02_cluster",
  recordTable_soReta$Station
)

# group_cols changes what counts as independent: the two real stations
# are now merged for this purpose
rad_clustered <- extract_radians(recordTable_soReta, threshold_min = 30,
                                   group_cols = c("Cluster", "Species"))

# group_col, unchanged in meaning: still just splits the finished
# result -- here, by the same clustered column, just to display it
rad_clustered_split <- extract_radians(recordTable_soReta, threshold_min = 30,
                                         group_cols = c("Cluster", "Species"),
                                         group_col  = "Cluster")
rad_clustered_split[["wolf"]][["S_01_S_02_cluster"]]

## -----------------------------------------------------------------------------
Rad_obj <- radians_to_env(rad_by_station, prefix = "Rad_", sep = "-")
Rad_obj
str(get(Rad_obj[1]))

## -----------------------------------------------------------------------------
pair_data <- build_species_pair_intervals(
  recordTable_soReta, camOp_soReta,
  speciesA = "wolf", speciesB = "wild boar",
  threshold_min = 30
)
head(pair_data)

## -----------------------------------------------------------------------------
# 1. linear model (log-transformed, as Niedballa et al. 2019 did to meet
#    linear model assumptions)
mod <- lm(log(delta_hours) ~ direction, data = pair_data[!pair_data$censored, ])
summary(mod)

# 2. Mann-Whitney U-test 
wilcox.test(delta_hours ~ direction, data = pair_data[!pair_data$censored, ])

# 3. permutation test (shuffle species labels, keeping real timestamps
#    and each species' total count fixed, recompute the AB/BA ratio
#    each time)
set.seed(1)
n_perm <- 999
rt_pair <- recordTable_soReta[recordTable_soReta$Species %in% c("wolf", "wild boar"), ]
obs_ratio <- median(pair_data$delta_hours[pair_data$direction == "AB" & !pair_data$censored]) /
  median(pair_data$delta_hours[pair_data$direction == "BA" & !pair_data$censored])
null_ratio <- replicate(n_perm, {
  rt_perm <- rt_pair
  rt_perm$Species <- sample(rt_perm$Species)
  out_perm <- build_species_pair_intervals(rt_perm, camOp_soReta, speciesA = "wolf", speciesB = "wild boar",
                                            threshold_min = 30)
  median(out_perm$delta_hours[out_perm$direction == "AB" & !out_perm$censored]) /
    median(out_perm$delta_hours[out_perm$direction == "BA" & !out_perm$censored])
})
mean(null_ratio >= obs_ratio, na.rm = TRUE)  # empirical p-value

## -----------------------------------------------------------------------------
interruptions <- build_species_pair_interruptions(
  recordTable_soReta,
  speciesA = "wolf", speciesB = "wild boar",
  threshold_min = 30
)
table(interruptions$type)

# permutation test on AA vs BB (same logic as above, applied to the
# interruption-type intervals instead of AB/BA)
obs_ratio_aabb <- median(interruptions$delta_hours[interruptions$type == "AA"]) /
  median(interruptions$delta_hours[interruptions$type == "BB"])
null_ratio_aabb <- replicate(n_perm, {
  rt_perm <- rt_pair
  rt_perm$Species <- sample(rt_perm$Species)
  out_perm <- build_species_pair_interruptions(rt_perm, speciesA = "wolf", speciesB = "wild boar",
                                                threshold_min = 30)
  median(out_perm$delta_hours[out_perm$type == "AA"]) /
    median(out_perm$delta_hours[out_perm$type == "BB"])
})
mean(null_ratio_aabb >= obs_ratio_aabb, na.rm = TRUE)

## ----eval = requireNamespace("GLMMadaptive", quietly = TRUE) && requireNamespace("mgcv", quietly = TRUE)----
diel_month <- build_diel_binomial_month(recordTable_soReta, camOp_soReta, bin_hours = 1, min_days = 10)
diel_wolf  <- diel_month[diel_month$Species == "wolf", ]
diel_wolf$Station <- factor(diel_wolf$Station)

# trigonometric GLMM (Iannarilli et al. 2024, section 3.3)
trig_model <- GLMMadaptive::mixed_model(
  fixed  = cbind(success, failure) ~ cos(2 * pi * Time / 24) + sin(2 * pi * Time / 24) +
                                      cos(2 * pi * Time / 12) + sin(2 * pi * Time / 12),
  random = ~ 1 | Station,
  data   = diel_wolf,
  family = binomial()
)
summary(trig_model)

## -----------------------------------------------------------------------------
ch <- build_cmr_block(recordTableIndividuals_soReta, camOp_soReta, block_days = 7, min_days = 4)
as_capture_strings(ch)