Using R for Photobiology - Functional analysis of spectra with ‘photobiology’ to ‘fda.usc’

Abstract

Explanations and example R code for applying functional data analysis to spectral data stored in the classes defined in package photobiology and plotting the obtained deepest curve with ggplot2 and ggspectra and aplying ANOVA-like methods to compare groups of spectra. The excahnge of the spectral data between packages is done with functions from package photobiologyInOut.

Tip

In this page code chunks are “folded” so as to decrease the clutter when searching for examples. Above each plot you will find a small triangle followed by “Code”. Clicking on the triangle “unfolds” the code chunk making visible the R code used to produce the text or graphic output.

Except for the loading of packages shown in section Preliminaries code examples continue through out each section.

All “words” defined in base R or in extension packages are linked to the corresponding HTML-rendered help pages.

The code in the chunks can be copied by clicking on the top right corner, where an icon appears when the mouse cursor hovers over the code listing.

Introduction

I have started measuring time series of spectra, and one approach to summarizing them and for removing outlier curves is to use functional data analysis (FDA). Sets of spectra can be also compared using an asymptotic method reminiscent of one-way ANOVA. Package ‘fda.usc’ and in particular objects of class fdata provide a unified entry point to many different methods for FDA.

Package ‘photobiologyInOut’ (>= 0.4.27) exports functions spct2fdata()and fdata2spct() making it easy to convert objects of classes source_spct, response_spct, filter_spct and reflector_spct into objects of class fdata and back.

In this page I am collecting examples of code that I have found useful. The examples use source_spct objects, but the code can be adapted to other types of spectral data.

In this first version of the page I use package ‘fda.usc’ for all examples. Newer approaches are available for FANOVA, such as those implemented in package ‘fdANOVA’.

Preliminaries

Code

library(photobiology)
library(photobiologyInOut)
library(ggspectra)
library(magrittr)
library(rlang)

Warning: package 'rlang' was built under R version 4.3.2

Code

library(fda.usc)

Warning: package 'RCurl' was built under R version 4.3.2

Code

library(lubridate)

theme_set(theme_bw())

# energy_as_default()
photon_as_default()

Deepest curve

The word “deepest” refers to the curve that is farthest from the extreme ones, based on some criterion. The examples below are for median and mean.

Examples in this chapter are for a time series of 45 spectra, measured at a frequency of one per minute close to noon on a summer day under broken clouds.

Code

load("sun.cosine.spct.Rda")

'rOmniDriver': Communication with spectrometers is disabled because environment variable OOI_HOME is not set to a path.
To acquire new data, please, make sure that the OmniDriver driver is installed and that the environment variable OOI_HOME is set to its path.

Code

sun.cosine.spct <- clean(sun.cosine.spct)
when_measured(sun.cosine.spct)[c(1, 45)]

$time.01
[1] "2023-05-29 09:24:00 UTC"

$time.45
[1] "2023-05-29 10:08:01 UTC"

Code

how_measured(sun.cosine.spct)

[1] "Acquired with MayaPro2000 (MAYP11278), R (4.3.0), 'ooacquire' (0.3.4.9000) in mode \"series-attr\",\n 'rOmniDriver' (0.1.19) and OmniDriver (2.71.0)."

As array spectrometers have a variable wavelength step between pixels, we re-express the spectra at 1 nm resolution. We also trim the shortest and longest wavelengths as these spectral irradiance values are noisy.

Code

subset2mspct(sun.cosine.spct) %>%
  clean() %>%
  trim_wl(c(280, 1000)) %>%
  interpolate_mspct(length.out = 721) -> sun.cosine.mspct

Highlight the “median” spectrum

Code

autoplot(sun.cosine.mspct, 
         annotations = c("-", "peaks")) +
  geom_line(data = . %>% spct2fdata() %>% func.med.FM() %>% fdata2spct(),
            colour = "red", linewidth = 1, alpha = 0.67)

Highlight the “mean” spectrum

Code

autoplot(sun.cosine.mspct, 
         annotations = c("-", "peaks")) +
  geom_line(data = . %>% spct2fdata() %>% func.mean() %>% fdata2spct(),
            colour = "blue", linewidth = 1, alpha = 0.67)

