---
title: "Fitting light response curves"
author: "Joseph R. Stinziano and Christopher D. Muir"
date: "`r Sys.Date()`"
output: rmarkdown::html_document
vignette: >
 %\VignetteIndexEntry{Fitting light response curves}
 %\VignetteEngine{knitr::rmarkdown}
 %\VignetteEncoding{UTF-8}
---

## Preferred version (**photosynthesis** >= 2.1.1)

This package currently only implements the Marshall & Biscoe (1980) non-rectangular 
hyperbola model of the photosynthetic light response.

### Fit the light-response curve using nonlinear least-squares (nls)

```{r, message = FALSE}

library(broom)
library(dplyr)
library(photosynthesis)

# Read in your data
dat = system.file("extdata", "A_Ci_Q_data_1.csv", package = "photosynthesis") |>
  read.csv() |>
  # Set grouping variable
  mutate(group = round(CO2_s, digits = 0)) |>
  # For this example, round sequentially due to CO2_s set points
  mutate(group = as.factor(round(group, digits = -1)))

# Fit one light-response curve
fit = fit_photosynthesis(
  .data = filter(dat, group == 600),
  .photo_fun = "aq_response",
  .vars = list(.A = A, .Q = Qabs),
)

# The 'fit' object inherits class 'nls' and many methods can be used

## Model summary:
summary(fit)

## Estimated parameters:
coef(fit)

## 95% confidence intervals:
confint(fit)

## Tidy summary table using 'broom::tidy()'
tidy(fit, conf.int = TRUE, conf.level = 0.95)

## Calculate light compensation point
coef(fit) |>
  t() |>
  as.data.frame() |>
  mutate(LCP = ((Rd) * (Rd * theta_J - k_sat) / (phi_J * (Rd - k_sat)))) |>

## Calculate residual sum-of-squares
sum(resid(fit) ^ 2)

```

## Plot model fit and raw data

