# forceR

## Introduction

The package `forceR` has originally been written for insect bite force data preparation and analysis, but it can be used for any kind of time series measurements. Functions include

- loading, plotting, and cropping of data
- correction of charge amplifier drifts
- correction of baseline drifts
- reduction of sampling frequency
- automatic extraction of single peaks
- rescaling (normalization) of curves
- reduction of curves to 100 time steps each
- finding of best polynomial fits to describe all curves

This vignette guides you through all functions of the package. A limited data set of 8 bite force measurement files is used to get the ideas across while allowing for fast calculations.

At the end of this vignette you will find a self-sufficient workflow example that runs all `forceR` commands after file loading and drift corrections. Instead of loading files, bite series are simulated and then analyzed.

Please cite the following publication when using the `forceR` package:

Rühr, PT & Blanke, A (2022): `forceX` and `forceR`: a mobile setup and R package to measure and analyze a wide range of animal closing forces. Methods in Ecology and Evolution XXX: pp.XX-XX. doi: 10.1111/2041-210X.13909

## Installation

### Official release

You can install the official release of `forceR` from CRAN:

    install.packages('forceR')

### Development version

The development version of `forceR` is available at GitHub:

    require(devtools)
    devtools::install_github("https://github.com/Peter-T-Ruehr/forceR")

## Detailed example workflow: Data loading and preparation

### Load packages to work with this vignette, including `forceR`

    library(magrittr)
    library(dplyr)
    library(stringr)
    library(purrr)
    library(ggplot2)
    library(readr)

    library(forceR)

### Download raw measurement files

From hereon, we will use data stored in a local folder. Please download and unzip the content of the example_data.zip file from Zenodo onto your hard drive, unzip the files and store the folder location in the variable `data.folder`. This folder now contains every file that is needed and has been produced during the creation of this vignette, so you can always restore the files from `Zenodo` in case you have accidentally overwritten anything you wanted to keep.

If you decide to use your own files, consider that `forceR` generally expects file names to start with leading digits specifying the measurement number (E.g. `0001_G_maculatus.csv`). The number (in this case `0001`) is used to keep data files, log files, and PDF files of the same measurement associated with each other. If not already the case, we strongly recommend renaming the raw data files accordingly before continuing with this `forceR` vignette.

First, we set the folder containing the example data as the `data.folder` (make sure this fits your local path, e.g. `"./example_data"`):

    data.folder <- "./example_data"

In case you want to use your own measurements, and if they were recorded by the LJStream software (LabJack Corporation, Lakewood, Colorado, US; tested for v1.19) as suggested in the forceX instructions, you can use the `convert_measurement()` function of `forceR` to convert your files and save them in the `data.folder`. See `?convert_measurement()` for details.

### Plot raw measurement

The measurement file should have the time steps (preferably in m.secs, since all code of this package expects this) in the first column, and the measurements in the second column. Column names do not matter, and everything but the first two columns will be ignored.

    file <- file.path(data.folder, "0982.csv")
    plot_measurement(file,
                     columns = c(1:2))

### Crop raw measurement

The function `crop_measurement()` can be used to crop a measurement. Since we want to work with the cropped file later, we define `path.data` when calling the function to store the result as a new file. Before running the function, create the folder on your hard drive, e.g. at `"./example_data/cropped"`.

    cropped.folder <- "./example_data/cropped"

If a measurement contains two distinct regions of bites with a lot of unnecessary data and/or measurement artefacts in-between (such as user-made peaks), I recommend to manually copy the RAW data files, give the copy a new measurement number (as if it was actually a separate measurement), and crop the distinct parts containing actual bites separately from the two copies of the file. For more distinct regions, create more copies.

I recommend to not crop the files too much in case baseline corrections are needed later, because then the `baseline_corr()` function will not be able to figure out where the actual baseline might be. Leaving several seconds before and after the first and last bite of a series will prevent such problems.

    file.cropped <- crop_measurement(file,
                                     path.data = cropped.folder)

Select two points at the beginning and end of the part of the measurement we want to keep. If more than two points are selected, only the last two points will be considered. This allows for corrections of erroneous clicks. The position of the points on the y-axis is irrelevant. Make sure that the start and end points are placed on the graph between peaks, i.e. on the baseline. The y-values of the graph at these positions will be used as reference points in the following `amp_drift_corr()` function.

The resulting file will be saved with the suffix `"_cropped"`. We can plot it again:

    cropped.folder <- "./example_data/cropped"

    file <- file.path(cropped.folder, "0982_cropped.csv")
    plot_measurement(file)

### Amplifier drift correction

To remove the systemic, asymptotical drift characteristic for charge amplifiers with resistor–capacitor (RC) circuits, we can use the `amp_drift_corr()` function. This function requires a few parameters. `?amp_drift_corr` brings up the full explanation of all values.

Basically, we only have to decide how we want to plot the data, define the correct folder with our measurement data, provide the time constant (tau, \(\tau\)) of our charge amplifier, and define a folder where we want to store the results. This function, in contrast to the previous functions, takes a folder as input because when the same charge amplifier was used for all measurements, the same settings can be used for all measurement files.

    # create file list with all cropped files that need amplifier drift corrections.
    #   Here we use the cropped.folder with the cropped measurements.
    file.list <- list.files(cropped.folder,
                            pattern = "csv$",
                            full.names = TRUE)

    # define folder where to save the amplifier drift corrected file.
    #   If this folder does not yet exist, it will be created.
    ampdriftcorr.folder <- "./cropped/ampdriftcorr"

    for(filename in file.list){
      print(filename)
      amp_drift_corr(filename = filename,
                     tau = 9400,
                     res.reduction = 10,
                     plot.to.screen = FALSE,
                     write.data = TRUE,
                     write.PDFs = TRUE,
                     write.logs = TRUE,
                     output.folder = ampdriftcorr.folder,
                     show.progress = FALSE)
      print("***********")
    }

Because we set the arguments `write.data`, `write.PDFs` and `write.logs` to `TRUE`, this creates the `output.folder` (`"./ampdriftcorr"`) and two new folders inside the `output.folder`:

- `"./ampdriftcorr/logs/"` where log files containing info on the `amp_drift_corr` used are stored
- `"./ampdriftcorr/pdfs/"` where the PDFs of the time series before and after the corrections are stored.

The resulting graph of one corrected measurement may look like this:

In gray, we see the raw measurements, and in green the corrected data. A linear drift in the corrected data may occur in case the first value of the measurement is not exactly `0`. This drift is also corrected for by the function. The final result for each file is, since `write.data` was set to `TRUE`, also saved in a new CSV file in the `"./ampdriftcorr/"` folder.

### Automatic or manual baseline correction of time series

If the baseline (zero-line) of a measurement is unstable (e.g. due to temperature fluctuations, wind, …), it needs to be continually adjusted throughout the measurement. Ideally, this is done automatically by the function `baseline.corr()`. This automatic adjustment of the baseline invokes a sliding window approach, during which the ‘minimum’ within each window is stored. A ‘minimum’ is defined by the parameter `quantile.size`. If `quantile.size` is set to `0.05`, the value below which 5% of the measurement data within the current window lies, is treated as the current window’s ‘minimum’. Not taking the actual minimum within each window prevents the treatment of short undershoots (artifacts created during the measurement) as minima. In a second iteration, another sliding window calculates the average of the ‘minima’ within each window. The resulting ‘minimal’ values are subtracted from the original time series.

