This vignette is designed to demonstrate how to use the curve fitting and sensitivity analysis tools Sections are named based on the set of methods to be used:

  1. Fitting light response curves

  2. Fitting CO2 response curves

  3. Fitting temperature response curves (Need data & to complete tutorial here)

  4. Fitting stomatal conductance models

  5. Fitting light respiration

  6. Fitting mesophyll conductance

  7. Fitting pressure-volume curves

  8. Fitting hydraulic vulnerability curves

  9. Sensitivity analyses (just need to think about measures of sensitivity & multi fits)

  10. Dependency checking

Components under construction:

  1. Full Gu-type CO2 response fitting

  2. alphag fitting

  3. Busch et al (2018) CO2 response model

Within each section, data will either be generated or used from an installed dataset within the package. For help with a given function, please consult the help file via: ?functionname in the console. If you want to know the fine details of the code, please go to:

https://github.com/jstinzi/photosynthesis

And look in the R folder to find the raw function files. These contain heavily annotated code that explains the why and how of their operation.

#Installing the package

You will need the following packages:

devtools - lets you install packages from Github and Bitbucket

minpack.lm - useful for nonlinear curve fitting that is more robust than base R

tidyverse - set of tools for manipulating data within R

FOR WINDOWS USERS

You will need to install Rtools, available at:

https://cran.r-project.org/bin/windows/Rtools/

#To install, run the following without comments
#library(devtools)
#install_github("jstinzi/photosynthesis")

#Load package
library(photosynthesis)
## Loading required package: ggplot2
## Loading required package: minpack.lm
## Loading required package: units
## udunits database from /usr/local/share/udunits/udunits2.xml
#To cite, use:
citation("photosynthesis")
## It is recommended to use 'given' instead of 'middle'.
## It is recommended to use 'given' instead of 'middle'.
## 
## To cite photosynthesis in publications use:
## 
##   Stinziano JR, Roback C, Gamble D, Murphy B, Hudson P, Muir CD.
##   (2020). photosynthesis: tools for plant ecophysiology & modeling. R
##   package version 2.0.1.
##   https://CRAN.R-project.org/package=photosynthesis.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Misc{,
##     title = {photosynthesis: tools for plant ecophysiology & modeling},
##     author = {Joseph R Stinziano and Cassaundra Roback and Demi Gamble and Bridget Murphy and Patrick Hudson and Christopher D Muir},
##     note = {R package version 2.0.1},
##     year = {2020},
##     url = {https://CRAN.R-project.org/package=photosynthesis},
##   }
#Load tidyr - needed for vignette manipulations
library(tidyr)

#Reading Li-Cor data

If you are trying to read in the raw data files of the Li-Cor 6400 or 6800 models, you can use the package RLicor by Erik Erhardt available on Github.

#library(devtools)
#install_github("erikerhardt/RLicor")
#library(RLicor)

#The following will detect and read Li-Cor 6400 and 6800 files
#?read_Licor

#To cite, use:
#citation("RLicor")

#1. Fitting light response curves

This package currently only implements the Marshall et al. 1980 non-rectangular hyperbola model of the photosynthetic light response.

#Read in your data
#Note that this data is coming from data supplied by the package
#hence the complicated argument in read.csv()
#This dataset is a CO2 by light response curve for a single sunflower
#Note that to read in your own data, you will need to delete the
#system.file() function, otherwise you will get an error
data <- read.csv(system.file("extdata", "A_Ci_Q_data_1.csv",
                             package = "photosynthesis"))

#Fit many AQ curves
#Set your grouping variable
#Here we are grouping by CO2_s and individual
data$C_s <-(round(data$CO2_s, digits = 0))

#For this example we need to round sequentially due to CO2_s setpoints
data$C_s <- as.factor(round(data$C_s, digits = -1))

#To fit one AQ curve
fit <- fit_aq_response(data[data$C_s == 600,],
                       varnames = list(A_net = "A",
                                         PPFD = "Qin",
                                       Q_cut = 250))

