Fitting light response curves
In photosynthesis: Tools for Plant Ecophysiology & Modeling

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)

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 and plot the fitted curve. This allows users the ability to more easily customize their figures.

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 that are already pretty good.

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 package which fits models in Stan. 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.

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)

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.

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.