This approach works well for time series with relatively short peaks. For longer peaks, the `window.size.mins` should be increased. However, the greater `window.size.mins` gets, the smaller becomes the baseline-correcting effect of the function.

For the file `"./example_data/cropped/ampdriftcorr/1068_ampdriftcorr.csv"`, we choose a `window.size.mins` of 2000:

    filename = file.path(ampdriftcorr.folder, "1068_ampdriftcorr.csv")

    plot_measurement(filename)

    baselinecorr.folder <- "./cropped/ampdriftcorr/baselinecorr"

    file.baselinecorr <- baseline_corr(filename = filename,
                                       corr.type = "auto",
                                       window.size.mins = 2000,
                                       window.size.means = NULL,
                                       quantile.size = 0.05,
                                       y.scale = 0.5,
                                       res.reduction = 10,
                                       Hz = 100,
                                       plot.to.screen = TRUE,
                                       write.data = TRUE,
                                       write.PDFs = TRUE,
                                       write.logs = TRUE,
                                       output.folder = baselinecorr.folder,
                                       show.progress = FALSE)

Again, new folders, PDFs and log files are created. The resulting plot shows the original data in gray, the ‘minima’ of the sliding windows in blue, the sliding mean of the ‘minima’ in bright green and the corrected time series in dark green:

If the automatic approach does not yield acceptable results, an interactive manual approach to correct the baseline can be performed instead. We do the manual correction for the file `"./example_data/cropped/ampdriftcorr/1174_ampdriftcorr.csv"`, which contains measurements with long, plateau-like peaks where the automatic baseline correction would fail.

    filename = file.path(ampdriftcorr.folder, "1174_ampdriftcorr.csv")

    plot_measurement(file)

    file.baselinecorr <- baseline_corr(filename = filename,
                                       corr.type = "manual",
                                       plot.to.screen = TRUE,
                                       write.data = TRUE,
                                       write.PDFs = TRUE,
                                       write.logs = TRUE,
                                       output.folder = baselinecorr.folder,
                                       show.progress = FALSE)

The green line represents a spline along the manually defined points and is subtracted from the original data (gray). The corrected data is plotted in dark green again.

We have used `baseline_corr()` on the files `1041_ampdriftcorr.csv`, `1068_ampdriftcorr.csv`, `1174_ampdriftcorr.csv`, `1440_ampdriftcorr.csv`.

### File sorting