Highlight the “trimmed mean” spectrum

Code

autoplot(sun.cosine.mspct, 
         annotations = c("-", "peaks")) +
  geom_line(data = . %>% spct2fdata() %>% func.trim.FM(trim = 0.25) %>% fdata2spct(),
            colour = "green", linewidth = 1, alpha = 0.67)

Functional one-way ANOVA

We compare measurements done in parallel with two spectrometers with different entrance optics. Ten spectra were acquired once every 40 s inside a mixed forest next to a forest edge. A flat diffuser with cosine response and a hemispherical (or dome-shaped) diffuser with increased response at low angles of incidence are compared.

Code

load("cos.A2.2m.series.spct.Rda")
load("dome.A2.2m.series.spct.Rda")

As the two spectrometers even with identical configuration have slightly different wavelength calibrations for the individual pixels in the array, we need to re-express the spectra on identical wavelengths. We also trim the shortest and longest wavelengths as these spectral irradiance values are noisier.

Code

when_measured(cos.A2.2m.series.spct)[c(1, 10)]

$time.01
[1] "2023-07-04 15:23:40 UTC"

$time.10
[1] "2023-07-04 15:29:40 UTC"

Code

cat(how_measured(cos.A2.2m.series.spct))

Acquired with MayaPro2000 (MAYP112785), with a cosine diffuser)
R (4.3.0), 'ooacquire' (0.4.1) in mode "series", 'rOmniDriver' (0.1.19.9000) and OmniDriver (2.72).

Code

shade.cosine.mspct <- subset2mspct(cos.A2.2m.series.spct) %>%
  clean() %>%
  trim_wl(c(280, 1000)) %>%
  interpolate_mspct(length.out = 721)

Code

when_measured(dome.A2.2m.series.spct)[c(1, 10)]

$time.01
[1] "2023-07-04 15:23:40 UTC"

$time.10
[1] "2023-07-04 15:29:40 UTC"

Code

cat(how_measured(dome.A2.2m.series.spct))

Acquired with MayaPro2000 (MAYP11278), with a hemispherical diffuser)
R (4.3.1), 'ooacquire' (0.4.1) in mode "series", 'rOmniDriver' (0.1.19.9000) and OmniDriver (2.72).

Code

shade.dome.mspct <- subset2mspct(dome.A2.2m.series.spct) %>%
  clean() %>%
  trim_wl(c(280, 1000)) %>%
  interpolate_mspct(length.out = 721)

We concatenate the collections of spectra.

Code

shade.mspct <- c(cosine = shade.cosine.mspct,
                 dome = shade.dome.mspct)
names(shade.mspct)

 [1] "cosine.time.01" "cosine.time.02" "cosine.time.03" "cosine.time.04"
 [5] "cosine.time.05" "cosine.time.06" "cosine.time.07" "cosine.time.08"
 [9] "cosine.time.09" "cosine.time.10" "dome.time.01"   "dome.time.02"  
[13] "dome.time.03"   "dome.time.04"   "dome.time.05"   "dome.time.06"  
[17] "dome.time.07"   "dome.time.08"   "dome.time.09"   "dome.time.10"

A plot of all 20 spectra.

Code

autoplot(shade.mspct) +
  aes(linetype = ifelse(grepl("cosine", spct.idx), "cosine", "dome")) +
  labs(linetype = "Entrance\noptics")

We test in a one-way ANOVA if the spectra differ, and they do differ significantly.

Code

shade.fda <- mspct2fdata(shade.mspct)
shade.anova <-
  fanova.onefactor(shade.fda,
                   group = factor(rep(c("cosine", "dome"), each = 10L)),
                   plot = TRUE)

Code

shade.anova$pvalue

[1] 0

To test if the shape of the spectra differs, we first scale them to equal area under the curve, and repeat the ANOVA analysis.

Code

shade_scaled.fda <- 
  mspct2fdata(fscale(shade.mspct, f = irrad))
shade_scaled.anova <-
  fanova.onefactor(shade_scaled.fda,
                   group = factor(rep(c("cosine", "dome"), each = 10L)),
                   plot = TRUE)

Code

shade_scaled.anova$pvalue

