---
title: "Getting started with FishDiveR"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with FishDiveR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 4,
  out.width = "100%"
)

# Include a figure saved by FishDiveR, copied next to the vignette's own figures
include_saved_figure <- function(file) {
  dest <- file.path(knitr::opts_chunk$get("fig.path"), basename(file))
  dir.create(dirname(dest), recursive = TRUE, showWarnings = FALSE)
  file.copy(file, dest, overwrite = TRUE)
  knitr::include_graphics(dest)
}
```

FishDiveR classifies the vertical movement behaviour of tagged aquatic animals from their depth records. Each day of a record is summarised, and days with similar diving patterns are grouped into behaviour clusters. This vignette walks through the whole workflow on the example tag that comes with the package, explaining what each step does and how to choose its settings.

If you'd rather not write code, the same workflow is available as a point-and-click app: see `vignette("app", package = "FishDiveR")`. For long or high-frequency records, also read `vignette("large-records", package = "FishDiveR")`.

## How the method works

The unit of analysis is the day: each 24-hour period, midnight to midnight in local time, becomes one row of data to cluster.

1. **Wavelet analysis.** A continuous wavelet transform finds how strongly the depth record oscillates at each period (for example every 20 minutes, every 3 hours, once a day) at every point in time. Unlike a Fourier transform, it copes with behaviour that changes through the record.
2. **Daily statistics.** The wavelet results are summarised for each day and each period (mean power, power variance, mean amplitude and so on). Separately, daily depth statistics are calculated from the depth record itself: mean and maximum depth, vertical speed, time near the surface and, optionally, differences between day and night.
3. **Principal component analysis (PCA).** Neighbouring wavelet periods carry very similar information, so the many wavelet statistics are condensed into a few principal components.
4. **K-means clustering.** The principal components and depth statistics are combined, and the days are grouped into *k* clusters of similar behaviour.

The method is described in detail in Beale et al. (2026), *Movement Ecology* 14:8, <https://doi.org/10.1186/s40462-025-00622-w>.

## Setting up

Every step saves its results in a project folder, with one sub-folder per tag (and a `Combined_tags` folder when several tags are analysed together). Later steps read the earlier results from there. Here we use a temporary folder; for your own data, use a permanent folder.

```{r setup}
# Use rgl's off-screen device, as k_clustering() can draw a 3D plot
options(rgl.useNULL = TRUE)
library(FishDiveR)

# Example data shipped with the package
example_folder <- system.file("extdata", package = "FishDiveR")