After all original files have been processed, we need to separate the resulting files that we want to keep from files that have been created along the way. The function `sort_files()` looks through a vector of folders (`data.folders`) to identify files with the same measurement name (i.e., the same character string in the file names before the first underscore (`_`) and keeps the version of the file that is found in the folder that is most closely to the end of the list of `data.folders`. Thus, it is important that the folders are parsed in the correct order. In our example, this order would look like this:

1. `"./example_data/"`
2. `"./example_data/cropped/"`
3. `"./example_data/cropped/ampdriftcorr/"`
4. `"./example_data/cropped/ampdriftcorr/baselinecorr"`

If a file is present in `"./example_data/cropped/ampdriftcorr/baselinecorr"`, this will be copied (`move = FALSE`) or moved (`move = TRUE`) to a folder specified in `path.corrected`, e.g., `"./example_data/corrected/"`. Versions of the same measurement in the other folders will be ingored. This is iteratively repeated for the rest of the folders in reverse order of the `data.folders` vector.

    data.folders <- c(data.folder,
                      file.path(data.folder, "/cropped"),
                      file.path(data.folder, "/cropped/ampdriftcorr"),
                      file.path(data.folder, "/cropped/ampdriftcorr/baselinecorr"))

    results.folder <- file.path(data.folder, "/corrected/")

    sort_files(data.folders = data.folders,
               results.folder = results.folder,
               move = FALSE)

### File loading

Now we will load the first file in the `results` folder using `load_single()`. The function removes all columns except the one defined in `columns = c(1:2)` and renames them to `t` and `y`. Then, it adds a `filename` column based on the file name of the data.

    file.list <- list.files(results.folder, pattern = "csv", full.names = TRUE)
    df.1 <- load_single(file = file.list[1],
                        columns = c(1:2))

    df.1
    #> # A tibble: 129,001 x 3
    #>        t       y filename
    #>    <dbl>   <dbl> <chr>
    #>  1     0 0.00447 0979_ampdriftcorr
    #>  2     1 0.00954 0979_ampdriftcorr
    #>  3     2 0.0146  0979_ampdriftcorr
    #>  4     3 0.0146  0979_ampdriftcorr
    #>  5     4 0.00951 0979_ampdriftcorr
    #>  6     5 0.00951 0979_ampdriftcorr
    #>  7     6 0.00950 0979_ampdriftcorr
    #>  8     7 0.00949 0979_ampdriftcorr
    #>  9     8 0.00948 0979_ampdriftcorr
    #> 10     9 0.00948 0979_ampdriftcorr
    #> # ... with 128,991 more rows

    unique(df.1$filename)
    #> [1] "0979_ampdriftcorr"

Alternatively, We can rapidly load all files in the `results` folder by using `load_mult()`:

    df.all <- load_mult(folder = results.folder,
                        columns = c(1:2),
                        show.progress = TRUE)

    df.all
    #> # A tibble: 615,724 x 3
    #>        t       y filename
    #>    <dbl>   <dbl> <chr>
    #>  1     0 0.00447 0979_ampdriftcorr
    #>  2     1 0.00954 0979_ampdriftcorr
    #>  3     2 0.0146  0979_ampdriftcorr
    #>  4     3 0.0146  0979_ampdriftcorr
    #>  5     4 0.00951 0979_ampdriftcorr
    #>  6     5 0.00951 0979_ampdriftcorr
    #>  7     6 0.00950 0979_ampdriftcorr
    #>  8     7 0.00949 0979_ampdriftcorr
    #>  9     8 0.00948 0979_ampdriftcorr
    #> 10     9 0.00948 0979_ampdriftcorr
    #> # ... with 615,714 more rows

    unique(df.all$filename)
    #> [1] "0979_ampdriftcorr"
    #> [2] "0980_ampdriftcorr"
    #> [3] "0981_ampdriftcorr"
    #> [4] "0982_ampdriftcorr"
    #> [5] "1041_baselinecorr"
    #> [6] "1068_baselinecorr"
    #> [7] "1174_baselinecorr"
    #> [8] "1440_baselinecorr"

The new tibble contains the measurement of all files in the selected folder.

## Detailed example workflow: Data analysis

To allow users to test the package and follow the analysis steps of this vignette without having to load actual files, we have simulated data using the `forceR` function `simulate_bites()` and stored them within `forceR`. We will replace our previous `df.all` tibble with the internal `forceR` tibble. You may, however, continue with your just created version of `df.all`.

    # create/replace df.all
    df.all <- forceR::df.all
    head(df.all)
    #> # A tibble: 6 × 3
    #>       t     y measurement
    #>   <int> <dbl> <chr>
    #> 1     1     0 m_01
    #> 2     2     0 m_01
    #> 3     3     0 m_01
    #> 4     4     0 m_01
    #> 5     5     0 m_01
    #> 6     6     0 m_01

Now let us plot the simulated measurements ot get a feeling for them:

    # plot simulated measurements
    ggplot(df.all,
           aes(x = t ,
               y = y,
               colour=measurement)) +
      geom_line()

We can see that some bites have more (reddish colors) or less (blueish colors) plateau-like peaks, whereas others have more sinusoidal peaks with early (orange) or late peaks (green).

### Sampling Frequency Reduction

`df.all` has a sampling frequency of 1000 Hz. This is more than is often needed to analyse a time series. To reduce the sampling frequency of all measurements that are above a desired frequency of e.g. `200 Hz`, we can use the function `reduce.frq()`. However, please not that the simulated bites are only 50 ms long, which is much shorter than, e.g., actual insect bites. This will keep all calculations quick, but will result in measurement series that are more downsampled than we would recommend. For real insect bites, though, a sampling frequency of 200 Hz is sufficient, so we will proceed as if we were working with such real bites:

    # reduce frequency to 200 Hz
    df.all.200 <- reduce_frq(df = df.all,
                             Hz = 200,
                             measurement.col = "measurement")

    head(df.all.200)
    #> # A tibble: 6 × 3
    #> # Groups:   t [6]
    #>       t      y measurement
    #>   <dbl>  <dbl> <chr>
    #> 1     0 0      m_01
    #> 2     5 0      m_01
    #> 3    10 0.0439 m_01
    #> 4    15 0.299  m_01
    #> 5    20 0.644  m_01
    #> 6    25 0.821  m_01

If `measurement.col` is not defined, the whole input data frame will be treated as if it was just one single time series. This is okay for data frames like `df.1` that only contain one measurement, but for data frames with multiple time series measurements (like `df.all`), we need to define a grouping column - in our case the column `measurement`.

### Creating a classifier with dataset info

For further analyses, we need a `classifier` in which data on the species, specimens and measurements of our dataset are stored:

    # create a classifier
    number_of_species <- 4
    number_of_specimens_per_species <- 3
    number_of_measurements_per_specimen <- 2
    number_of_rows <- number_of_species *
      number_of_specimens_per_species *
      number_of_measurements_per_specimen

    species <- sort(rep(paste0("species_", LETTERS[1:number_of_species]),
                        length=number_of_rows))

    specimens <- sort(rep(paste0("speciemen_", letters[1:(number_of_species*number_of_specimens_per_species)]),
                          length=number_of_rows))

    classifier <- tibble(species = species,
                         specimen = specimens,
                         measurement = paste0("m_",  str_pad(string= 1:number_of_rows, width = 2, pad = "0")),
                         amp = c(rep(0.5, number_of_rows/2), rep(2, number_of_rows/2)),
                         lever.ratio = rep(0.5, number_of_rows))
    head(classifier)
    #> # A tibble: 6 × 5
    #>   species   specimen    measurement   amp lever.ratio
    #>   <chr>     <chr>       <chr>       <dbl>       <dbl>
    #> 1 species_A speciemen_a m_01          0.5         0.5
    #> 2 species_A speciemen_a m_02          0.5         0.5
    #> 3 species_A speciemen_b m_03          0.5         0.5
    #> 4 species_A speciemen_b m_04          0.5         0.5
    #> 5 species_A speciemen_c m_05          0.5         0.5
    #> 6 species_A speciemen_c m_06          0.5         0.5

### Convert measurement to force

To convert a voltage measurement to a force measurement [Newton], we need an amplification value and, depending on the measurement setup, the lever ratio of the lever forwarding the force from the force plate to the sensor.

Now we can convert the y-values of `df.all.200` to force data and add the `specimen` and `species` info from the `classifier` using the `y_to_force()` function:

    df.all.200.tax <- y_to_force(df = df.all.200,
                                 classifier = classifier,
                                 measurement.col = "measurement")
    head(df.all.200.tax)
    #> # A tibble: 6 × 4
    #> # Groups:   t [6]
    #>   measurement specimen        t  force
    #>   <chr>       <chr>       <dbl>  <dbl>
    #> 1 m_01        speciemen_a     0 0
    #> 2 m_01        speciemen_a     5 0
    #> 3 m_01        speciemen_a    10 0.0439
    #> 4 m_01        speciemen_a    15 0.299
    #> 5 m_01        speciemen_a    20 0.644
    #> 6 m_01        speciemen_a    25 0.821

### Minimum, maximum and standard deviation of force data per measurement and species

#### Initial reduction to one row per specimen

Finally, we can start with the actual analysis of the data. First, we want to calculate the minimum, maximum and standard deviation of force for each measurement and, in case there were several measurements per specimen, also for each specimen.

From the previous steps, df.all.200.tax already contains the columns `measurement`, `specimen` and `species`. Now, we can use the function `summarize_measurements()` to create a summary tibble by adding the column names for which we want a summary in `var1` and `var2`. `var1` must name the column that stores the unique measurement IDs, e.g. measurement number, to calculate minimal and maximal force values per measurement. As `var2` we will use ‘specimen’ to calculate the summary data per specimen from the min. and max. values of the measurements of that specimen. Thus, the resulting tibble will contain summaries of all measurements that were performed with the same specimen.

    # add species to df.all.200.tax
    df.all.200.tax <- df.all.200.tax %>%
      left_join(classifier %>%
                  select(species, measurement))

    var1 = "measurement"
    var2 = "specimen"
    df.summary.specimen <- summarize_measurements(df.all.200.tax,
                                                  var1,
                                                  var2)
    head(df.summary.specimen)
    #> # A tibble: 6 × 8
    #>   measurement specimen    species   max.F.meas…¹ mean.…² max.F…³ sdv.m…⁴ n.mea…⁵
    #>   <chr>       <chr>       <chr>            <dbl>   <dbl>   <dbl>   <dbl>   <int>
    #> 1 m_01        speciemen_a species_A         1.24    1.24    1.24 0.00116       2
    #> 2 m_02        speciemen_a species_A         1.24    1.24    1.24 0.00116       2
    #> 3 m_03        speciemen_b species_A         1.13    1.13    1.13 0.00658       2
    #> 4 m_04        speciemen_b species_A         1.12    1.13    1.13 0.00658       2
    #> 5 m_05        speciemen_c species_A         1.09    1.18    1.26 0.117         2
    #> 6 m_06        speciemen_c species_A         1.26    1.18    1.26 0.117         2
    #> # … with abbreviated variable names ¹​max.F.measurement, ²​mean.F.specimen,
    #> #   ³​max.F.specimen, ⁴​sdv.max.F.specimen, ⁵​n.measurements.in.specimen

    # boxplot of maximum force in specimens
    ggplot(data = df.summary.specimen, mapping = aes(x=specimen,y=max.F.measurement)) +
      geom_jitter(aes(color='blue'),alpha=0.7, width = 0.2, height = 0.0) +
      geom_boxplot(fill="bisque",color="black",alpha=0.3) +
      # scale_y_log10() +
      labs(y="max(F)/specimen") +
      guides(color="none") +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust=1))

