---
title: "Analysis of Real-World Gas-Exchange Data using Gasanalyzer"
output:
  rmarkdown::html_vignette:
    df_print: paged
    fig_caption: yes
  html_document:
    df_print: paged
  pdf_document:
    fig_caption: yes
vignette: >
  %\VignetteIndexEntry{gasanalyzer} 
  %\VignetteEncoding{UTF-8} 
  %\VignetteEngine{knitr::rmarkdown}
---

```{r setup, include = FALSE, eval = knitr::is_latex_output()}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "cairo_pdf")
# if you want to build vignettes, set NOT_CRAN to true in your .Renviron!
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")
```
```{r setup2, include = FALSE, eval = !knitr::is_latex_output()}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
# if you want to build vignettes, set NOT_CRAN to true in your .Renviron!
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")
```

```{r libs, message = FALSE}
library(gasanalyzer)
library(units)
```
```{r morelibs, message = FALSE, eval = NOT_CRAN, purl = NOT_CRAN}
library(ggplot2)
library(gridExtra)
```
```{r theming, include=FALSE, eval = NOT_CRAN, purl = NOT_CRAN}
old_options <- options(digits = 2)

theme_pub <- theme_bw() + theme(
  plot.title = element_text(size = 14),
  plot.subtitle = element_text(size = 12),
  axis.text = element_text(size = 10),
  axis.ticks.x.top = element_blank(),
  axis.title = element_text(size = 12),
  legend.text = element_text(size = 11),
  legend.title = element_text(size = 12),
  panel.grid.minor = element_blank())
theme_set(theme_pub)
```

## GasanalyzeR
This R package provides methods to read and process data from different portable photosynthesis systems. Variable names used in vendor-specific data files are translated into uniform names. Metadata, such as the instrument name and filename, are also added. For files that do not store timezone information, the local timezone is used by default.

Files from different instruments can be read. For the LI-COR 6800 gas-analyzer, both text files and the XML-based `.xlsx` files can be imported. In case `.xlsx` files are used, the equations stored in the files are converted to `R` and can be used to recalculate the data.

The three examples discussed below demonstrate the use of the package.

## 1. LI-6800 data and recalculations
The data used in this example is a CO~2~ response curve measured at low oxygen. Before loading the data files, the calibration factors for the instrument used to acquire the data must be imported. This is only necessary for some advanced functionality used later, and currently only available for calibration factors stored by the LI-6800 instrument.

The code below shows how to load the data from the `.txt` and `.xlsx` files.

```{r o2data, message=FALSE}
exampleFiles <- system.file("extdata", package = "gasanalyzer")
# import calibration factors for the 6800 used here:
import_factory_cals(exampleFiles)
# load data
xlsxfile <- read_6800_xlsx(paste(exampleFiles, "lowo2.xlsx", sep = "/"))
txtfile <- read_6800_txt(paste(exampleFiles, "lowo2", sep = "/"))
```

The LI-6800 confusingly stores the results of divisions by zero as zero in the `txt` file. After excluding columns containing the results of divisions by zero and the column with equations (a list column), there are very few differences between both files:

```{r}
# exclude list columns with equations, and divisions by zero
excludeCols <- vapply(xlsxfile, 
                      FUN=\(x) is.list(x) || 
                        all(is.infinite(x)) ||
                        all(is.nan(x)), TRUE)

# Interestingly, Leak.Fan reported in the txt file differs slightly
# from that in the xlsx. It is likely caused by rounding/averaging in the
# instrument firmware.
all.equal(txtfile[!excludeCols], xlsxfile[!excludeCols], tolerance = 1e-6)
```
Using plotting facilities available with `R`, relationships in the data can be  visualized. In the following example, `var2label()` was used to replace the variable names in the data with commonly-used symbols for the plot. The relation between net photosynthesis ($A$) and intercellular CO~2~ ($C_i$) is shown.

```{r plot1, fig.height = 4, fig.width = 5.5, eval = NOT_CRAN, purl = NOT_CRAN}
# default points:
update_geom_defaults("point", list(shape = 1, size = 3))

# convenience function to change plot labels automatically:
nice_labels <- function(plt) {
  plt$labels <- var2label(plt$labels)
  plt
}

#generate some descriptive stats  and plotmath it with a call to var2label:
descr <- colMeans(xlsxfile[c("LeafQ.Qin", "Meas.Tleaf", "Const.Oxygen")])
dlist <- var2label(rep(names(descr), 2),
                   use_units = rep(c(FALSE, TRUE), each = 3),
                   val = c(rep(NA, 3), as.character(descr)),
                   group = c("", "")) |>
    setNames(letters[1:6]) |> unlist() |>
    substitute(a==d*","~b==e*","~c==f, env = _ ) |> as.expression()

AvsCi1 <- xlsxfile |> ggplot(aes(GasEx.Ci, GasEx.A)) + geom_point() +
  labs(title = dlist)

nice_labels(AvsCi1)
```

