processing.qmd

# Processing

In this section of the training workshop, we dive into the basic processing operations of soil spectroscopy data. This includes importing spectral data, making tabular and element-wise operations, visualization, resampling, preprocessing, and compression.

At the end, we will produce two files using a specific processing pipeline (`train.csv` and `test.csv`) that will be employed in the next [Machine Learning](#sec-ml) section.

It is important to mention that all files are being saved to an external working directory, i.e., we are not storing big files inside this GitHub and book repository.

You can set an external folder as your working directory and RStudio project, and copy/paste the code chunks of this session into plain R scripts.

Please, set the working directory (or create an RStudio project) in your local machine:
```{r working_directory}
my.wd <- "~/projects/local-files/soilspec_training"
# Or within an RStudio project
# my.wd <- getwd()
```

A list of all required packages for this section is provided in the following code chunk. You will see that some specific/special functions are highlighted in the text by linking back the function with the original package using the `package::function()` syntax:
```{r setup, message=FALSE, warning=FALSE}
library("tidyverse")
library("asdreader")
library("opusreader2")
library("prospectr")
library("qs")
library("moments")
library("resemble")
```

## Importing spectra

At the beginning of a project, we need to import the spectral measurements as raw binary files (like `ASD` and `OPUS` files) or other text/data file formats that are commonly used across different software, for instance, a `CSV` file.

For learning some of the common operations of this part, let's use some of the datasets shared through the [Open Soil Spectral Library (OSSL)](https://soilspectroscopy.github.io/ossl-manual/). The OSSL has data for different spectral regions: visible and near-infrared [`VisNIR`, 350-2500 nm], near-infrared [`NIR`, 1350-2550 nm], and middle-infrared (`MIR`, 4000-600 cm^-1^).

A common raw format for `VisNIR` measurements is the `.asd`. This format is used across the Malvern Panalytical instruments, like the [ASD FieldSpec](https://www.malvernpanalytical.com/en/products/product-range/asd-range/fieldspec-range) models. In R, we can import `ASD` files using the **asdreader** package. After downloading an example `.asd` file, we can use the function `asdreader::get_spectra()` by indicating the local file path. The imported spectra is a matrix with 1 row and 2551 columns, with column names as integer numbers (wavelength, nm) ranging from 350 to 2500 nm.

```{r asdreader, message=FALSE, warning=FALSE}
# Downloading an .asd file
visnir.spectra.url <- "https://github.com/soilspectroscopy/ossl-models/raw/main/sample-data/101453MD01.asd"

visnir.spectra.path <- file.path(my.wd, "file1.asd")

download.file(url = visnir.spectra.url,
              destfile = visnir.spectra.path,
              mode = "wb")

# Reading asd file
visnir.spectra <- asdreader::get_spectra(visnir.spectra.path)

# Inspecting the file
class(visnir.spectra)
dim(visnir.spectra)
visnir.spectra[1,1:5]

# Spectral range
range(as.numeric(colnames(visnir.spectra)))
```

The same operation can be done with some MIR measurements. The `OPUS` file (`.0`) is a common binary file format used across the instruments of Bruker Optics GmbH & Co. According to the original producers of **opusreader2**:

> (...) **opusreader2** is a state-of-the-art [opus] binary reader. We recommend the package as a solid foundation for your spectroscopy workflow. It is modular and has no hard dependencies apart from base R. (...) The Bruker corporation manufactures reliable instruments but there is no official documentation of the OPUS file format.

As the `OPUS` format is proprietary, [spectra-cockipit.space](https://spectral-cockpit.space/) and the open source community have cracked and reverse-engineered the binary files to be open directly in R. After downloading an example `.0` file from the OSSL project, we can use the function `opusreader2::read_opus_single()` by just indicating the local file path. The imported spectra is a list with several information (metadata and spectral data), with the spectra having column names formatted as floating numbers (wavenumbers, cm^-1^) due to the Fourier Transformation, and ranging from 599.7663 to 7498.0428 cm^-1^.

```{r opusreader2, message=FALSE, warning=FALSE}
# Downloading an .0 file
mir.spectra.url <- "https://github.com/soilspectroscopy/ossl-models/raw/main/sample-data/235157XS01.0"

mir.spectra.path <- file.path(my.wd, "file2.0")

download.file(url = mir.spectra.url,
              destfile = mir.spectra.path,
              mode = "wb")

# Reading asd file
mir.spectra <- opusreader2::read_opus_single(dsn = mir.spectra.path)

# Inspecting the file
class(mir.spectra)
names(mir.spectra)

# Spectra is stored in file$ab$data
class(mir.spectra$ab$data)
dim(mir.spectra$ab$data)

# Spectral range
range(as.numeric(colnames(mir.spectra$ab$data)))
```

In many cases, instead of importing the raw binary files, we can directly import a `CSV` exported from those spectral instruments and their accompanying software. For example, using an example `CSV` file exported from a Neospectra device that is also available through the OSSL project, the imported spectra is a table with scans in the rows, the first column as ID, and the spectra having column names formatted as floating numbers (wavelength, nm) due to the Fourier Transformation, ranging in this case, from 1350 to 2550 nm.

```{r read_csv, message=FALSE, warning=FALSE}
# Downloading an csv output from Neospectra
nir.spectra.url <- "https://github.com/soilspectroscopy/ossl-models/raw/main/sample-data/sample_neospectra_data.csv"

nir.spectra.path <- file.path(my.wd, "file3.csv")

download.file(url = nir.spectra.url,
              destfile = nir.spectra.path,
              mode = "wb")

# Reading csv file
nir.spectra <- readr::read_csv(nir.spectra.path)

# Inspecting the file
class(nir.spectra)
nir.spectra[1:5,1:5]

# Spectral range after removing first column
range(as.numeric(colnames(nir.spectra[,-1])))
```

## Tabular operations

When we upload spectra into R, we usually need to perform some operations across rows, columns, or in an element-wise mode. For example, the measurement unit may differ across instruments and spectral ranges, so when we integrate different datasets, we need to harmonize the original datasets to a specific format.

Using the Neospectra example dataset imported in the previous subsection, we can see that the scale of the measurements is provided in percent units. Rather than having the data in 0-100% percent units, we can transform it to reflectance factor units in the 0-1 interval and keep 5 decimal places of precision.

```{r nir_reflectance}
# Original data
nir.spectra[1:5,1:5]

# Spectra column names
spectra.column.names <- nir.spectra %>%
  select(-sample_id) %>%
  names()

# Transforming reflectance (%) to reflectance factor (decimal, 0-1)
# Also, rounding to 5 decimal places of precision
nir.spectra.rf <- nir.spectra %>%
  mutate(across(all_of(spectra.column.names), ~round(.x/100, 5)))

nir.spectra.rf[1:5,1:5]
```

We can also use the same element-wise operations to transform data in absorbance (`A`, in log10 units) or reflectance (`R`, reflectance factor 0-1). This is not run here, but the following equations and R code can be used:

- Absorbance from reflectance: $$A=\log_{10}\left(\frac{1}{R}\right)$$ or `mutate(across(all_of(spectra.column.names), ~round(log10(1/.x), 5)))`.  
- Reflectance from absorbance: $$R=\frac{1}{10^{A}}$$ or `mutate(across(all_of(spectra.column.names), ~round(1/(10^.x), 5)))`.

## Visualization

Another important operation is to to be able to visualize the spectra. For this task, we can use the **ggplot2** package after pivoting the [wide] table to a long version that stores the data in two new columns: wavelength (`x` variable) and reflectance (`y` variable).

```{r visualization}
## Pivot to long format
nir.spectra.rf.long <- nir.spectra.rf %>%
  pivot_longer(all_of(spectra.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance))

head(nir.spectra.rf.long)

## Visualization
ggplot(data = nir.spectra.rf.long) +
  geom_line(aes(x = wavelength, y = reflectance,
                group = sample_id),
            alpha = 0.5, linewidth = 0.5) +
  theme_light()
```

## Resampling spectra

From the previous table views, we saw the spectral column headers are being represented by an uneven interval with floating numbers. We can resample or harmonize the spectra (using the **prospectr** R package) to a defined range with an even interval (e.g. 2 nm) using spline interpolation. For this, we need to use the spectra stored as a wide table

::: {.callout-note}
Linear and spline interpolation do not work outside of the original range. If you have missing data out of the original range, you must use a different approach like imputation to fill the gaps.
:::

```{r interpolation_columns}
# Old columns, reversed, and as numeric
old.wavelength <- as.numeric(rev(spectra.column.names))
head(old.wavelength)

# New columns, increasing order, spaced 2 nm
new.wavelength <- seq(1350, 2550, by = 2)
head(new.wavelength)
```

::: {.callout-tip}
## Tip
Dot (`.`) is used as a placeholder for the previous output of the pipe. We can create quick internal pipes with brackets `{}`.
:::

```{r interpolation}
# Selecting old spectra in increasing order
# Parse to matrix (input of prospectr::resample)
# Resample
# Parse to tibble
# Bind to original sample ids
nir.spectra.rf.int <- nir.spectra.rf %>%
  select(all_of(rev(spectra.column.names))) %>%
  as.matrix() %>%
  prospectr::resample(X = .,
                      wav = old.wavelength,
                      new.wav = new.wavelength,
                      interpol = "spline") %>%
  as_tibble() %>%
  bind_cols({nir.spectra.rf %>%
      select(sample_id)}, .)

nir.spectra.rf.int[1:5,1:5]
```

The harmonization/resample results in the same spectral patterns, but now the data is well formatted.

::: {.callout-tip}
## Tip
We can pipe together pivot, mutate, and ggplot operations.
:::

```{r visualization_resample}
new.spectra.column.names <- as.character(new.wavelength)

nir.spectra.rf.int %>%
  pivot_longer(all_of(new.spectra.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance, group = sample_id),
            alpha = 0.5, linewidth = 0.5) +
  theme_light()
```

## Preprocessing

**prospectr** is a very useful package for signal processing and chemometrics as it contains various utilities for working with spectral data. As stated in the package vignette:

> The aim of spectral preprocessing is to enhance signal quality before modeling as well as to remove physical information from the spectra. Applying a pre-treatment can increase the repeatability/reproducibility of the method, model robustness and accuracy, although there are no guarantees this will actually work.

There are several algorithms available through **prospectr** and elsewhere, e.g. listed in [Table 1](https://cran.r-project.org/web/packages/prospectr/vignettes/prospectr.html#signal-processing) of its vignette. However, in this section we are going to showcase only a few that are more common: [Savitzky–Golay (SG) smoothing and derivatives](https://cran.r-project.org/web/packages/prospectr/vignettes/prospectr.html#savitzky-golay-filtering), and [Standard Normal Variate (SNV)](https://cran.r-project.org/web/packages/prospectr/vignettes/prospectr.html#scatter-and-baseline-corrections).

[Savitzky–Golay](https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter) is an algorithm that fits a moving local polynomial regression that can smooth and/or derive the spectra (using the fitted polynomial) to enhance the signal quality and absorption features. The parameters are the polynomial order (`p`), the half-window size to sample and fit the spectra (`w`), the derivative order (`m`, where `m = 0` is used for smoothing and `m > 1` are the respective derivatives), and the spacing interval (`delta.wav`). When we use the SG algorithm, the edges of the spectral range are reduced by one half-window size minus the center.

Except for `delta.wav`, we can fine tune all of these parameters to find the best preprocessing combination, although several studies have already explored and we can adopt what has been recommended [@Dotto2018; @Seybold2019; @Barra2021]. In this example, we are going to use the first-derivative of a second-order polynomial regression with a half-window size of 11 nm.

```{r sg}
# Select spectra columns
# Parse to matrix (input of prospectr::savitzkyGolay)
# Apply preprocessing
# Parse to tibble
# Bind the spectra to id column
nir.spectra.sg <- nir.spectra.rf.int %>%
  select(all_of(new.spectra.column.names)) %>%
  as.matrix() %>%
  savitzkyGolay(X = ., p = 2, w = 11, m = 1, delta.wav = 2) %>%
  as_tibble() %>%
  bind_cols({nir.spectra.rf %>%
      select(sample_id)}, .)

nir.spectra.sg[1:5,1:5]
```
Moving-window first-derivatives are great because they preserve the orientation of the absorption features and enhance only the important regions of the spectra. For more complex spectra, however, this may hamper the interpretation.

::: {.callout-tip}
## Tip
**dplyr** allows us to use column selectors like `first()`, `everything()`, `all_of` etc. As we now have a shorter spectral range, rather than using `all_of(new.spectra.column.names)` we can replace it with `any_of(new.spectra.column.names)` to select any of the available columns.
:::

```{r sg_visualization}
nir.spectra.sg %>%
  pivot_longer(any_of(new.spectra.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance, group = sample_id),
            alpha = 0.5, linewidth = 0.5) +
  theme_light()
```
Another very common preprocessing is the Standard Normal Variate (SNV). SNV is a normalization algorithm (centers to mean 0 and rescales to 1 standard deviation) that works across the spectrum (row-wise). This preprocessing changes both the range of values and the amplitude of the curves and is intended to correct the scattering of light.

SNV was originally proposed to deal with multiplicative effects of particle size, light scatter, and multicollinearity issues in diffuse reflectance spectroscopy [@Barnes1989]. Although the first derivative has been routinely used in soil spectroscopy studies, SNV is in many cases preferred because it does not affect the interpretation of the spectral features.

```{r snv}
# Select spectra columns
# Parse to matrix (input of prospectr::savitzkyGolay)
# Apply preprocessing
# Parse to tibble
# Bind the spectra to id column
nir.spectra.snv <- nir.spectra.rf.int %>%
  select(all_of(new.spectra.column.names)) %>%
  as.matrix() %>%
  prospectr::standardNormalVariate(X = .) %>%
  as_tibble() %>%
  bind_cols({nir.spectra.rf %>%
      select(sample_id)}, .)

nir.spectra.sg[1:5,1:5]
```

```{r snv_visualization}
nir.spectra.snv %>%
  pivot_longer(any_of(new.spectra.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance, group = sample_id),
            alpha = 0.5, linewidth = 0.5) +
  theme_light()
```

## Preparing machine learning dataset {#subsec-preparing}

For the [Machine Learning](#sec-ml) task and compression/projection, we will use the NIR Neospectra Handheld database shared as part of the OSSL that contains more than 2000 unique soil samples scanned in replicates with different devices (identified by serial number).

We will prepare a case-study dataset using several operations that include reading, filtering, preprocessing with **SNV**, and compressing the spectra with **PCA** to predict Soil Organic Carbon (SOC).

```{r neospectra_read, eval=TRUE}
## Internet configuration for downloading big datasets
options(timeout = 10000)

## Reading serialized files
neospectra.soil <- qread_url("https://storage.googleapis.com/soilspec4gg-public/neospectra_soillab_v1.2.qs")
dim(neospectra.soil)
neospectra.soil[1:5,1:5]

neospectra.site <- qread_url("https://storage.googleapis.com/soilspec4gg-public/neospectra_soilsite_v1.2.qs")
dim(neospectra.site)
neospectra.site[1:5,1:5]

neospectra.nir <- qread_url("https://storage.googleapis.com/soilspec4gg-public/neospectra_nir_v1.2.qs")
dim(neospectra.nir)
neospectra.nir[1:5,1:5]
```

Once we download the database, we can see that we have 2,106 unique soil samples and more than 8,000 measurements. We need to explore it a bit further and calculate the average spectra across the different devices, so each soil sample is linked to a unique measurement. We also select only the relevant information for using in further processing:

- Soil sample id: `id.sample_local_c`.  
- Origin country: `location.country_iso.3166_txt`.  
- Organic carbon content (%): `oc_usda.c729_w.pct.`  
- Average NIR spectra: columns starting with `scan_nir`.  

```{r neospectra_filter, eval=TRUE}
# Column names
neospectra.site %>%
  names()

neospectra.soil %>%
  names()

neospectra.nir %>%
  select(1:20) %>%
  names()

# How many distinct id.sample_local_c in the NIR file?
neospectra.nir %>%
  distinct(id.sample_local_c) %>%
  nrow()

# How many distinct countries and samples from?
# 2016 samples from USA, other from African countries
# We can use this columns to split and block datasets
neospectra.site %>%
  count(location.country_iso.3166_txt)

# Selecting relevant site data
neospectra.site <- neospectra.site %>%
  select(id.sample_local_c, location.country_iso.3166_txt)

# Selecting relevant soil data
neospectra.soil <- neospectra.soil %>%
  select(id.sample_local_c, oc_usda.c729_w.pct)

# Selecting relevant NIR data and taking average across repeats
neospectra.nir <- neospectra.nir %>%
  select(id.sample_local_c, starts_with("scan_nir")) %>%
  group_by(id.sample_local_c) %>%
  summarise_all(mean, .group = "drop")

# Renaming spectral columns headers to numeric integers
neospectra.nir %>%
  select(starts_with("scan_nir")) %>%
  names() %>%
  head()

old.names <- neospectra.nir %>%
  select(starts_with("scan_nir")) %>%
  names()

new.names <- gsub("scan_nir.|_ref", "", old.names)

neospectra.nir <- neospectra.nir %>%
  rename_with(~new.names, all_of(old.names))

spectral.column.names <- new.names

neospectra.nir[1:5,1:5]
```

**Spectral visualization.**

```{r neospectra_view, eval=TRUE}
neospectra.nir %>%
  pivot_longer(any_of(spectral.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance, group = id.sample_local_c),
            alpha = 0.25, linewidth = 0.25) +
  theme_light()
```

**Preprocessing with Standard Normal Variate (SNV).**

```{r neospectra_snv, eval=TRUE}
neospectra.nir.snv <- neospectra.nir %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix() %>%
  prospectr::standardNormalVariate(X = .) %>%
  as_tibble() %>%
  bind_cols({neospectra.nir %>%
      select(id.sample_local_c)}, .)

neospectra.nir.snv %>%
  pivot_longer(any_of(spectral.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance, group = id.sample_local_c),
            alpha = 0.25, linewidth = 0.25) +
  theme_light()
```

**Joining data and spliting for train/test**

```{r neospectra_join, eval=TRUE}
# Joining data
neospectra <- left_join(neospectra.site,
                        neospectra.soil,
                        by = "id.sample_local_c") %>%
  left_join(., neospectra.nir.snv, by = "id.sample_local_c")

# Filtering out samples without SOC values
neospectra <- neospectra %>%
  filter(!is.na(oc_usda.c729_w.pct))

neospectra[1:5,1:5]

# Preparing train and test split
neospectra.train <- neospectra %>%
  filter(location.country_iso.3166_txt == "USA")

neospectra.test <- neospectra %>%
  filter(location.country_iso.3166_txt != "USA")
```

## Spectral compression and projection

Another interesting operation that is routinely employed in soil spectroscopy is the compression and projection of the spectra into new coordinate systems.

![Source: carpentries-incubator.github.io](https://carpentries-incubator.github.io/high-dimensional-stats-r/fig/pendulum.gif)

**Partial Least Squares (PLS)** and **Principal Component Analysis (PCA)** are popular orthogonal projection methods that helps to build regression models or explore the associations between spectra and variables of interest. Within this scope, orthogonalization means the production of uncorrelated features (predictors) such as the correlation between any pair of the new variables is 0 [@chang2001; @deSantana2023].

**PLS** is preferred if the goal is to find orthogonal features of the spectra that are more predictive of a soil property. On the other hand, if the goal is to find orthogonal features that maximally explain the variability of the spectra (irrespective of the variable of interest), then **PCA** can be adopted.

Both methods have the advantage of producing uncorrelated features (control colinearity) and reduce the dimensionality of the original data. A few orthogonal features obtained from the spectra can be employed to replace the the original spectra matrix when training a model, with the downside of sacrificing some of the original spectral information

This opens doors for the optimization between the balance of retained representation and compression magnitude. Similarly, each spectral range must have its own compression and projection model, thus the orthogonal features of each spectral range account for different amounts of the total original variance.

The values of the new uncorrelated features are usually named `scores`.

We can use the package **resemble** to perform [**Orthogonal Projection**](https://cran.r-project.org/web/packages/resemble/vignettes/resemble.html#5_Dimensionality_reduction) of high dimensional data, i.e., spectra using either **PCA** and **PLS** with different optimization criteria.

```{r compression, eval=TRUE}
train.pca <- neospectra.train %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix() %>%
  resemble::ortho_projection(Xr = .,
                             Xu = NULL,
                             Yr = NULL,
                             method = "pca",
                             pc_selection = list(method = "cumvar",
                                                 value = 0.99),
                             center = TRUE, scale = TRUE)

names(train.pca)

# How many components where retained?
plot(train.pca, col = "#D42B08CC")

# We can get the scores and exp. variance back and make a plot
train.pca.scores <- train.pca$scores %>%
  as_tibble() %>%
  bind_cols({neospectra.train %>%
      select(id.sample_local_c,
             location.country_iso.3166_txt,
             oc_usda.c729_w.pct)}, .)

train.pca$variance

train.pca.expvar <- round(train.pca$variance$x_var["explained_var",]*100, 2)

p.pca <- ggplot(train.pca.scores) +
  geom_point(aes(x = pc_1, y = pc_2, color = log1p(oc_usda.c729_w.pct)),
             alpha = 0.5, size = 0.5) +
  scale_colour_gradient(low = "gold",
                        high = "darkred",) +
  labs(title = "Neospectra PCA compression",
       x = paste0("PC1 (", train.pca.expvar[1], "%)"),
       y = paste0("PC2 (", train.pca.expvar[2], "%)")) +
  theme_light() +
  theme(legend.position = "bottom"); p.pca

# Predicting the test samples and plotting
test.pca.scores <- neospectra.test %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix() %>%
  predict(train.pca, newdata = .)
  
test.pca.scores <- neospectra.test %>%
    select(id.sample_local_c,
           location.country_iso.3166_txt,
           oc_usda.c729_w.pct) %>%
  bind_cols(test.pca.scores)

# Adding test samples to PCA plot
p.pca.test <- p.pca +
  geom_point(data = test.pca.scores,
             aes(x = pc_1, y = pc_2),
             size = 0.75) +
  labs(subtitle = "Black dots represent testing points"); p.pca.test
```

**Saving files**

```{r save, eval=TRUE}
# Saving the PCA plot
ggsave(file.path(my.wd, "plot_pca_train_test.png"), p.pca.test,
       dpi = 300, width = 4, height = 3, units = "in", scale = 1.5)

# Saving the train set
write_csv(train.pca.scores, file.path(my.wd, "train.csv"))

# Saving the test set
write_csv(test.pca.scores, file.path(my.wd, "test.csv"))
```

## Representation flag

We can leverage the PCA and any other compression algorithms to identify unrepresentative spectra when predicting unknown samples.

As only the first most important components are used for calibration, the remaining farther components can be used to test if a new sample is underrepresented in respect of the feature space of the calibration set. This underrepresentation flag might happen due to unique spectral features that are still present in farther components.

In chemometrics, this is usually named as Control Chart. A good overview of this task can be found in @deSantana2023.

In this task, each new spectrum to be predicted is back-transformed using only the first retained PCs of the model calibration. With the back-transformed version, the residual (difference between the original and back-transformed spectra) is employed to calculate the Q-statistics. The result is therefore compared against a critical value estimated across the calibration setting a predefined confidence interval (e.g., 99%).

```{r backtransform, message=FALSE, error=FALSE}
# Getting the necessary parameters for backtransformation
train.loadings <- train.pca$X_loadings
train.center <- train.pca$center
train.scale <- train.pca$scale

# Getting scores
train.scores <- predict(train.pca,
                       newdata = {
                         neospectra.train %>%
                           select(all_of(spectral.column.names)) %>%
                           as.matrix()})

test.scores <- predict(train.pca,
                       newdata = {
                         neospectra.test %>%
                           select(all_of(spectral.column.names)) %>%
                           as.matrix()})

# Backtransforming
train.bt <- train.scores %*% train.loadings
test.bt <- test.scores %*% train.loadings

# Rescaling
train.bt <- sweep(x = train.bt,
                  MARGIN = 2, FUN = "*", STATS = train.scale)

test.bt <- sweep(x = test.bt,
                 MARGIN = 2, FUN = "*", STATS = train.scale)

# Recentering
train.bt <- sweep(x = train.bt,
                  MARGIN = 2, FUN = "+", STATS = train.center)

test.bt <- sweep(x = test.bt,
                 MARGIN = 2, FUN = "+", STATS = train.center)

# Binding id columns
train.bt <- as_tibble(train.bt) %>%
  bind_cols({neospectra.train %>%
      select(-all_of(spectral.column.names))}, .)

test.bt <- as_tibble(test.bt) %>%
  bind_cols({neospectra.test %>%
      select(-all_of(spectral.column.names))}, .)

# Quick visualization
selected.id <- "65140"

plot.data <- bind_rows(
  {neospectra.test %>%
      filter(id.sample_local_c == selected.id) %>%
      mutate(source = "original", .before = 1)},
  {test.bt %>%
      filter(id.sample_local_c == selected.id) %>%
      mutate(source = "compressed", .before = 1)})

plot.data %>%
  pivot_longer(any_of(spectral.column.names),
               names_to = "wavelength",
               values_to = "reflectance") %>%
  mutate(wavelength = as.numeric(wavelength),
         reflectance = as.numeric(reflectance)) %>%
  ggplot(data = .) +
  geom_line(aes(x = wavelength, y = reflectance,
                group = source, color = source),
            alpha = 0.5, linewidth = 1) +
  theme_light() + theme(legend.position = "bottom")
```

We can see there are some subtle mismatch between the original and back-transformed spectra, and depending on the importance of those underrepresented features in the calibration models, the new spectra might not be properly compressed and yield unsatisfactory prediction. We use this approach to flag potential outliers or underrepresented samples before model prediction.

With the original and backtransformed spectra, we can subtract each other and calculate the Q-statistics (sum of squared differences)
```{r q_stats, message=FALSE, error=FALSE}
neospectra.train.spectra <- neospectra.train %>%
  arrange(id.sample_local_c) %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix()

train.bt.spectra <- train.bt %>%
  arrange(id.sample_local_c) %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix()

neospectra.test.spectra <- neospectra.test %>%
  arrange(id.sample_local_c) %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix()

test.bt.spectra <- test.bt %>%
  arrange(id.sample_local_c) %>%
  select(all_of(spectral.column.names)) %>%
  as.matrix()

# Q stats - sum of squared differences - only for test
test.q.stats <- apply((neospectra.test.spectra-test.bt.spectra)^2,
                       MARGIN = 1, sum)

# Binding q.stats to compressed data
test.pca.scores <- test.pca.scores %>%
  mutate(q_stats = test.q.stats, .before = pc_1)
```

The critical value (1% significance level) comes from the training set (used for calibration) using the equation provided in @Jackson1979, @Laursen2011, and @deSantana2023.

```{r q_critical, message=FALSE, error=FALSE}
E <- cov(neospectra.train.spectra-train.bt.spectra)
teta1 <- sum(diag(E))^1
teta2 <- sum(diag(E))^2
teta3 <- sum(diag(E))^3
h0 <- 1-((2*teta1*teta3)/(3*teta2^2))
Ca <- 2.57 # 1% significance level
Qa <- teta1*(1-(teta2*h0*((1-h0)/teta1^2))+((sqrt(Ca*(2*teta2*h0^2)))/teta1))^(1/h0)
Qa
```

Now we can compare the test samples with this critical value and flag them in the PCA plot.

```{r q_flag, message=FALSE, error=FALSE}
test.pca.scores <- test.pca.scores %>%
  mutate(represented = ifelse(q_stats <= Qa, TRUE, FALSE), .after = q_stats)

p.pca +
  geom_point(data = test.pca.scores,
             aes(x = pc_1, y = pc_2, fill = represented),
             shape = 21, size = 1.5) +
  labs(fill = "Represented?")
```