#### Further summary calculations for species-wise info

To get a summary on species-wise info, we have to calculate with the summarized specimen info. We are not using the `summarize_measurements()` functions because this would ignore the fact the some measurements of one species may come from the same specimen, but we only want to consider one maximum force value per specimen and not one per measurement.

    df.summary.species <- df.summary.specimen %>%
      # find max Fs of species
      group_by(species) %>%
      # calculate force values for each species
      mutate(max.F.species = max(max.F.specimen),
             mean.F.species = round(mean(max.F.specimen),6),
             sdv.max.F.species = sd(max.F.specimen)) %>%
      ungroup() %>%
      # count specimens / species
      group_by(species) %>%
      mutate(n.specimens.in.species = length(unique(specimen))) %>%
      ungroup()
    df.summary.species
    #> # A tibble: 24 × 12
    #>    measurement specimen  species max.F…¹ mean.…² max.F…³ sdv.m…⁴ n.mea…⁵ max.F…⁶
    #>    <chr>       <chr>     <chr>     <dbl>   <dbl>   <dbl>   <dbl>   <int>   <dbl>
    #>  1 m_01        specieme… specie…    1.24    1.24    1.24 0.00116       2    1.26
    #>  2 m_02        specieme… specie…    1.24    1.24    1.24 0.00116       2    1.26
    #>  3 m_03        specieme… specie…    1.13    1.13    1.13 0.00658       2    1.26
    #>  4 m_04        specieme… specie…    1.12    1.13    1.13 0.00658       2    1.26
    #>  5 m_05        specieme… specie…    1.09    1.18    1.26 0.117         2    1.26
    #>  6 m_06        specieme… specie…    1.26    1.18    1.26 0.117         2    1.26
    #>  7 m_07        specieme… specie…    2.19    2.10    2.19 0.125         2    2.44
    #>  8 m_08        specieme… specie…    2.01    2.10    2.19 0.125         2    2.44
    #>  9 m_09        specieme… specie…    2.27    2.28    2.29 0.0125        2    2.44
    #> 10 m_10        specieme… specie…    2.29    2.28    2.29 0.0125        2    2.44
    #> # … with 14 more rows, 3 more variables: mean.F.species <dbl>,
    #> #   sdv.max.F.species <dbl>, n.specimens.in.species <int>, and abbreviated
    #> #   variable names ¹​max.F.measurement, ²​mean.F.specimen, ³​max.F.specimen,
    #> #   ⁴​sdv.max.F.specimen, ⁵​n.measurements.in.specimen, ⁶​max.F.species

    # boxplot of maximum force in species
    ggplot(data = df.summary.species, mapping = aes(x=species,y=max.F.specimen)) +
      geom_jitter(aes(color='blue'),alpha=0.7, width = 0.2, height = 0.0) +
      geom_boxplot(fill="bisque",color="black",alpha=0.3) +
      # scale_y_log10() +
      labs(x='species', y="max(F)/specimen") +
      guides(color="none") +
      theme_minimal()

### Identify bites

Now we come to one of the core functions of `forceR`: `find_strongest_peaks()`. This function identifies the strongest peaks of each measurement. But first, we need to create some folders to store the data and plots that are produced in all subsequent analyses of this example:

    # create folders to save df and results
    path.plots <- paste0(data.folder, "/plots/")
    ifelse(!dir.exists(path.plots), dir.create(path.plots), "./plots already exists")

    path.plots.initial_peak_finding <- paste0(data.folder, "/plots/initial_peak_finding/")
    ifelse(!dir.exists(path.plots.initial_peak_finding), dir.create(path.plots), "./plots/initial_peak_finding already exists")

    path.data <- paste0(data.folder, "/data/")
    ifelse(!dir.exists(path.data), dir.create(path.data), "./data already exists")

Now we let `find_strongest_peaks()` identify all actual peak curves in all measurements:

    peaks.df <- find_strongest_peaks(df = df.all.200.tax,
                                     no.of.peaks = 5,
                                     path.data = path.data,
                                     path.plots = path.plots,
                                     show.progress = TRUE)

    head(peaks.df)
    #> # A tibble: 4 × 4
    #>   species   measurements                 starts                 ends
    #>   <chr>     <chr>                        <chr>                  <chr>
    #> 1 species_A m_06; m_01; m_02; m_01; m_06 265; 265; 265; 335; 70 325; 325; 325; …
    #> 2 species_B m_11; m_11; m_12; m_10; m_09 140; 5; 205; 5; 340    200; 65; 265; 6…
    #> 3 species_C m_17; m_16; m_18; m_16; m_18 5; 70; 70; 200; 5      65; 130; 130; 2…
    #> 4 species_D m_19; m_22; m_24; m_21; m_20 135; 0; 265; 330; 330  200; 65; 330; 3…

The initial peak finding in the first iteration of the function `find_strongest_peaks()` identifies peaks via a the `initial.threshold` which is by default set to `0.05` (= 5%) of the maximum force of the respective measurement. The threshold is indicated as a green line in the plots. Bite starts are indicated as blue points, peak ends as orange points:

In a second iteration, the function `find_strongest_peaks()` optimizes the peak starts and ends found in the first iteration. Starting from each start, a sliding window of size `slope.length.start` (in time steps) moves backwards in time and calculates the slope within the current window. The parameter `slope.thresh.start` defines below which slope the process is stopped. The time point where the threshold is reached is treated as the actual start of that peak. The same is done in forward direction with the peak ends. Here, `slope.length.end` (in time steps) defines the sliding window size, and `slope.thresh.end` defines the slope threshold. All these settings have default values which proved as good choices for insect bites. The peak optimization will not be done for all peaks, but only for the strongest peaks per species (not per measurement!), specified by `no.of.peaks`. If fewer than `no.of.peaks` peaks are found for a species, the results table will contain fewer peaks for this species. No warning is issued.

The results are not automatically plotted. This can be done with the function `plot_peaks()`:

    plot_peaks(df.peaks = peaks.df,
               df.data = df.all.200.tax,
               additional.msecs = 2000,
               path.plots = path.plots)

The first page of the resulting PDF looks like this, showing the five strongest peaks of `species_A`, which was measured in measurements `m_01 - m_06`, and the strongest peak of `species_A` which occured in measurement `m_11`. Plots of further peaks of `species_A` and the peaks of the other species follow on the next pages.

We strongly recommend saving `df.peaks` as a csv file on your computer and keep a copy of it save. The following procedures overwrite values in this table, so the peak-finding step needs to be repeated to restore original values.

### Manual peak corrections

