photosynthesis-curve-fitting-sensitivity-analyses
Joseph R. Stinziano
2020-09-11
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:
Fitting light response curves
Fitting CO2 response curves
Fitting temperature response curves (Need data & to complete tutorial here)
Fitting stomatal conductance models
Fitting light respiration
Fitting mesophyll conductance
Fitting pressure-volume curves
Fitting hydraulic vulnerability curves
Sensitivity analyses (just need to think about measures of sensitivity & multi fits)
Dependency checking
Components under construction:
Full Gu-type CO2 response fitting
alphag fitting
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 system database from /usr/local/share/udunits#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).
#View data with modeled parameters attached
#fit[[3]]
#Fit many curves
fits <- fit_many(data = data,
varnames = list(A_net = "A",
T_leaf = "T_leaf",
C_i = "Ci",
PPFD = "Qin"),
funct = fit_aci_response,
group = "Q_2")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%#Print the parameters
#First set of double parentheses selects an individual group value
#Second set selects an element of the sublist
fits[[3]][[1]]## Num V_cmax V_cmax_se J_max J J_se V_TPU V_TPU_se R_d
## 6 0 8.94862 0.5509706 47.01527 16.63315 0.08692268 1000 NA -0.1565895
## R_d_se cost citransition1 citransition2 V_cmax_pts J_max_pts
## 6 0.1264438 0.1194886 441.2967 1442.493 8 4
## V_TPU_pts alpha alpha_g gamma_star25 Ea_gamma_star K_M25 Ea_K_M g_mc25
## 6 0 0.24 0 42.75 37830 718.4 65508.28 0.08701
## Ea_g_mc Oconc theta_J
## 6 0 21 0.85#Print the graph
fits[[3]][[2]]## Warning: Removed 12 row(s) containing missing values (geom_path).
#Compile graphs into a list for plotting
fits_graphs <- compile_data(fits,
list_element = 2)
#Print graphs to pdf.
#print_graphs(data = fits_graphs,
# path = tempdir(),
# output_type = "pdf",
# pdf_filename = "mygraphs.pdf")
#Compile parameters into dataframe for analysis
fits_pars <- compile_data(fits,
output_type = "dataframe",
list_element = 1)#3. Fitting temperature response curves
This package provides support for multiple temperature response functions (Arrhenius 1915; Medlyn et al. 2002; Kruse & Adams. 2006; Heskel et al. 2016; Liang et al. 2018).
#Read in data
data <- read.csv(system.file("extdata", "A_Ci_T_data.csv",
package = "photosynthesis"),
stringsAsFactors = FALSE)
#Round temperatures to group them appropriately
#Use sequential rounding
data$T2 <- round(data$Tleaf, 1)
data$T2 <- round(data$Tleaf, 0)
#Look at unique values to detect rounding issues
unique(data$T2)