[1] 0

References

Code

citation('fda.usc')

To cite fda.usc in publications use:

  Manuel Febrero-Bande, Manuel Oviedo de la Fuente (2012). Statistical
  Computing in Functional Data Analysis: The R Package fda.usc. Journal
  of Statistical Software, 51(4), 1-28. URL
  https://www.jstatsoft.org/v51/i04/.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {Statistical Computing in Functional Data Analysis: The {R} Package {fda.usc}},
    author = {Manuel Febrero-Bande and Manuel {Oviedo de la Fuente}},
    journal = {Journal of Statistical Software},
    year = {2012},
    volume = {51},
    number = {4},
    pages = {1--28},
    url = {https://www.jstatsoft.org/v51/i04/},
  }

Code

citation('photobiology')

To cite package ‘photobiology’ in publications use:

  Aphalo, Pedro J. (2015) The r4photobiology suite. UV4Plants Bulletin,
  2015:1, 21-29. DOI:10.19232/uv4pb.2015.1.14

A BibTeX entry for LaTeX users is

  @Article{,
    author = {Pedro J. Aphalo},
    title = {The r4photobiology suite},
    journal = {UV4Plants Bulletin},
    volume = {2015},
    number = {1},
    pages = {21-29},
    year = {2015},
    doi = {10.19232/uv4pb.2015.1.14},
  }

Reuse

CC BY-SA 4.0