For some measurements, the automatic finding identification of peak starts and ends may fail. Changing the parameters of `find_strongest_peaks()` might help, but in some cases, the starts and ends have to be manually defined. For these cases, the function `correct_peak()` has been written. Set `path.data` to the folder where you stored `peaks.df` before. This will overwrite the info of the bite you want to correct (here: `bite 1` of measurement `m_01`!

    peaks.df <- correct_peak(df.peaks = peaks.df,
                             df.data = df.all.200.tax,
                             measurement = "m_01",
                             peak = 1,
                             additional.msecs = 100,
                             path.data = path.data)

    [1] "Select new peak start and end. If more than two points are selected, the operation is automatically terminated."

This functions prompts the user to define a new start and end of the peak. Define the `measurement` ID (as a character string) in which the peak to correct is located, and the `peak` number (as a numerical value). Both can be found in the respective title of the peak curve plot created by `plot_peaks()`. If the desired start and/or end lies outside of the plot window, increase `additional.msecs`.

The function overwrites the start and end values of the selected peak in the selected measurement and plots the corrected peak curve. Again, note that this example is taken from real measurements and not the simulated data we use in the other analysis steps.

Note that the peak in the figures above looks differently from the one in your plots, but the principles are the same.

### Normalize (rescale) peaks

In the next step, we will use `rescale_peaks()` to rescale the time and force data so that they both range from 0 to 1. This function takes the data frame that contains the individual peak starts and ends, and the measurement data. Both data frames are linked by the ‘measurement’ column, so this column needs to be present in both data frames.

    peaks.df.norm <- rescale_peaks(df.peaks = peaks.df,
                                   df.data = df.all.200.tax,
                                   plot.to.screen = TRUE,
                                   path.data = path.data,
                                   show.progress = TRUE)

    head(peaks.df.norm)
    #> # A tibble: 6 × 6
    #>   measurement  peak t.norm force.norm specimen    species
    #>   <chr>       <int>  <dbl>      <dbl> <chr>       <chr>
    #> 1 m_06            1 0          0      speciemen_c species_A
    #> 2 m_06            1 0.0833     0      speciemen_c species_A
    #> 3 m_06            1 0.167      0.0764 speciemen_c species_A
    #> 4 m_06            1 0.25       0.467  speciemen_c species_A
    #> 5 m_06            1 0.333      0.836  speciemen_c species_A
    #> 6 m_06            1 0.417      1      speciemen_c species_A

Now we will plot the resulting data. You will note that the peaks are not very smooth. This is because the simulated peaks of `df.all` where so short that the reduction to 200 Hz in one of the previous steps resulted in very few time steps per bite. It is, thus, important to choose the sampling frequency wisely. For our example, however, these few time steps are sufficient.

    # plot all normalized peaks
    ggplot(peaks.df.norm %>%
             mutate(color.column = paste0(measurement, "__bite_", peak)),
           aes(x = t.norm ,
               y = force.norm,
               colour=color.column)) +
      geom_line()

### Reduce to 100 time steps per peak

`peaks.df.norm` now contains the same data as `df.all.200.tax`, only with the additional columns of `t.norm` and `force.norm` which contain the rescaled time and force data, both ranging from 0 to 1.

To reduce the data to 100 observations (time steps), we will use the function `red_peaks_100()`. Note that, because the peaks of the example dataset are so short, we now actually increase the number of time steps. In real data, however, `red_peaks_100()` usually boils down the number of observations.

    peaks.df.norm.100 <- red_peaks_100(df = peaks.df.norm,
                                       plot.to.screen = TRUE,
                                       path.data = path.data,
                                       path.plots = path.plots,
                                       show.progress = TRUE)
    head(peaks.df.norm.100)

    head(peaks.df.norm.100)
    #> # A tibble: 6 × 6
    #>   species   measurement specimen     peak index force.norm.100
    #>   <chr>     <chr>       <chr>       <int> <dbl>          <dbl>
    #> 1 species_A m_06        speciemen_c     1     1        0
    #> 2 species_A m_06        speciemen_c     1     2        0.00568
    #> 3 species_A m_06        speciemen_c     1     3        0.00891
    #> 4 species_A m_06        speciemen_c     1     4        0.0101
    #> 5 species_A m_06        speciemen_c     1     5        0.00972
    #> 6 species_A m_06        speciemen_c     1     6        0.00814

Again, we will plot the result:

    head(peaks.df.norm.100)
    #> # A tibble: 6 × 6
    #>   species   measurement specimen     peak index force.norm.100
    #>   <chr>     <chr>       <chr>       <int> <dbl>          <dbl>
    #> 1 species_A m_06        speciemen_c     1     1        0
    #> 2 species_A m_06        speciemen_c     1     2        0.00568
    #> 3 species_A m_06        speciemen_c     1     3        0.00891
    #> 4 species_A m_06        speciemen_c     1     4        0.0101
    #> 5 species_A m_06        speciemen_c     1     5        0.00972
    #> 6 species_A m_06        speciemen_c     1     6        0.00814

    # plot normalized peaks: 5 bites per measurement
    ggplot(peaks.df.norm.100 %>%
             mutate(color.column = paste0(measurement, "-bite", peak)),
           aes(x = index ,
               y = force.norm.100,
               colour=color.column)) +
      geom_line()

The resulting tibble consists of 100 rows per peak, and within each peak, the `index` column runs from 1 to 100, and the `force.norm.100` column consists of the rescaled peak curves, each of whose values range from 0 to 1. Since we included 4 species in our dataset, the resulting tibble is 2,000 rows long (`species * peaks/species * time steps per peak = 4 * 5 * 100 = 2,000`).

### Average peak curve per species

Next, we will calculate the average peak curve per species and normalize the curves again:

    peaks.df.100.avg <- avg_peaks(df = peaks.df.norm.100,
                                  path.data = path.data)

    head(peaks.df.100.avg)
    #> # A tibble: 6 × 3
    #> # Groups:   species [1]
    #>   species   index force.norm.100.avg
    #>   <chr>     <dbl>              <dbl>
    #> 1 species_A     1           0.00531
    #> 2 species_A     2           0.00267
    #> 3 species_A     3           0.000855
    #> 4 species_A     4           0
    #> 5 species_A     5           0.000247
    #> 6 species_A     6           0.00173

    # plot averaged normalized curves per species
    ggplot(peaks.df.100.avg, aes(x = index ,
                                 y = force.norm.100.avg,
                                 colour=species)) +
      geom_line()

### Fit polynomial functions to the force curves

To find out the minimum number of coefficients that well describe all curves, we use the function `find_best_fits()`. It fits polynomial functions with 1 to 20 coefficients and uses the Akaike Information Criterion (AIC) to evaluate the goodness of the fits. A model is considered a good fit when the percentage of change from one model to the next (e.g. a model with 6 coefficients to a model with 7 coefficients) is `<= 5%`. The first four models meeting this criterion are plotted as colored graphs and the AICs of these models are visualized in a second plot for each curve. All first four coefficients per curve that fulfill the criterion are stored and in the end, a histogram of how often which coefficients were good fits is plotted as well. The function returns the numerical value of the coefficient that fulfilled the criterion of a good fit in most curves. If `write.PDFs == TRUE`, then the plots are saved as PDFs in `path.plots`. Resulting data is saved in `path.data`.

    best.fit.poly <- find_best_fits(df = peaks.df.100.avg,
                                    plot.to.screen = TRUE,
                                    path.data = path.data,
                                    path.plots = path.plots)

    best.fit.poly
    #> [1] 10

Both 4 and 8 coefficients most often were under the first 4 models that were followed by an AIC-change below `5%`:

The first image shows the polynomial fits to the average force curves of species_A and species_B. The second image shows a frequency distribution of how often a polynomial function with n coefficients was followed by a function with n+1 coefficients that only brought an AIC-change below `5%`.

Given the fact that we only analyzed four species in this example, the histogram looks a bit scarce (`n = 4 species * 4 coefficients = 16`). If more species are analyzed, this may rather look like this:

#### Convert curves to polynomial models

In the last step of this vignette, we convert the average peak curve shape into polynomial models in which all have the same amount of coefficients (defined by the previously identified `best.fit.poly`):

    models <- peak_to_poly(df = peaks.df.100.avg,
                           coeff = best.fit.poly,
                           path.data = path.data,
                           show.progress = TRUE)

    # create tibble with model data
    models.df <- NULL
    for(i in 1:length(models)){
      model.df <- tibble(species = rep(names(models)[i], 100),
                         index = 1:100,
                         y = predict(models[[i]]))
      models.df <- rbind(models.df, model.df)
    }

    # plot all polynomial models
    ggplot(models.df,
           aes(x = index ,
               y = y,
               colour=species)) +
      geom_line()

## Well done!

Congratulations - you have made it to the end of this vignette. We hope this package helps you with your data. In case of any remarks or questions, please feel free to contact Peter T. Rühr at `peter.ruehr at gmail.com`.

## `forceR` workflow example

Below we have compiled a self-sufficient `R` script to run all functions of `forceR` after file loading and drift corrections. Instead of loading data from in vivo measurements, we simulate bite series with the function `simulate_bites()` and store it in `df.all`. The variable names used are the same as used in the rest of this vignette. They are also stored in the `forceR` package internally and can be restored by calling, e.g. `df.all <- forceR::df.all`.

    ## ---------------------------
    ##
    ## PURPOSE OF SCRIPT:
    ##      Simulate data and test various functions of the forceR package.
    ##
    ##
    ## AUTHOR:
    ##      Peter T. RÃ¼hr
    ##
    ## DATE CREATED:
    ##      2022-04-07
    ##
    ## Copyright (c) Peter T. RÃ¼hr, 2022
    ## Email: peter.ruehr@gmail.com
    ##
    ## ---------------------------
    ##
    ## NOTES:
    ##
    ##
    ##
    ##  -------------------------------------------------------------------------
    ##  DEPENDENCIES

    # to use viewport()
    require(grid)

    # forceR
    install.packages("C:/Users/pruehr.EVOLUTION/Documents/forceR",
                     repos = NULL,
                     type = "source")
    library(forceR)

    # various plotting functions
    require(ggplot2)

    # load tidyverse for its various conveniences
    require(tidyverse)

    ##  -------------------------------------------------------------------------
    ##  FUNCTIONS
    # get data as string for saving files
    today <- function(){
      date.string <- gsub("-", "_", substring(as.character(as.POSIXct(Sys.time())), 1, 10))
      return(date.string)
    }

    # plot linear regression and return relevant values
    plot.linear.regression <- function(x, y,
                                       logarithmic = F,
                                       x.axis.label = "x",
                                       title = NULL,
                                       x.lim = NULL, y.lim = NULL){
      if(logarithmic == "10"){x <- log10(x); y <- log10(y)}
      if(logarithmic == "e"){x <- log(x); y <- log(y)}
      lin.mod <- lm(y ~ x)
      lin.mod.sum <- summary(lin.mod)
      lin.mod.r2 <- lin.mod.sum$adj.r.squared
      lin.mod.p <- lin.mod.sum$coefficients[2,4]
      lin.mod.intercept <- lin.mod$coefficients[1]
      lin.mod.slope <- lin.mod$coefficients[2]
      lin.mod.label.r2 <- bquote(italic(R)^2 == .(format(lin.mod.r2, digits = 3)))
      if(lin.mod.p > 0.05) {lin.mod.p.ast <- lin.mod.p}
      if(lin.mod.p <= 0.05 & lin.mod.p > 0.01) {lin.mod.p.ast <- "*"}
      if(lin.mod.p <= 0.01 & lin.mod.p > 0.001) {lin.mod.p.ast <- "**"}
      if(lin.mod.p <= 0.001) {lin.mod.p.ast <- "***"}

      lin.mod.label.p <- bquote(italic(p) == .(lin.mod.p.ast))

      if(is.null(xlim)){
        x.lim = c(min(x), max(x))
      }
      if(is.null(ylim)){
        y.lim = c(min(y), max(y))
      }

      if(lin.mod.p >= 0.001) p.print <- paste0("p = ", round(lin.mod.p, 3))
      if(lin.mod.p < 0.001) p.print <- "p < 0.001"
      plot(x, y, pch = 16, xlab = paste0(x.axis.label, ": ", p.print,
                                         "; R2 = ", round(lin.mod.r2,3),
                                         "; m = ",  round(lin.mod.slope,3),
                                         "; y = ",  round(lin.mod.intercept,3)),
           xlim = x.lim, ylim = y.lim)

      abline(lin.mod, col = "red")

      if(!is.null(title)){
        title(main = title)
      }
      print(x.axis.label)
      print(paste0(p.print))
      print(paste0("R2 = ", round(lin.mod.r2,3)))
      print(paste0("m = ", round(lin.mod.slope,3)))
      print(paste0("y = ", round(lin.mod.intercept,3)))

      res <- tibble(p = lin.mod.p,
                    r2 = lin.mod.r2,
                    intercept = lin.mod.intercept,
                    slope = lin.mod.slope)
      return(res)
    }

    # add leading zeros to number and return as string
    add_eading_zeros <- function(number, length){
      require(stringr)
      str_pad(number, length, pad = "0")
    }

    ## -------------------------------------------------------------------------
    ## Simulate bite force data
    # set seed for randomization so results are reproducible
    set.seed(1)

    ## -------------------------------------------------------------------------
    ## Create classifier to store data (1 row per measurement)
    classifier <- tibble(bite.type = c(rep("sinosoidal", 5), rep("plateau", 3), rep("inter", 4)),
                                                peak.position = c("early","early","center","late","late","center","center","center","center","center","early","late"),
                                                species = NA,
                                                specimen = NA,
                                                measurement = NA,
                                                body_mass = NA,
                                                force_in = NA,
                                                length.of.bite = 4000,
                                                peak.pos = c(20, 25, 50, 65, 70, rep(NA, 7)),
                                                slope.perc.starts = c(rep (NA, 5), 10,15,20,30,40,20,50),
                                                slope.perc.ends = c(rep (NA, 5), 10,15,20,30,40,50,20),
                                                type = c(rep("sin", 5), rep("plat", 7)),
                                                no.of.bites = 7,
                                                amp = 1,
                                                lever.ratio = 1,
                                                length.of.series = 35000)

    # amend classifier
    classifier <- classifier %>%
      mutate(no = row_number())
    classifier

    # define maximum forces per bite simulation type and add additional rows to classifier
    forces <- c(1, 5, 15)
    for(i in 1:nrow(classifier)){
      classifier$force_in[i] <- forces[1]
      classifier <- classifier %>%
        add_row(classifier[i, ])
      classifier$force_in[nrow(classifier)] <- forces[1]
      for(f in 2:3){
        classifier <- classifier %>%
          add_row(classifier[i, ])
        classifier$force_in[nrow(classifier)] <- forces[f]
        classifier <- classifier %>%
          add_row(classifier[i, ])
        classifier$force_in[nrow(classifier)] <- forces[f]
      }
    }

    # arrange classifier by original bite.type and remove bite.typeing column
    classifier <- classifier %>%
      arrange(no) %>%
      select(-no)

    # create vector of species names and add to classifier
    species.names <- NULL
    for(i in 1:(nrow(classifier)/2)){
      species.names <- c(species.names,
                         rep(paste0("S", add_eading_zeros(i, 2)), 2))
    }
    classifier$species <- species.names

    classifier$specimen <- paste0("s", add_eading_zeros(1:nrow(classifier), 3))
    classifier$measurement <- paste0("m", add_eading_zeros(1:nrow(classifier), 3))

    # assign body mass according to the maximum force
    classifier$body_mass[classifier$force_in == forces[1]] <- 1
    classifier$body_mass[classifier$force_in == forces[2]] <- 6
    classifier$body_mass[classifier$force_in == forces[3]] <- 25

    # add jitter to force and body mass to replicate biological variation
    for(i in 1:nrow(classifier)){
      classifier$force_in[i] <- round(classifier$force_in[i] +
                                        ((rnorm(1, -0.2, 0.2)) * classifier$force_in[i]), 2)
      classifier$body_mass[i] <- round(classifier$body_mass[i] +
                                         ((rnorm(1, -0.2, 0.2)) * classifier$body_mass[i]), 2)
    }

    # -------------------------------------------------------------------------
    # data processing

    # get overview of input data before simulating bite series
    BFQ.regression_in <- plot.linear.regression(x = classifier$body_mass,
                                                y = classifier$force_in,
                                                logarithmic = "10",
                                                x.axis.label = "body mass")

    # jitter for variation in maximum bite force within a bite series
    # this was set to 0 when checking if the package finds the correct max. force values
    # and to 15 to increase bite shape diversity when checking if the package can tell the different
    # bite shapes apart
    max.y.jit = 0 # 0 15

    # jitter to make the bite curve more unstable
    # this was set to 0 when checking if the package finds the correct max. force values
    # and to 1 to increase bite shape diversity when checking if the package can tell the different
    # bite shapes apart
    jit = 0 # 0 1

    # create tibble with simulated time series with different
    # bite characteristics for each measurement, specimen and species

    # path.plots <- "Z:/PAPERS/PTR_Bite force METHODS/R/package_tests/plots/"
    # print(paste0("Saving plots at ", path.plots, "/", today(),"_bite_series.pdf..."))
    # pdf(paste0(path.plots, "/", today(),"_bite_series.pdf..."), onefile = TRUE, paper = "a4", height = 14)
    par(mfrow = c(3,2))
    df.all <- NULL
    for(i in 1:nrow(classifier)){
      df.curr <- simulate_bites(no.of.bites = 5,
                                length.of.bite = classifier$length.of.bite[i],
                                length.of.series = 5*classifier$length.of.bite[i] + 5*1000,
                                max.y = classifier$force_in[i],
                                max.y.jit = max.y.jit,
                                jit = jit,
                                peak.pos = classifier$peak.pos[i],
                                slope.perc.start <- classifier$slope.perc.starts[i],
                                slope.perc.end <- classifier$slope.perc.ends[i],
                                bite.type = classifier$type[i],
                                plot = TRUE)

      # add measurement number to df.curr
      df.curr <- df.curr %>%
        mutate(measurement = classifier$measurement[i])

      # add current sumulated bite series to df.all
      df.all <- rbind(df.all, df.curr)
    }
    # dev.off()

    # remove columns from classifier that were only used during bite series simulation
    classifier <- classifier %>%
      select(-c(jit, length.of.bite, peak.pos,
                slope.perc.starts, slope.perc.ends,
                type, no.of.bites))

    # reduce sampling frequency to 200 Hz
    df.all.200 <- reduce_frq(df = df.all,
                             Hz = 200,
                             measurement.col = "measurement")

    # convert y values to force and add measurement columns from classifier info (df.all)
    df.all.200.tax <- y_to_force(df = df.all.200,
                                 classifier = classifier,
                                 measurement.col = "measurement")

    # add species to df.all.200.tax
    df.all.200.tax <- df.all.200.tax %>%
      left_join(classifier %>%
                  select(species, measurement))

    # summarize force data per specimen
    df.summary.specimen <- summarize_measurements(df = df.all.200.tax,
                                                  var1 = "measurement",
                                                  var2 = "specimen")

    # add body mass from classifier:
    df.summary.specimen <- df.summary.specimen %>%
      left_join(classifier %>%
                  select(measurement, body_mass),
                by = "measurement")

    # boxplot of maximum force of all specimens
    ggplot(data = df.summary.specimen, mapping = aes(x=specimen,y=max.F.measurement)) +
      geom_jitter(color='blue',alpha=0.5, width = 0.2) +
      geom_boxplot(fill="blue",color="black",alpha=0.1) +
      # scale_y_log10() +
      labs(x='specimen', y="max(F)/specimen") +
      guides(color="none") +
      theme_minimal() +
      ggtitle("max. bite force per measurement") +
      xlab("specimen") +
      ylab("max. force (N)") +
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1, size = 5))

    # Summarize to species-wise info
    # We are not using the summarize_measurements() functions because this would ignore
    # the fact the some measurements may come from the same specimen, but we only want
    # to consider on maximum force value per specimen and not on per measurement.
    df.summary.species <- df.summary.specimen %>%
      # find max Fs of species
      group_by(species) %>%
      # calculate force values for each species
      mutate(max.F.species = max(max.F.specimen),
             mean.F.species = round(mean(max.F.specimen),6),
             sdv.max.F.species = sd(max.F.specimen)) %>%
      ungroup() %>%
      # count specimens / species
      group_by(species) %>%
      mutate(n.specimens.in.species = length(unique(specimen))) %>%
      # add body mass from classifier
      left_join(classifier %>%
                  select(measurement),
                by = c("measurement")) %>%
      # calculate mean body mass per species
      group_by(species) %>%
      mutate(body_mass.species = mean(body_mass)) %>%
      ungroup()

    # boxplot of maximum force in species
    ggplot(data = df.summary.species, mapping = aes(x=species,y=mean.F.species)) +
      geom_jitter(color='blue',alpha=0.5, width = 0.2) +
      geom_boxplot(fill="blue",color="black",alpha=0.1) +
      # scale_y_log10() +
      labs(x='species', y="mean(F)/species") +
      guides(color="none") +
      theme_minimal() +
      ggtitle("max. bite force per species") +
      xlab("species") +
      ylab("max. force (N)") +
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1, size = 5))

    # calculate and plot the regressions of known (simulation inputs) and extracted forces over body length
    # pdf(file = paste0(path.plots, today(), "_regressions.pdf"),
    #     paper = "special", width = 5, height = 12)
    par(mfrow=c(3,1))
    # Specimen-wise regression of known maximum force values over body mass
    BFQ.regression_in <- plot.linear.regression(x = classifier$body_mass,
                                                y = classifier$force_in,
                                                logarithmic = "10",
                                                x.axis.label = "body mass")

    # Specimen-wise regression of extracted maximum force data over body mass
    BFQ.regression.specimen <- plot.linear.regression(x = df.summary.specimen$body_mass,
                                                      y = df.summary.specimen$max.F.specimen,
                                                      logarithmic = "10",
                                                      x.axis.label = "body mass")

    # Species-wise regression of extracted maximum force data over body mass
    BFQ.regression.species <- plot.linear.regression(x = df.summary.species$body_mass.species,
                                                     y = df.summary.species$mean.F.species,
                                                     logarithmic = "10",
                                                     x.axis.label = "body mass")

    # dev.off()
    par(mfrow=c(1,1))

    # calculate differences between known and extracted maximum forces
    force.comp <- classifier %>%
      select(measurement, force_in) %>%
      left_join(df.summary.specimen %>%
                  select(measurement, max.F.measurement),
                by = "measurement") %>%
      mutate(diff.abs = max.F.measurement - force_in,
             diff.perc = diff.abs*100/force_in) %>%
      arrange(diff.perc)

    # violin plot of differences between known and extracted maximum for per specimen [in %]
    ggplot(data = force.comp, mapping = aes(x= 1, y=diff.perc)) +
      geom_jitter(color='blue',alpha=0.7, width = 0.1) +
      geom_violin(fill="blue",color="black",alpha=0.1) +
      # scale_y_log10() +
      labs(x='', y="diff. pred/actual [%]") +
      guides(color="none") +
      theme_minimal()

    # extract coefficients of species-wise regression
    # to calculate bite force quotient (BFQ; Wroe et al. 2005, Christiansen & Wroe 2007)
    regression.m <- BFQ.regression.species$slope
    regression.y <- BFQ.regression.species$intercept

    # calculate BFQ per species
    df.summary.species$BFQ.body_mass <- 100*(df.summary.species$mean.F.species/
                                               10^(regression.m * log10(df.summary.species$body_mass) + regression.y))

    # plot species-wise BFQ as histogram
    hist(df.summary.species$BFQ.body_mass, breaks = 25)

    # INDIVIDUAL BITE CURVE FINDING
    # find five strongest peaks per species
    peaks.df <- find_strongest_peaks(
      df = df.all.200.tax,
      no.of.peaks = 5,
      path.data = NULL,
      path.plots = NULL,
      show.progress = TRUE)

    # save plots of every peak in a PDF
    plot_peaks(df.peaks = peaks.df,
               df.data = df.all.200.tax,
               additional.msecs = 2000,
               path.plots = NULL,
               show.progress = TRUE)

    # rescale bites
    peaks.df.norm <- rescale_peaks(df.peaks = peaks.df,
                                   df.data = df.all.200.tax,
                                   plot.to.screen = FALSE,
                                   path.data = NULL,
                                   show.progress = TRUE)

    # check if rescaling worked: both following lines should print 1
    max(peaks.df.norm$t.norm)
    max(peaks.df.norm$force.norm)

    # reduce to 100 observations per bite
    peaks.df.norm.100 <- red_peaks_100(df = peaks.df.norm,
                                       plot.to.screen = TRUE,
                                       path.plots = NULL,
                                       path.data = NULL,
                                       show.progress = TRUE)

    # get average bite curve per species
    peaks.df.100.avg <- avg_peaks(df = peaks.df.norm.100,
                                  path.data = NULL)

    # find best polynomial degree to describe all average curves
    best.fit.poly <- find_best_fits(df = peaks.df.100.avg,
                                    plot.to.screen = FALSE,
                                    path.data = NULL,
                                    path.plots = NULL,
                                    show.progress = TRUE)

    # convert species-wise average curves to polynomial models
    models <- peak_to_poly(df = peaks.df.100.avg,
                           coeff = best.fit.poly,
                           path.data = NULL,
                           show.progress = TRUE)

    # convert models to PCA input data
    pca.data <- sapply(models, function(x){
      cbind(x[[1]]) # extract the common coefficients
    })

    # transpose PCA data (coefficients of polynomial models)
    pca.data <- t(pca.data)

    # perform PCA
    PCA.bite.shape <- prcomp(pca.data)
    summary(PCA.bite.shape)

    # store and principal component scores in a tibble
    PCA.res <- as_tibble(PCA.bite.shape$x[,1:3]) %>%
      mutate(species = rownames(PCA.bite.shape$x)) %>%
      left_join(classifier %>%
                  select(bite.type, peak.position, species),
                by = "species") %>%
      select(bite.type, peak.position, species, PC1, PC2) %>%
      distinct(species, .keep_all = TRUE) %>%
      dplyr::rename(peak.position = peak.position,
                    bite.type = bite.type)

    # plot PC1 against PC2
    ggplot(data = PCA.res, aes(x = PC1, y = PC2, col = bite.type)) +
      geom_point() +
      theme_minimal()

    # pdf(file = paste0(path.plots, today(), "_PCA_bite_shape.pdf"),
    #     paper = "a4r", width = 29, height = 21) # , height = 14
    ggplot(data = PCA.res, aes(x = PC1, y = PC2, col = peak.position)) +
      geom_point() +
      theme_minimal()
    # dev.off()

    # plot PC1 against PC1 with bite shapes as insets
    # pdf(file = paste0(path.plots, today(), "_PCA_bite_shape_w_curves.pdf"),
    #     paper = "a4r", width = 20, height = 21) # , height = 14
    # create main PCA plot
    main_plot <- ggplot(data = PCA.res, aes(x = PC1, y = PC2, col = peak.position)) +
      geom_point() +
      theme(axis.title.x=element_blank(),
            axis.text.x=element_blank(),
            axis.ticks.x=element_blank(),
            axis.title.y=element_blank(),
            axis.text.y=element_blank(),
            axis.ticks.y=element_blank(),
            legend.position = "none")
    print(main_plot)

    # create an and plot an inlet with the bite shape
    # of the first bite of the first
    # simulation settings of each specimen
    species_to_plot <- paste0("S", seq(1, nrow(PCA.res), 3))
    for(i in seq(1, nrow(PCA.res), 3)){
      curr.PC1 <- PCA.res$PC1[i]
      curr.PC2 <- PCA.res$PC2[i]
      curr.species <- PCA.res$species[i]
      curr.bite.data <- peaks.df.100.avg %>%
        filter(species == curr.species)

      inset_plot <- ggplot(curr.bite.data, aes(index, force.norm.100.avg)) +
        geom_line() +
        theme(axis.title.x=element_blank(),
              axis.text.x=element_blank(),
              axis.ticks.x=element_blank(),
              axis.title.y=element_blank(),
              axis.text.y=element_blank(),
              axis.ticks.y=element_blank(),
              plot.margin=grid::unit(c(0,0,0,0), "null"),
              panel.background = element_rect(fill = 'white', colour = 'white')) +
        theme(aspect.ratio=1)

      #A viewport taking up a fraction of the plot area
      vp <- viewport(width = 0.1, height = 0.1,
                     x = (curr.PC1-min(PCA.res$PC1))/diff(range(PCA.res$PC1)),
                     y =(curr.PC2-min(PCA.res$PC2))/diff(range(PCA.res$PC2)))

      print(inset_plot, vp = vp)
    }
    # dev.off()