If a researcher would like to make changes to some of constants used by the instrument for calculating derived variables, all data should be recalculated to take the effect into account. For example, an application note by the instrument vendor^[LI-COR Biosciences Inc (2022) Photosynthesis | Complementary Leaf Physiology Measurements: The LI-600 and LI-6800. [Link](https://www.licor.com/env/support/Photosynthesis/topics/complementary-photosynthesis-measurements.html)] shows how to use complimentary measurements with a porometer to determine the ratio of stomatal conductances between both sides of the leaves. For the leaf in this experiment, the ratio was 0.22. Such corrections typically have only very small effects because of the high boundary layer conductance in the leaf chambers that are commonly used for such experiments.

```{r}
gswOld <- xlsxfile$GasEx.gsw
xlsxfile <- xlsxfile |> 
  transform(Const.K = 0.22) |>
  recalculate()
#max effect:
set_units(max(1 - as.numeric(xlsxfile$GasEx.gsw / gswOld)) * 100, "%")
```
Gas-exchange measurements can be combined with measurements of chlorophyll fluorescence to obtain estimates of the electron transport rate  ($J_F$) through photosystem 2 (PSII). However, this requires knowledge of the amount of light energy absorbed by the leaf ($\alpha$), and the fraction of energy partitioned to PSII ($\beta$^[Genty B, Briantais J-M, Baker NR (1989) The relationship between the quantum yield of photosynthetic electron transport and quenching of chlorophyll fluorescence. Biochimica et Biophysica Acta (BBA) - General Subjects 990:87–92. [Link](https://doi.org/10.1016/S0304-4165(89)80016-9)]):

$$J_F = \alpha \beta \Phi_{PSII} Q_{in}$$

where $\Phi_{PSII}$ is the operating quantum efficiency of PSII, and $Q_{in}$ is the light intensity incident on the leaf.

Calculations implemented by the vendor typically assume $\beta$ is 0.5 and $\alpha$ is around 0.85 (possibly adjusted for the spectrum of the light source used). An interesting application of this electron transport estimate is to to back-calculate the CO~2~ concentration at the site of carboxylation, and by extension estimate the mesophyll conductance ($g_m$^[Harley PC, Loreto F, Di Marco G, Sharkey TD (1992) Theoretical Considerations when Estimating the Mesophyll Conductance to CO~2~ Flux by Analysis of the Response of Photosynthesis to CO~2~. Plant Physiol 98:1429–1436. [Link](https://doi.org/10.1104/pp.98.4.1429)]).

However, leaves vary in the amount of light absorbed. Moreover, not all of the absorbed light is used for photosynthetic carbon reduction. To address such concerns, an additional estimate for the electron transport rate ($J_C$) may be obtained by measuring gas exchange under non-photorespiratory conditions (i.e. at low oxygen). In this case, it follows from the C~3~ biochemical model^[Farquhar GD, von Caemmerer S, Berry JA (1980) A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149:78–90. [Link](https://doi.org/10.1007/BF00386231)] that:

$$J_C \approx 4(A + R_L) \to A \approx J_C/4 - R_L$$

where $R_{L}$ represents the dark respiration rate in the light. Assuming that alternative electron flows are minimal, the following relationship can be derived (Yin et al. 2009^[Yin X, van der Putten PEL, Belay D, Struik PC (2020) Using photorespiratory oxygen response to analyse leaf mesophyll resistance. Photosynth Res 144:85–99. [Link](https://doi.org/10.1007/s11120-020-00716-z)]):

$$A \approx J/4 - R_L = \alpha \beta \Phi_{PSII} Q_{in} / 4 - R_L$$

However, applying this assumptions on the data used here shows that either this model is too simple, or the assumed values for $\alpha$ and/or $\beta$ are incorrect:
```{r}
#confirm that the slope of A vs ETR is not 1:
xlsxfile |> 
    lm(GasEx.A ~ I(FLR.ETR/4), data = _) |> coef() |> _[[2]]
```

If the goal is to estimate $g_m$, it has been recommended (e.g. Yin et al. 2009) to not rely on the $\alpha \beta$ outputted by the instrument (`r xlsxfile$LeafQ.alpha[1]/2` in this case), but instead use a linear regression to find an "effective" value:

```{r}
# a linear regression:
lm1 <-
  xlsxfile |> transform(FLR.PhiQin_4 = FLR.phiPS2 * LeafQ.Qin / 4) |>
  lm(GasEx.A ~ FLR.PhiQin_4, data = _)

alphabeta <- lm1$coef[[2]]
```

The following plots the result of the regression:

```{r plot2, fig.height = 4, fig.width = 5.5, eval = NOT_CRAN, purl = NOT_CRAN}

AvsPhiQin1 <- ggplot(lm1$model, aes(x = .data[[names(lm1$model)[2]]],
                    y = .data[[names(lm1$model)[1]]])) +
    geom_point() +
    stat_smooth(method = "lm", formula = y ~ x) +
    labs(subtitle = bquote(alpha ~ beta == .(lm1$coef[[2]]) ~~~
                               R[L] == .(-lm1$coef[[1]])))

nice_labels(AvsPhiQin1)

```

The fitted value for the slope can now be used to adjust the constant value used by the instrument. The default equations used by the LI-6800 recalculate $\alpha$ based on the spectral composition of the light source. These equations are stored as `R` calls in the `gasanalyzer.Equations` column. They can also be manually extracted from a given `xlsx` file using `extract_6800_equations()` and edited using `modify_equations()`:
```{r}
# a long equation estimating the light absorption:
xlsxEqs <- xlsxfile$gasanalyzer.Equations[[1]]
print(xlsxEqs$LeafQ.alpha)

# add our custom alpha by dividing it through the used beta:
xlsxfile$Custom.alpha <- alphabeta / xlsxfile$Const.fPS2
# and change the equation to use the new value:
xlsxEqs <- modify_equations(xlsxEqs, LeafQ.alpha = \() {Custom.alpha})

# make a copy of the data and apply the new equations:
xlsxfileCopy <- xlsxfile |> 
  transform(gasanalyzer.Equations = list(xlsxEqs)) |>
  recalculate()

# confirm that the slope of A vs ETR is now 1:
xlsxfileCopy |> 
  lm(GasEx.A ~ I(FLR.ETR / 4), data = _) |> coef() |> _[[2]]

```
The `xlsx` equations imported from the file are typically sufficient for simple recalculations, but limited to the specific set of equations defined in the file. Moreover, in case only text-based data files are available, an alternative way to recalculate data is required. To handle calculations more comprehensively for all types of data, `gasanalyzer` includes a several sets of equations related to gas exchange, chlorophyll fluorescence and carbon isotope discrimination. To ensure greater reliability, units of the used quantities in these equations are automatically enforced by the `units` package^[Pebesma E, Mailund T, Hiebert J (2016) Measurement Units in R. R J 8:486–494. [Link](https://doi.org/10.32614/RJ-2016-061)]. In the code below, the data is recalculated using the default set, with additional instrument specific computations for the LI-6800. In addition, a several equations are added to take the oxygen concentration stored in the data into account when recalculating the gas mole fractions. Finally, $\alpha$ is set to the previously fitted value.

```{r tryeqs}
RecalEqs <- create_equations(c("li6800", "default", "O2_correction"), 
                             LeafQ.alpha = \() {Custom.alpha})
xlsxfileRecal <- recalculate(xlsxfile, RecalEqs)

all.equal(xlsxfileRecal[names(xlsxfile)], xlsxfile, tol = 0.001)
```
Small differences are found in variables as a result of the new $\alpha$ and in a value related to the energy balance. The latter is a result of a more precise conversion from Celsius to Kelvin units. Equations can be checked and modified before use, and custom equations can be added as `R` functions:

```{r eqs}
# providing a bad name will show a warning with valid names:
try <- create_equations("help")

# an equation set for gm, with a custom equation for resistance:
create_equations("gm_fluorescence", FLR.rm = \(x) {1 / FLR.gm})
```

The calibration of infra-red gas analyzers is sensitive to the gas-mixture used for the measurements. Thus, when doing measurements are done under low oxygen, this has to be taken into account. Analyzers produced by LI-COR Inc. or Heinz Walz GmbH allow an oxygen concentration to be specified or measured and adjust the reported concentrations based on several calibration factors. These equations are unfortunately not included in the stored data files, but are included when either the `O2_correction` or `raw` flags are passed to `create_equations()` as done above. The latter allows for recalculation from raw absorptance values but requires such data to be stored during the measurement, which is not done by default. The `O2_correction` flag back-calculates such raw values based on the original O~2~ concentration specified in the data, and then recalculates gas mole fractions based on a potentially updated O~2~ concentration. It must be added that since the data reported by the instrument may be averaged before storing, recalculating the data at different oxygen levels may not always be exact.

The equations depend on calibration values that depend on the instrument type, and it may be a good idea to verify the calibration factors. For the current example, the calibration factors for the used instrument need to be imported before reading the data (as was mentioned at the top of this section).

To investigate how much the results would differ when a user forgot to enter the correct O~2~ concentration during the measurement, I have recalculated the data using an oxygen concentration of 21 \%:
```{r}
xlsxfile21 <- xlsxfileRecal |> 
  transform(Const.Oxygen = set_units(21, "%")) |>
  recalculate(RecalEqs)
```

And plot the results:
```{r plot3, fig.height = 3, fig.width = 6, eval = NOT_CRAN, purl = NOT_CRAN}
AvsCi21 <- xlsxfile21 |> ggplot() + 
  geom_point(aes(GasEx.Ci, GasEx.A, shape = "21% O2")) +
  geom_point(aes(xlsxfileRecal$GasEx.Ci,
                 xlsxfileRecal$GasEx.A, shape = "1% O2")) +
  scale_shape_manual(name = "Calculated at:",
                     values = c("1% O2" = 1, "21% O2" = 2)) +
  theme(legend.position = "inside", legend.position.inside = c(0.7, 0.3))

AvsPhiQin21 <-
  xlsxfile21 |>
  lm(GasEx.A ~ FLR.PhiQin_4, data = _) |> (\(x) {
    ggplot(x$model, aes(x = .data[[names(x$model)[2]]],
                        y = .data[[names(x$model)[1]]])) +
      geom_point() +
      stat_smooth(method = "lm", formula = y ~ x) +
      labs(subtitle = bquote(alpha ~ beta == .(x$coef[[2]]) ~~~
                            R[L] == .(-x$coef[[1]]))) +
      geom_point(data = xlsxfileRecal, shape = 2) 
    })()

grid.arrange(nice_labels(AvsCi21), nice_labels(AvsPhiQin21), ncol = 2)
```

The figure confirms that since the correction mainly affects water concentrations, variables such as $g_{sw}$ and therefore also $C_i$ are most affected. However the effect on $C_i$ remained below 5\%. It is noteworthy that the effect of the oxygen correction is negligible for the assimilation rate, and therefore also does not affect the relation between $A$ and $J_F$:
```{r}
xlsxfileRecal |> 
    lm(GasEx.A ~ I(FLR.ETR / 4), data = _) |> coef()

xlsxfile21 |> 
    lm(GasEx.A ~ I(FLR.ETR / 4), data = _) |> coef() 
```

## 2. GFS-3000, loading multiple files and doing a sensitivity analysis
Using base R functions, multiple files can easily be concatenated and used with `gasanalyzer`:
```{r loadACi}
aciFiles <- list.files(exampleFiles, "aci\\d\\.csv", full.names = TRUE)

ACi <- lapply(aciFiles, \(x) {print(basename(x)); read_gfs(x, tz = "CEST")}) |>
   do.call("rbind", args = _)

data.frame(aggregate(list(RH = ACi$GasEx.RHcham), 
                     list(file = ACi$SysObs.Filename), mean),
           RH.SD = aggregate(list(ACi$GasEx.RHcham), 
                     list(ACi$SysObs.Filename), sd)[, 2]) |>
  knitr::kable()

```
The $C_i$ estimates provided by the GFS-3000 instrument do not take boundary layer conductance into account. `gasanalyzer` provides the possibility to recalculate the data using predefined sets of equations. The boundary layer conductance ($g_{bw}$) of the default chamber at different fan speeds was provided by the instrument vendor (Heinz Walz GmbH). However, the effect on $C_i$ is small when the stomatal conductance ($g_s$) is much lower than $g_{bw}$:

```{r}
gfsEqs <- create_equations(c("default", "gfs3000"))
ACi2 <- recalculate(ACi, gfsEqs)
all.equal(ACi$GasEx.Ci, ACi2$GasEx.Ci)
```
The relationship between $A_n$ and $C_i$ can be used to obtain *in vivo* estimates for parameters related to biochemical processes underlying photosynthesis. The R package `photosynthesis`^[Stinziano JR, Roback C, Sargent D, et al (2021) Principles of resilient coding for plant ecophysiologists. AoB PLANTS 13:plab059. [Link](https://doi.org/10.1093/aobpla/plab059)] implements methods to obtain such estimates. Recently it has been suggested^[Márquez DA, Stuart-Williams H, Farquhar GD (2021) An improved theory for calculating leaf gas exchange more precisely accounting for small fluxes. Nat Plants 7:317–326. [Link](https://doi.org/10.1038/s41477-021-00861-w)] that the $C_i$ values calculated using equations implemented in most gas-exchange instruments are biased because  the underlying theory did not account for gas fluxes through the cuticula. Since estimates for cuticular conductance are not widely available, it can be useful to perform a sensitivity analysis on gas-exchange data to investigate the possible consequences of limiting cuticular conductance. A first assessment of the importance of including cuticular conductance was recently published^[Lamour J, Davidson KJ, Ely KS, et al (2022) New calculations for photosynthesis measurement systems: what’s the impact for physiologists and modelers? New Phytol 233:592–598. [Link]( https://doi.org/10.1111/nph.17762)]. Here, an example is shown that includes a similar analysis, but also test for the effect of potential bias in assumed values for the ratio of stomatal conductances between both sides of the leaf (`Const.K`), cuticular conductance (`Const.gcw`), assumed relative humidity in the stomatal cavities (`Const.RHi`), and the measured values for the leaf temperature (`Meas.Tleaf`) and H~2~O mol fraction reported by the instrument (`Meas.H2Or`). The `permutate()` method was used to create a range of values for these quantities while keeping all other measured and assumed variables constant. Subsequently all derived variables are recomputed:
```{r}
# we define some constants that specify the relative value of the variables
# on which we will do a sensitivity analysis:
SAset <- c("Const.calk", "Const.calg", "Const.calrhi", "Const.calt",
           "Const.calh")
# symbols for the plot header:
SAlabs <- c("italic(K)", "italic(g)[cw]", "RH[i]", 
            "italic(T)[leaf]", "'['*H[2]*O*']'[r]")
# set defaults to 1 (100%) and define sequences over which to permutate
ACi[SAset] <- 1
# a 5% range:
seq5 <- seq(-0.05, 0.05, length.out = 11)
seq5 <- seq5[seq5 != 0]

# make equations allowing to vary the variables by a relative value:
MSFeqs <- create_equations(c("default", "gfs3000", "cuticular_conductance"), 
                           Const.K = function() { Const.calk * Const.K },
                           Meas.Tleaf = function() { Const.calt * Meas.Tleaf },
                           Const.gcw = function() { Const.calg * Const.gcw },
                           Const.gcc = function() { Const.gcw / 20 },
                           Const.RHi = function() { Const.calrhi * Const.RHi },
                           Meas.H2Or = function() { Const.calh * Meas.H2Or },
                           Meas.H2Os = function() { Const.calh * Meas.H2Os })

# set gcw and use steady-state equations, then use the sequences
ACi[["Const.gcw"]] <- set_units(25, "mmol*m^-2*s^-1")
ACi[["Const.RHi"]] <- set_units(95, "%")
ACi[["SysConst.UseDynamic"]] <- FALSE

ACi.SA <- rbind(permutate(ACi, Const.calk = c(1, 1 + 20 * seq5)),
              permutate(ACi, Const.calg = 1 + 20 * seq5),
              permutate(ACi, Const.calrhi = 1 + seq5),
              permutate(ACi, Const.calt = 1 + seq5),
              permutate(ACi, Const.calh = 1 + seq5)) |>
  recalculate(MSFeqs)
```

In the next step, the `photosynthesis` package was used to obtain estimates for the carboxylation capacity ($V_{cmax}$), electron transport rate ($J$), triose-phosphate limitation ($T_p$) and the respiration rate in the light ($R_L$). To facilitate  analysis of the results, the following function takes care of some intermediate steps with the curve-fitting and provides convenient tabular output:

```{r fitaci, eval = NOT_CRAN, purl = NOT_CRAN}
library(photosynthesis)

aPars <- c("V_cmax", "J_max", "V_TPU", "R_d")

fitaci_per_group <- function(df, aPars, g1, ...) {
  grp <- c(g1, ...)
  outCols <- c(grp, aPars, "group", "Plots")
  splitOn <- paste("~", paste(grp, collapse = " + "))
  # drop category name and Filename for shorter labels:
  lbls <- vapply(strsplit(grp, ".", fixed = TRUE), \(x) rev(x), c("nm","gp")) |>
    paste0(": ") |>
    gsub("Filename: ", "", x = _, fixed = TRUE) 
  # layout:
  plotAdd <- list(do.call(labs, var2label(list(y = "GasEx.A", x = "GasEx.Ci"),
                                         TRUE)),
                  theme(legend.position = "none", plot.title.position = "plot",
                        plot.title = element_text(size = 8)))
  # translate names to what the photosynthesis package uses (sadly not essdive):
  out <- df[c("GasEx.A", "GasEx.Ci", "LeafQ.Qin", "Meas.Tleaf", 
              "Meas.Pa", "gasanalyzer.Permutate", grp)] |>
    setNames(c("A_net", "C_i", "PPFD", "T_leaf", "Pa", "group", grp)) |>
    drop_units() |> transform(T_leaf = T_leaf + 273.15) |>
    split(as.formula(splitOn), drop = TRUE) |>
    lapply(\(x) {
      lbl <- paste0(lbls, x[1, grp], collapse=", ")
      fit <- try(fit_aci_response(x, Pa = mean(x$Patm)))
      if (!inherits(fit, "try-error"))
        data.frame(x[1, c("group", grp)], 
                   fit[["Fitted Parameters"]],
                   Plots = I(list(fit$Plot + ggtitle(lbl) + plotAdd))
                   )[outCols]
      else 
        data.frame()
    }) 
  do.call(rbind, c(out, make.row.names = FALSE)) 
}
```

Now, we will call this function and transform the output to relative values and average them over the 3 plants:

```{r eval = NOT_CRAN, purl = NOT_CRAN}
# call the function to fit the curves:
CO2curves.SA <- ACi.SA |>
  fitaci_per_group(aPars, "SysObs.Filename", SAset)

# reference values are where the SAset variables equal 1:
refids <- rowSums(CO2curves.SA[SAset] == 1) == length(SAset)
# combine all relative values in one column based on group:
CO2curves.SA$gval <- unlist(CO2curves.SA[cbind(seq_len(nrow(CO2curves.SA)),
                                        match(CO2curves.SA$group,
                                              names(CO2curves.SA)))])
# value as relative change:
CO2curves.SA$gval <- 100 * (as.numeric(CO2curves.SA$gval) - 1)

# calc mean over the 3 reps:
CO2curves.SAm <- 
  lapply(names(CO2curves.SA[aPars]),
         \(x) {
           # relative change:
           rval <- 100 * as.numeric((CO2curves.SA[[x]] - 
                                       CO2curves.SA[[x]][refids]) /
                                      abs(CO2curves.SA[[x]][refids]))
           # grouping:
           byl <- as.list(CO2curves.SA[SAset])
           with(CO2curves.SA, 
                data.frame(nm = x, 
                           aggregate(cbind(gval), byl, head, n = 1),
                           group = aggregate(group, byl, head, n = 1)[["x"]],
                           val = aggregate(rval, byl, mean)[["x"]],
                           sd = aggregate(rval, byl, sd)[["x"]]))}) |>
  do.call(rbind, args = _)
```

The model fits computed by the `photosynthesis` package are stored in the data and can be plotted:
```{r fig.height = 2.5, fig.width = 7, warning = FALSE, eval = NOT_CRAN, purl = NOT_CRAN}
grid.arrange(grobs = CO2curves.SA$Plots[refids], ncol = 3)
```

The code below creates a combined plot showing how a relative change in the 5 different variables affects the parameters of the C~3~ biochemical model:

```{r eval=NOT_CRAN, fig.height = 4.5, fig.width = 7, purl=NOT_CRAN}
# make lbls factors for good order and plotmath
CO2curves.SAm$group <- factor(CO2curves.SAm$group,
                              levels = SAset,
                              labels = SAlabs)
CO2curves.SAm$nm <- factor(CO2curves.SAm$nm,
                           levels = c("V_cmax", "J_max", "V_TPU", "R_d"),
                           labels = c("italic(V)[cmax]", "italic(J)[max]",
                                      "italic(V)[TPU]", "italic(R)[L]"))
# and plot:
CO2curves.SAm |>
  ggplot(aes(x = gval, y = val, group = nm)) + 
  geom_vline(aes(xintercept = 0), linetype = "dotdash", color = "darkgrey") +
  geom_point() +
  geom_ribbon(aes(ymin = val - sd, ymax = val + sd), alpha = 0.1) +
  facet_grid(cols = vars(group), rows = vars(nm),
             scales="free", labeller = label_parsed) +
  xlab("Relative change in variable (%)") + 
  ylab("Relative change in parameter (%)") +
  theme(plot.title = element_text(size = 14),
        panel.grid.minor = element_blank(),
        panel.spacing.x = unit(14, "pt"))


```

The results show that for these measurements, accounting for stomatal or cuticular conductance has only very small effects on the estimates. Bias in leaf temperature measurements have much larger effects.

## 3. Combining gas-exchange data with isotope discrimination

The ability to quick analyze the sensitivity of a result to different constants is also valuable for calculations of mesophyll conductance ($g_m$) from concurrent measurements of gas exchange and carbon isotope discrimination. As an example, a data file from an experiment combining gas exchange with carbon isotope measurements was included. This file contains additional columns with the $\delta^{13}\textrm{C}$ value of the reference and sample air, and of the air in the glasshouse where the plants were grown. 

The dataset can be used to estimate the mesophyll conductance. Such calculations require several fractionation factors describing how the abundance ratio of the carbon isotopes will change using a physiological or physical process. Since uncertainty exist surrounding the exact value of such factors, it is useful to analyze the data under several different assumptions. In addition, the effect of different assumptions in the isotopic model (see Busch et al. 2020^[Busch FA, Holloway-Phillips M, Stuart-Williams H, Farquhar GD (2020) Revisiting carbon isotope discrimination in C~3~ plants shows respiration rules when photosynthesis is low. Nat Plants 6:245–258. [Link](https://doi.org/10.1038/s41477-020-0606-6)]) can be tested.

In the following example, `gasanalyzer` is used to analyze the effect of the fractionation factor associated with photorespiration ($f$) using both the connected and disconnected isotopic models as described by Busch et al. 2020.

```{r d13Crecalc}
# load data
isoFile <- file.path(exampleFiles, "d13C.tsv")
# read and recalculate
isotopes <- read_gasexchange(isoFile) |>
  recalculate(create_equations(c("default", "li6400")))

# vary f between 0 and 20, and use 2 sets of equations:
isotopes$d13CConst.e <- set_units(0, "permille")
isotopes$gasanalyzer.Equations <- list(NA)
isotopes <- isotopes |>
  permutate(d13CConst.f = seq(0, 20, 0.5)) |> 
  permutate(gasanalyzer.Equations = 
              c(connected = list(create_equations("d13C")),
                disconnected = list(create_equations(c("d13C", 
                                                       "d13C_dis"))))) |>
  recalculate() |>
  transform(factor.f = factor(d13CConst.f))

# point method gm is already calculated, average the results:
gmEst <- aggregate(list(gm = as.numeric(isotopes$d13C.gm)),
                   list(d13CConst.f = isotopes$d13CConst.f,
                        Model = isotopes$gasanalyzer.PermutateLabel),
                   FUN = function(x) 
                     c(d13C.gm = mean(x), d13C.gm.sem = sd(x)
                       / sqrt(length(x))))
gmEst <- do.call(data.frame, c(gmEst[1:2], unname(gmEst[3])))
gmEst$method <- "point"

```

Before plotting the results, it is useful to also derive $g_m$ estimates using the so-called slope method (see e.g. Tazoe et al. 2009^[Tazoe Y, von Caemmerer S, Badger MR, Evans JR (2009) Light and CO2 do not affect the mesophyll conductance to CO2 diffusion in wheat leaves. J Exp Bot 60:2291–2301. [Link](https://doi.org/10.1093/jxb/erp035)]. For this approach, the relationship between the difference between modeled and observed discrimination is plotted against the ratio of net assimilation and the partial pressure of CO~2~ in the ambient air ($A/p_{C_a}$):
```{r d13Cslopes, fig.height = 4, fig.width = 5.5, eval = NOT_CRAN, purl = NOT_CRAN}

d13Cslopes <- isotopes |>
  subset(gasanalyzer.PermutateLabel == "connected" &
           as.numeric(d13CConst.f) %in% c(7, 11, 16)) |>
  ggplot(aes(y = d13C.DeltaiDeltao, x = d13C.A_pCa, color = factor.f)) + 
  geom_point() +  
  geom_smooth(method = "lm", formula = y ~ x, aes(fill = factor.f)) +
    guides(fill = "none", 
         colour = guide_legend(override.aes = list(fill = NA))) +
  labs(color = var2label("d13CConst.f", TRUE)[[1]]) +
  theme(legend.position = "inside", legend.position.inside = c(0.75, 0.10),
        legend.direction = "horizontal")

nice_labels(d13Cslopes)

```

The mesophyll conductance is inversely related to the slope of this relationship. The slopes appear quite similar, suggesting that the slope method is much less sensitive to assumptions on the value taken for $f$. However, a more thorough analysis requires a nonlinear fitting approach.

Such an approach has an additional advantage, because it allows for the estimation of the effective fractionation associated with mitochondrial respiration ($e'$).
Respiration typically releases carbon from substrate that was assimilated before the leaf was clamped in the gas-exchange chamber, and since the isotopic composition of the tank air typically differs drastically from that of the air during growth conditions, this may result in a strong isotopic signal that needs to be accounted for. The correction proposed by Busch et al. 2020 could not be used for this analysis because the discrimination against ${}^{13}\textrm{C}$ under growth conditions was not quantified for these plants. However, for the slope method, it is possible to fit both $e'$ and $g_m$ in one analysis (Tazoe et al. 2009), making the results free of any assumptions regarding the actual effect of respiratory fractionations.

Below we will use the isotope equations for the connected and disconnected model provided by `gasanalyzer` for the estimation of $g_m$ using the slope method:

```{r eval=NOT_CRAN, purl=NOT_CRAN}
# The equation for the slope method is already in the data, so this function
# just evaluates that. Its not the fastest way, because recalculate adds
# overhead and unneccessary calcs, but simple:
fitfunc <- function(gm, ep) {
    # Di in the equation doesn't take gm into account, and assumes Ci=Cc
    # to get the true Delta we so we modify Ci to be Cc again:
    df$GasEx.Ci <-  df$GasEx.Ci - df$GasEx.A / 
      (set_units(gm, "mol*m^-2*s^-1*bar^-1") * df$Meas.Pa)
    # substitute estimated ep in data and recalc:
    df$d13C.ep <- set_units(ep, "permille")
    as.numeric(recalculate(df)$d13C.Deltai)
}

#this will take a while:
gmEstSlope <- isotopes |> 
  #we want to fit ep, so modify the equation list to not recalculate it:
  transform(gasanalyzer.Equations = 
              lapply(gasanalyzer.Equations,
                     \(x) modifyList(x, list(d13C.ep = NULL)))) |>
  #split by f and equation version, and call nls: 
  split(~ isotopes$d13CConst.f + isotopes$gasanalyzer.PermutateLabel) |>
  lapply(function(dat) {
           # data is moved to environment of fitfunc
           environment(fitfunc) <- list2env(list(df = dat))
           nlsfit <- nls(drop_units(d13C.Deltao) ~ fitfunc(gm, ep),
                         start = c(gm = 0.3, ep = -10), dat)
           cfs <- summary(nlsfit)$coefficients
           data.frame(d13CConst.f = dat$d13CConst.f[1],
                      Model = dat$gasanalyzer.PermutateLabel[1],
                      d13C.gm = cfs[1],
                      d13C.gm.sem = cfs[3], d13C.ep = cfs[2], 
                      d13C.ep.sem = cfs[4])}) |>
  do.call(rbind, args = _) |>
  transform(method = "slope")

gmEst <- rbind(gmEst, gmEstSlope[names(gmEstSlope) %in% names(gmEst)],
               make.row.names = FALSE)

#re-adds units:
units(gmEst$d13C.gm) <- deparse_unit(isotopes$d13C.gm)
units(gmEst$d13C.gm.sem) <- deparse_unit(isotopes$d13C.gm)

```

Plot the results of the analysis:

```{r fsens, fig.height = 4, fig.width = 7, eval = NOT_CRAN, purl = NOT_CRAN}

fsens <- ggplot(gmEst, aes(d13CConst.f, d13C.gm, 
                           ymin = d13C.gm - d13C.gm.sem,
                           ymax = d13C.gm + d13C.gm.sem, 
                           shape = method, color = Model)) + 
  geom_point(color = "black") +
  geom_ribbon(alpha = 0.1, color = NA) + guides(color = "none") +
  scale_shape_manual("Method:", values = c(1, 2)) + facet_wrap(vars(Model)) +
  theme(legend.position = "inside", legend.position.inside = c(0.12, 0.85),
        panel.grid.minor = element_blank()) 
nice_labels(fsens)
```

The results indicate that the disconnected model resulted in small differences between the slope and point methods. Although all approaches showed some increase in the $g_m$ estimate when higher values for photorespiratory fractionation $f$ were used, the connected model lead to an unusual strong increase in the estimates when higher values for $f$ were used.


```{r include=FALSE, eval = NOT_CRAN, purl = NOT_CRAN}
options(old_options)
```