#Print model summary
summary(fit[[1]])
## 
## Formula: A_net ~ aq_response(k_sat, phi_J, Q_abs = data$Q_abs, theta_J) - 
##     Rd
## 
## Parameters:
##                 Estimate Std. Error t value Pr(>|t|)    
## k_sat          21.167200   0.158332  133.69 1.88e-08 ***
## phi_J.Q_abs     0.051907   0.001055   49.18 1.02e-06 ***
## theta_J         0.775484   0.014920   51.98 8.20e-07 ***
## Rd.(Intercept)  0.668495   0.065235   10.25 0.000511 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05535 on 4 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 1.49e-08
#Print fitted parameters
fit[[2]]
##         A_sat      phi_J   theta_J        Rd      LCP  resid_SSs
## k_sat 21.1672 0.05190746 0.7754836 0.6684953 12.97289 0.01225491
#Print graph
fit[[3]]

#Fit many curves
fits <- fit_many(data = data,
                 varnames = list(A_net = "A",
                                         PPFD = "Qin",
                                         group = "C_s"),
                 funct = fit_aq_response,
                 group = "C_s")
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |======================================================================| 100%
#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]])
## 
## Formula: A_net ~ aq_response(k_sat, phi_J, Q_abs = data$Q_abs, theta_J) - 
##     Rd
## 
## Parameters:
##                Estimate Std. Error t value Pr(>|t|)    
## k_sat          7.347423   0.141931  51.768 8.33e-07 ***
## phi_J.Q_abs    0.027192   0.001511  17.994 5.61e-05 ***
## theta_J        0.837778   0.030608  27.371 1.06e-05 ***
## Rd.(Intercept) 0.615283   0.086994   7.073  0.00211 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06799 on 4 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 1.49e-08
#Print the parameters
fits[[2]][[2]]
##          A_sat      phi_J   theta_J        Rd      LCP  resid_SSs
## k_sat 2.637157 0.01458002 0.8858892 0.5951635 42.17813 0.02446394
#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 dataframe for analysis
fits_pars <- compile_data(fits,
                          output_type = "dataframe",
                          list_element = 2)

#2. Fitting CO2 response curves

This package currently implements a Gu-type fitting procedure for CO2 response curves similar to the Duursma (2015) implementation. There is ongoing work to implement a full Gu-type method whereby mesophyll conductance, Km, and GammaStar could all be fit (Gu et al 2010). There is also ongoing work to implement a procedure to fit alphag for the TPU-limited region and to incorporate the Sharkey (2019) suggestion of using chlorophyll fluorescence data to inform TPU limitations.

#Read in your data
#Note that this data is coming from data supplied by the package
#hence the complicated argument in read.csv()
#This dataset is a CO2 by light response curve for a single sunflower
data <- read.csv(system.file("extdata", "A_Ci_Q_data_1.csv", 
                             package = "photosynthesis"))

#Define a grouping factor based on light intensity to split the ACi
#curves
data$Q_2 <- as.factor((round(data$Qin, digits = 0)))

#Convert data temperature to K
data$T_leaf <- data$Tleaf + 273.15

#Fit ACi curve. Note that we are subsetting the dataframe
#here to fit for a single value of Q_2
fit <- fit_aci_response(data[data$Q_2 == 1500, ],
                        varnames = list(A_net = "A",
                                      T_leaf = "T_leaf",
                                      C_i = "Ci",
                                      PPFD = "Qin"))

#View fitted parameters
fit[[1]]
##   Num V_cmax V_cmax_se    J_max        J      J_se V_TPU V_TPU_se        R_d
## 6   0 62.797  2.176227 110.3051 103.9718 0.1847135  1000       NA -0.3470509
##      R_d_se     cost citransition1 citransition2 V_cmax_pts J_max_pts V_TPU_pts
## 6 0.3947545 1.063979      427.6839      1450.485          8         4         0
##   alpha alpha_g gamma_star25 Ea_gamma_star K_M25   Ea_K_M  g_mc25 Ea_g_mc Oconc
## 6  0.24       0        42.75         37830 718.4 65508.28 0.08701       0    21
##   theta_J
## 6    0.85
#View graph
fit[[2]]
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 12 row(s) containing missing values (geom_path).