# Project folder for this vignette's results
project <- file.path(tempdir(), "FishDiveR-vignette")
dir.create(project, showWarnings = FALSE)
```

## 1. Import the depth record

`import_tag_data()` reads a csv archive with a date-time column and a depth column, crops it to the deployment, and then to whole days (midnight to midnight) in the local time zone. The cropped record is saved as `archive_days.rds` in the tag's folder.

```{r import}
archive_days <- import_tag_data(
  tag_ID = "data",
  tag_deploy_UTC = "2000-01-01 00:00:00",
  tag_release_UTC = "2000-01-11 23:59:00",
  archive = file.path(example_folder, "data", "data-Archive.csv"),
  date_time_col = 1,
  depth_col = 2,
  time_zone = "UTC",
  output = TRUE,
  output_folder = project
)
head(archive_days)
```

A few things to check with your own data:

- **Time zone.** Days are split at local midnight, so give the time zone the animal was in. Time zones with daylight saving are not supported, as the clock changes create days of 23 and 25 hours.
- **Sampling interval.** The record needs a constant sampling interval. `import_tag_data()` warns if it finds more than one interval, listing each and how many records have it.
- **Date-time format.** Dates are currently read in the format `00:01:00 01-Jan-2000`.

`plot_TDR()` plots the depth record. Its `minmax_bins` option plots only the shallowest and deepest record in each of that many time slices: much faster than plotting every record, while keeping every dive's extremes.

```{r plot-tdr, fig.alt = "Depth record of the example tag over ten days, with depth increasing downwards."}
tdr <- plot_TDR(
  rds_file = "data/archive_days.rds",
  data_folder = project,
  minmax_bins = 2000,
  Y_lim = c(0, 300, 50),
  date_breaks = "2 day"
)
attr(tdr, "TDR_plot")
```

## 2. Wavelet analysis

`create_wavelet()` runs the continuous wavelet transform over a range of periods, from `lower_period_mins` up to `upper_period_hours`, with `suboctaves` periods for each doubling of period.

- **Lower period.** The shortest cycle analysed. 5 minutes captures short dives, where the sampling interval allows; it must be at least twice the sampling interval. This example uses 30 minutes to keep the vignette quick.
- **Upper period.** The longest cycle analysed. 24 hours focuses on daily behaviour.
- **Sub-octaves.** 12 gives periods that include 24, 12 and 6 hours exactly, with 12 periods between each.

Use the same settings for all tags you will cluster together.

```{r wavelet, results = "hide"}
wavelet <- create_wavelet(
  archive = archive_days,
  tag_ID = "data",
  lower_period_mins = 30,
  upper_period_hours = 24,
  suboctaves = 12,
  output = TRUE,
  output_folder = project,
  plot_wavelet = TRUE,
  interactive_mode = FALSE
)
```

With `plot_wavelet = TRUE`, the power spectrum is saved in the tag's `2_Wavelet_Figures` folder. Red shows strong oscillation at that period and time.

```{r wavelet-figure, echo = FALSE, fig.alt = "Wavelet power spectrum of the example tag: time in days along the bottom, period from 30 minutes to 24 hours up the side."}
include_saved_figure(file.path(project, "data", "2_Wavelet_Figures", "data_wavelet_spectrum_no.p.png"))
```

Wavelets of long, high-frequency records can need a lot of memory: a year at 1-second sampling can need tens of gigabytes. `vignette("large-records", package = "FishDiveR")` explains how to estimate this and how to reduce it with `thin_archive()`.

## 3. Daily statistics

`create_wavelet_stats()` summarises the wavelet for each day and period: mean and variance of the amplitude, phase and power, and the mean squared power.

```{r wavelet-stats}
wave_stats <- create_wavelet_stats(wavelet, tag_ID = "data", output = TRUE, output_folder = project)
dim(wave_stats)
```

`create_depth_stats()` calculates daily statistics from the full-resolution depth record: mean, standard deviation, minimum and maximum depth, mean vertical speed, maximum descent and ascent speeds, skewness and kurtosis of depth, and the proportion of time near the surface.

With `diel = TRUE` it adds the differences between night and day in these statistics, standardised by each day's maximum depth. Sunrise and sunset come from a GPS file with daily positions (columns `date`, `lat` and `lon`), or from fixed `sunrise_time` and `sunset_time`.

```{r depth-stats}
depth_stats <- create_depth_stats(
  archive = archive_days,
  tag_ID = "data",
  diel = TRUE,
  GPS = file.path(example_folder, "data", "GPS.csv"),
  sunset_type = "civil",
  output = TRUE,
  output_folder = project
)
names(depth_stats)
```

## 4. Principal component analysis

`pca_data()` collects the chosen wavelet statistics for the tags to analyse. Mean power, power variance and mean amplitude are the defaults: they capture the most variance in the fewest principal components.

```{r pca-data}
pc_data <- pca_data(
  tag_vector = "data",
  data_folder = project,
  power_mean = TRUE,
  power_variance = TRUE,
  amplitude_mean = TRUE
)
```

Before choosing how many components to keep, look at the eigenvalues and the cumulative variance they explain. Asking for up to 50 components returns all of their eigenvalues.

```{r eigenvalues}
all_pcs <- pca_results(pc_data, No_pcs = 50, plot_eigenvalues = FALSE, interactive_mode = FALSE)
round(head(all_pcs$eig, 8), 2)
```

Common choices are the components with eigenvalues of 1 or more, or enough components to explain at least 70% of the variance. Here three components explain over 90%.

```{r pca}
pc_results <- pca_results(pc_data, No_pcs = 3, plot_eigenvalues = FALSE, interactive_mode = FALSE)
pc_scores <- pca_scores(pc_results, plot_loadings = FALSE)
```

`combine_data()` joins each day's principal component scores with its depth statistics, ready for clustering.

```{r combine}
kmeans_data <- combine_data(tag_vector = "data", data_folder = project, pc_scores = pc_scores)
dim(kmeans_data)
```

## 5. Choose the number of clusters and cluster the days

`select_k()` runs k-means for each number of clusters up to `Max.k`, and plots the within-cluster sum of squares (A, the "elbow" plot) and the mean silhouette width (B). Look for a number near the elbow, where more clusters stop reducing the sum of squares by much, and where the silhouette width is high. K-means starts from random cluster centres, so set a seed for reproducible results.

```{r select-k, fig.alt = "Elbow plot of within-cluster sum of squares and plot of silhouette width against the number of clusters, with a dashed line at four clusters."}
set.seed(123)
select_k(kmeans_data, Max.k = 8, v_line = 4)
```

`k_clustering()` assigns each day to one of `k` clusters. With `output = TRUE` it also saves a 3D plot of the first three principal components and figures of each cluster's mean variable values.

```{r kmeans}
set.seed(123)
kmeans_result <- k_clustering(kmeans_data, k = 4, nstart = 50, output = FALSE)
table(Cluster = kmeans_result$cluster)
```

With only ten days, these clusters are just for illustration: real analyses need many more days for meaningful clusters.

## 6. Plot the results

`plot_cluster_TDR()` plots a tag's depth record with each day coloured by its cluster.

```{r cluster-tdr, fig.alt = "Depth record of the example tag with each day coloured by its behaviour cluster."}
plot_cluster_TDR(
  tag_ID = "data",
  data_folder = project,
  kmeans_result = kmeans_result,
  minmax_bins = 2000,
  Y_lim = c(0, 300, 50),
  date_breaks = "2 day"
)
```

`plot_clusters()` plots the most typical days of each cluster: those closest to the cluster's centre. With `diel_shade = TRUE`, night is shaded grey. It saves one figure per cluster and combined figures in the tag's `6_Cluster-plots` folder.

```{r example-days}
plots <- plot_clusters(
  tag_vector = "data",
  data_folder = project,
  kmeans_result = kmeans_result,
  No_days = 1,
  Y_lim = c(0, 300, 50),
  diel_shade = TRUE,
  dpi = 72,
  output = TRUE,
  output_folder = project
)
```

```{r example-days-figure, echo = FALSE, out.width = "70%", fig.alt = "The most typical day of each of the four clusters, stacked, with night shaded grey."}
include_saved_figure(file.path(project, "data", "6_Cluster-plots.K=4_shaded", "Cluster_facet.png"))
```

## Analysing several tags together

Clustering several tags together finds behaviours shared between individuals. Run steps 1 to 3 for each tag with the same project folder and wavelet settings, then pass all the tag IDs to `pca_data()` and `combine_data()`:

```{r several-tags, eval = FALSE}
tags <- c("tag_A", "tag_B", "tag_C")
pc_data <- pca_data(tag_vector = tags, data_folder = project, output = TRUE, output_folder = project)
# ... pca_results() and pca_scores() as above ...
kmeans_data <- combine_data(tag_vector = tags, data_folder = project, pc_scores = pc_scores)
```

Results for several tags are saved in the project's `Combined_tags` folder, and `plot_clusters()` takes the same `tag_vector`. Tags analysed together should have the same wavelet settings, and ideally the same sampling interval (see `thin_archive()`).

## Citation

If you use FishDiveR, please cite:

```{r citation, echo = FALSE, comment = ""}
print(citation("FishDiveR"), style = "text")
```