The deprecated function `fit_aq_response()` generated a figure automatically, but it used `geom_smooth()` rather than plotting the model fit. We now prefer to use generic methods from the package [**ggplot2**](https://ggplot2.tidyverse.org/) and plot the fitted curve. This allows users the ability to more easily customize their figures.

```{r}
#| fig.alt: >
#|   Example of fitted light-response curve

library(ggplot2)

b = coef(fit)

df_predict = data.frame(Qabs = seq(0, 0.84 * 1500, length.out = 100)) |>
  mutate(
    A = marshall_biscoe_1980(
      Q_abs = Qabs,
      k_sat = b["k_sat"],
      b["phi_J"],
      b["theta_J"]
    ) - b["Rd"]
  )

ggplot(mapping = aes(Qabs, A)) +
  geom_line(data = df_predict) +
  geom_point(data = filter(dat, group == 600)) +
  labs(
    x = expression("Irradiance (" * mu * mol ~ m^{-2} ~ s^{-1} * ")"),
    y = expression(A[net] ~ "(" * mu * mol ~ m^{-2} ~ s^{-1} * ")")
  ) +
  theme_bw()

```

## Fit multiple curves with **photosynthesis** and **purrr**

In the previous version, we used `fit_many()` to fit many light-response curves simultaneously. We now prefer to use generic methods from the package [**purrr**](https://purrr.tidyverse.org/) that are already pretty good.

```{r}
library(purrr)

fits = dat |>
  split(~ group) |>
  map(fit_photosynthesis, .photo_fun = "aq_response", .vars = list(.A = A, .Q = Qabs))

## Estimated parameters:
fits |>
  map(coef) |>
  map(t) |>
  map(as.data.frame) |>
  imap_dfr(~ mutate(.x, CO2_s = .y))

```

## Fit Bayesian light-response curves with **brms** and *Stan*

Traditional model fitting use a nonlinear least-squares approach, but Bayesian methods have some advantages, especially with more complex data sets. We added an option to fit a single Bayesian light-response curve using the amazing [**brms**](https://CRAN.R-project.org/package=brms) package which fits models in [*Stan*](https://mc-stan.org/). We have not implemented more complex approaches (e.g. multilevel light-response models) because once you are doing that, it's probably easier to code the model directly into **brms** functions. Hopefully this code can get you started though. We have not run the example below, but copy-and-paste into your *R* Console to try.

```{r, eval = FALSE}

fit = fit_photosynthesis(
  .data = filter(dat, group == 600),
  .photo_fun = "aq_response",
  .vars = list(.A = A, .Q = Qabs),
  .method = "brms",
  brm_options = list(chains = 1)
)

# The 'fit' object inherits class 'brmsfit' and many methods can be used
summary(fit)

```



### Fit the light-response curve with photoinhibition using nonlinear least-squares (nls)

```{r, message = FALSE}
#| fig.alt: >
#|   Example of fitted light-response curve with photoinhibition

library(broom)
library(dplyr)
library(tibble)
library(photosynthesis)

# Simulate data
set.seed(202411123)

## Parameters 
k_sat = 30
phi_J = 0.05
theta_J = 0.8
Rd = 3
b_inh = 0.003

## Vector of light values
Q = c(0, 25, 50, 75, 100, 125, 150, 200, 300, 400, 500, 600, 700, 800, 1000, 1200, 1400, 1600, 1800, 2000)

dat = tibble(
  Qabs = 0.84 * Q,
  A = photoinhibition(Qabs, k_sat, phi_J, theta_J, b_inh) - Rd + rnorm(length(Qabs), 0, 0.1)
)

# Fit one light-response curve
fit = fit_photosynthesis(
  .data = dat,
  .photo_fun = "aq_response",
  .model = "photoinhibition",
  .vars = list(.A = A, .Q = Qabs),
)

# The 'fit' object inherits class 'nls' and many methods can be used

## Model summary:
summary(fit)

## Estimated parameters:
coef(fit)

## 95% confidence intervals:
confint(fit)

## Tidy summary table using 'broom::tidy()'
tidy(fit, conf.int = TRUE, conf.level = 0.95)

```

## Plot model fit and raw data

```{r}

library(ggplot2)

b = coef(fit)

df_predict = data.frame(Qabs = seq(0, 0.84 * 2000, length.out = 100)) |>
  mutate(
    A = photoinhibition(
      Q_abs = Qabs,
      k_sat = b["k_sat"],
      b["phi_J"],
      b["theta_J"],
      b["b_inh"]
    ) - b["Rd"]
  )

ggplot(mapping = aes(Qabs, A)) +
  geom_line(data = df_predict) +
  geom_point(data = dat) +
  labs(
    x = expression("Irradiance (" * mu * mol ~ m^{-2} ~ s^{-1} * ")"),
    y = expression(A[net] ~ "(" * mu * mol ~ m^{-2} ~ s^{-1} * ")")
  ) +
  theme_bw()

```

## Deprecated version (**photosynthesis** <= 2.1.1)

The `fit_aq_response()` function is the original version, but we are no longer updating it and may phase it out of future releases. Use `fit_photosynthesisi(..., .photo_fun = 'aq_response')` instead.

```{r, message = FALSE}

library(dplyr)
library(photosynthesis)
# Read in your data
dat = system.file("extdata", "A_Ci_Q_data_1.csv", package = "photosynthesis") |>
  read.csv() |>
  # Set grouping variable
  mutate(group = round(CO2_s, digits = 0)) |>
  # For this example, round sequentially due to CO2_s setpoints
  mutate(group = as.factor(round(group, digits = -1))) |>
  rename(A_net = A, PPFD = Qin)

# To fit one AQ curve
fit = fit_aq_response(filter(dat, group == 600))

# Print model summary
summary(fit[[1]])

# Print fitted parameters
fit[[2]]

# Print graph
fit[[3]]

# Fit many curves
fits = fit_many(
  data = dat,
  funct = fit_aq_response,
  group = "group",
  progress = FALSE
)

# Look at model summary for a given fit
# First set of double parentheses selects an individual group value
# Second set selects an element of the sublist
summary(fits[[3]][[1]])

# Print the parameters
fits[[2]][[2]]

# Print the graph
fits[[3]][[3]]

#Compile graphs into a list for plotting
fits_graphs = compile_data(fits, list_element = 3)

# Print graphs to jpeg
# print_graphs(data = fits_graphs, path = tempdir(), output_type = "jpeg")

#Compile parameters into data.frame for analysis
fits_pars = compile_data(fits, output_type = "dataframe", list_element = 2)

```

# References

Marshall B, Biscoe P. 1980. A model for C3 leaves describing the dependence of net photosynthesis on irradiance. *Journal of Experimental Botany* 31:29-39.