--- title: "Functional analysis of spectra with 'photobiology' to 'fda.usc'" subtitle: "From 'photobiology' to 'fda.usc' and back." author: "Pedro J. Aphalo" date: "2023-06-02" date-modified: "2023-07-20" keywords: [photobiology pkg, photobiologyInOut pkg, fda.usc pkg, ggspectra pkg, functional data analysis] categories: [Photobiology examples] abstract: | Explanations and example R code for applying functional data analysis to spectral data stored in the classes defined in package [photobiology](https://docs.r4photobiology.info/photobiology/) and plotting the obtained deepest curve with [ggplot2](https://ggplot2.tidyverse.org/) and [ggspectra](https://docs.r4photobiology.info/ggspectra/) and aplying ANOVA-like methods to compare groups of spectra. The excahnge of the spectral data between packages is done with functions from package [photobiologyInOut](https://docs.r4photobiology.info/photobiologyInOut/). editor: markdown: wrap: 72 code-fold: true format: html: code-link: true code-tools: true --- ::: callout-tip In this page code chunks are "folded" so as to decrease the clutter when searching for examples. Above each plot you will find a small triangle followed by "Code". Clicking on the triangle "unfolds" the code chunk making visible the R code used to produce the text or graphic output. Except for the loading of packages shown in section **Preliminaries** code examples continue through out each section. All "words" defined in base R or in extension packages are linked to the corresponding HTML-rendered help pages. The code in the chunks can be copied by clicking on the top right corner, where an icon appears when the mouse cursor hovers over the code listing. ::: ## Introduction I have started measuring time series of spectra, and one approach to summarizing them and for removing outlier curves is to use _functional data analysis_ (FDA). Sets of spectra can be also compared using an asymptotic method reminiscent of one-way ANOVA. Package 'fda.usc' and in particular objects of class `fdata` provide a unified entry point to many different methods for FDA. Package 'photobiologyInOut' (>= 0.4.27) exports functions `spct2fdata()`and `fdata2spct()` making it easy to convert objects of classes `source_spct`, `response_spct`, `filter_spct` and `reflector_spct` into objects of class `fdata` and back. In this page I am collecting examples of code that I have found useful. The examples use `source_spct` objects, but the code can be adapted to other types of spectral data. In this first version of the page I use package 'fda.usc' for all examples. Newer approaches are available for FANOVA, such as those implemented in package 'fdANOVA'. ## Preliminaries ```{r, message=FALSE} library(photobiology) library(photobiologyInOut) library(ggspectra) library(magrittr) library(rlang) library(fda.usc) library(lubridate) theme_set(theme_bw()) # energy_as_default() photon_as_default() ``` ## Deepest curve The word "deepest" refers to the curve that is farthest from the extreme ones, based on some criterion. The examples below are for median and mean. Examples in this chapter are for a time series of 45 spectra, measured at a frequency of one per minute close to noon on a summer day under broken clouds. ```{r} load("sun.cosine.spct.Rda") sun.cosine.spct <- clean(sun.cosine.spct) when_measured(sun.cosine.spct)[c(1, 45)] how_measured(sun.cosine.spct) ``` As array spectrometers have a variable wavelength step between pixels, we re-express the spectra at 1 nm resolution. We also trim the shortest and longest wavelengths as these spectral irradiance values are noisy. ```{r} subset2mspct(sun.cosine.spct) %>% clean() %>% trim_wl(c(280, 1000)) %>% interpolate_mspct(length.out = 721) -> sun.cosine.mspct ``` ### Highlight the "median" spectrum ```{r} autoplot(sun.cosine.mspct, annotations = c("-", "peaks")) + geom_line(data = . %>% spct2fdata() %>% func.med.FM() %>% fdata2spct(), colour = "red", linewidth = 1, alpha = 0.67) ``` ### Highlight the "mean" spectrum ```{r} autoplot(sun.cosine.mspct, annotations = c("-", "peaks")) + geom_line(data = . %>% spct2fdata() %>% func.mean() %>% fdata2spct(), colour = "blue", linewidth = 1, alpha = 0.67) ``` ### Highlight the "trimmed mean" spectrum ```{r} autoplot(sun.cosine.mspct, annotations = c("-", "peaks")) + geom_line(data = . %>% spct2fdata() %>% func.trim.FM(trim = 0.25) %>% fdata2spct(), colour = "green", linewidth = 1, alpha = 0.67) ``` ## Functional one-way ANOVA We compare measurements done in parallel with two spectrometers with different entrance optics. Ten spectra were acquired once every 40 s inside a mixed forest next to a forest edge. A flat diffuser with cosine response and a hemispherical (or dome-shaped) diffuser with increased response at low angles of incidence are compared. ```{r} load("cos.A2.2m.series.spct.Rda") load("dome.A2.2m.series.spct.Rda") ``` As the two spectrometers even with identical configuration have slightly different wavelength calibrations for the individual pixels in the array, we need to re-express the spectra on identical wavelengths. We also trim the shortest and longest wavelengths as these spectral irradiance values are noisier. ```{r} when_measured(cos.A2.2m.series.spct)[c(1, 10)] cat(how_measured(cos.A2.2m.series.spct)) shade.cosine.mspct <- subset2mspct(cos.A2.2m.series.spct) %>% clean() %>% trim_wl(c(280, 1000)) %>% interpolate_mspct(length.out = 721) ``` ```{r} when_measured(dome.A2.2m.series.spct)[c(1, 10)] cat(how_measured(dome.A2.2m.series.spct)) shade.dome.mspct <- subset2mspct(dome.A2.2m.series.spct) %>% clean() %>% trim_wl(c(280, 1000)) %>% interpolate_mspct(length.out = 721) ``` We concatenate the collections of spectra. ```{r} shade.mspct <- c(cosine = shade.cosine.mspct, dome = shade.dome.mspct) names(shade.mspct) ``` A plot of all 20 spectra. ```{r} autoplot(shade.mspct) + aes(linetype = ifelse(grepl("cosine", spct.idx), "cosine", "dome")) + labs(linetype = "Entrance\noptics") ``` We test in a one-way ANOVA if the spectra differ, and they do differ significantly. ```{r} shade.fda <- mspct2fdata(shade.mspct) shade.anova <- fanova.onefactor(shade.fda, group = factor(rep(c("cosine", "dome"), each = 10L)), plot = TRUE) shade.anova$pvalue ``` To test if the shape of the spectra differs, we first scale them to equal area under the curve, and repeat the ANOVA analysis. ```{r} shade_scaled.fda <- mspct2fdata(fscale(shade.mspct, f = irrad)) shade_scaled.anova <- fanova.onefactor(shade_scaled.fda, group = factor(rep(c("cosine", "dome"), each = 10L)), plot = TRUE) shade_scaled.anova$pvalue ``` ## References ```{r} citation('fda.usc') ``` ```{r} citation('photobiology') ```