Generate an error function for C3 Variable J curve fitting
error_function_c3_variable_j.Rd
Creates a function that returns an error value (the sum of squared residuals)
representing the amount of agreement between modeled and measured An
values.
Usage
error_function_c3_variable_j(
replicate_exdf,
fit_options = list(),
sd_A = 1,
atp_use = 4.0,
nadph_use = 8.0,
curvature_cj = 1.0,
curvature_cjp = 1.0,
a_column_name = 'A',
ci_column_name = 'Ci',
j_norm_column_name = 'J_norm',
kc_column_name = 'Kc',
ko_column_name = 'Ko',
oxygen_column_name = 'oxygen',
phips2_column_name = 'PhiPS2',
qin_column_name = 'Qin',
rl_norm_column_name = 'RL_norm',
total_pressure_column_name = 'total_pressure',
vcmax_norm_column_name = 'Vcmax_norm',
cj_crossover_min = NA,
cj_crossover_max = NA,
hard_constraints = 0,
require_positive_gmc = 'positive_a',
gmc_max = Inf
)
Arguments
- replicate_exdf
An
exdf
object representing one CO2 response curve.- fit_options
A list of named elements representing fit options to use for each parameter. Values supplied here override the default values (see details below). Each element must be
'fit'
,'column'
, or a numeric value. A value of'fit'
means that the parameter will be fit; a value of'column'
means that the value of the parameter will be taken from a column inreplicate_exdf
of the same name; and a numeric value means that the parameter will be set to that value. For example,fit_options = list(alpha_g = 0, Vcmax_at_25 = 'fit', Tp = 'column')
means thatalpha_g
will be set to 0,Vcmax_at_25
will be fit, andTp
will be set to the values in theTp
column ofreplicate_exdf
.- sd_A
The standard deviation of the measured values of the net CO2 assimilation rate, expressed in units of
micromol m^(-2) s^(-1)
. Ifsd_A
is not a number, then there must be a column inexdf_obj
calledsd_A
with appropriate units. A numeric value supplied here will overwrite the values in thesd_A
column ofexdf_obj
if it exists.- atp_use
The number of ATP molecules used per C3 cycle.
- nadph_use
The number of NADPH molecules used per C3 cycle.
- curvature_cj
A dimensionless quadratic curvature parameter greater than or equal to 0 and less than or equal to 1 that sets the degree of co-limitation between
Wc
andWj
. A value of 1 indicates no co-limitation.- curvature_cjp
A dimensionless quadratic curvature parameter greater than or equal to 0 and less than or equal to 1 that sets the degree of co-limitation between
Wcj
andWp
. A value of 1 indicates no co-limitation.- a_column_name
The name of the column in
replicate_exdf
that contains the net assimilation inmicromol m^(-2) s^(-1)
.- ci_column_name
The name of the column in
replicate_exdf
that contains the intercellular CO2 concentration inmicromol mol^(-1)
.- j_norm_column_name
The name of the column in
replicate_exdf
that contains the normalizedJ
values (with units ofnormalized to J at 25 degrees C
).- kc_column_name
The name of the column in
replicate_exdf
that contains the Michaelis-Menten constant for rubisco carboxylation inmicromol mol^(-1)
.- ko_column_name
The name of the column in
replicate_exdf
that contains the Michaelis-Menten constant for rubisco oxygenation inmmol mol^(-1)
.- oxygen_column_name
The name of the column in
exdf_obj
that contains the concentration of O2 in the ambient air, expressed as a percentage (commonly 21% or 2%); the units must bepercent
.- phips2_column_name
The name of the column in
replicate_exdf
that contains values of the operating efficiency of photosystem II (dimensionless).- qin_column_name
The name of the column in
replicate_exdf
that contains values of the incident photosynthetically active flux density inmicromol m^(-2) s^(-1)
.- rl_norm_column_name
The name of the column in
replicate_exdf
that contains the normalizedRL
values (with units ofnormalized to RL at 25 degrees C
).- total_pressure_column_name
The name of the column in
replicate_exdf
that contains the total pressure inbar
.- vcmax_norm_column_name
The name of the column in
replicate_exdf
that contains the normalizedVcmax
values (with units ofnormalized to Vcmax at 25 degrees C
).- cj_crossover_min
The minimum value of
Cc
(in ppm) whereAj
is allowed to become the overall rate-limiting factor. Ifcj_crossover_min
is set toNA
, this restriction will not be applied.- cj_crossover_max
The maximim value of
Cc
(in ppm) whereWj
is allowed to be smaller thanWc
. Ifcj_crossover_max
is set toNA
, this restriction will not be applied.- hard_constraints
To be passed to
calculate_c3_assimilation
andcalculate_c3_variable_j
; see those functions for more details.- require_positive_gmc
A character string specifying the method to be used for requiring positive values of mesophyll conductance. Can be
'none'
,'all'
, or'positive_a'
. See below for more details.- gmc_max
The maximum value of mesophyll conductance that should be considered to be acceptable. See below for more details.
Details
When fitting A-Ci + chlorophyll fluorescence curves using the Variable J
method, it is necessary to define a function that calculates the likelihood
of a given set of alpha_g
, Gamma_star
, J_at_25
,
RL_at_25
, tau
, Tp
, and Vcmax_at_25
values by
comparing a model prediction to a measured curve. This function will be passed
to an optimization algorithm which will determine the values that produce the
smallest error.
The error_function_c3_variable_j
returns such a function, which is
based on a particular replicate and a set of fitting options. It is possible
to just fit a subset of the available fitting parameters; by default, the
fitting parameters are J_at_25
, RL_at_25
, Tp
,
tau
, and Vcmax_at_25
. This behavior can be changed via the
fit_options
argument.
For practical reasons, the function actually returns values of -ln(L)
,
where L
is the likelihood. The logarithm of L
is simpler to
calculate than L
itself, and the minus sign converts the problem from
a maximization to a minimization, which is important because most optimizers
are designed to minimize a value.
Sometimes an optimizer will choose biologically unreasonable parameter values that nevertheless produce good fits to the supplied assimilation values. There are several options for preventing an optimizer from choosing such parameter values:
Enforcing Rubisco limitations
: A common problem is that the fit result may not indicate Rubisc-limited assimilation at low CO2 values, which should be the case for any A-Ci curves measured at saturating light. In this case, the optionalcj_crossover_min
andcj_crossover_max
can be used to constrain the range ofCc
values (in ppm) whereWj
is allowed to be the overall rate limiting factor. If the crossover from Rubisco-limited to RuBP-regeneration limited carboxylation occurs outside these bounds (when they are supplied), a heavy penalty will be added to the error function, preventing the optimizer from choosing those parameter values.Requiring positive gmc
: The Variable J method sometimes predicts negative values of the mesophyll conductance (gmc
). Such values are probably not realistic, especially whenCc
is above the CO2 compensation point. Therequire_positive_gmc
input argument can be used to penalize negative values ofgmc
. Whenrequire_positive_gmc
is'all'
, a heavy penalty will be added to the error function if there are any negativegmc
values. Whenrequire_positive_gmc
is'positive_a'
, a heavy penalty will be added to the error function if there are any negativegmc
values whenA
is positive; negativegmc
for negativeA
will be allowed. Whenrequire_positive_gmc
is'none'
, these restrictions are disabled and no penalties are added for negativegmc
.Preventing large values of gmc
: The Variable J method sometimes produces unreasonably high values ofgmc
. Thegmc_max
argument can be used to address this. If any predictedgmc
values exceedgmc_max
whenA
is positive, a heavy penalty will be added to the error function.
A penalty is also added for any parameter combination where An
is not a
number, or where calculate_c3_variable_j
or
calculate_c3_assimilation
produce an error.
Value
A function with one input argument guess
, which should be a numeric
vector representing values of the parameters to be fitted (which are specified
by the fit_options
input argument.) Each element of guess
is the
value of one parameter (arranged in alphabetical order.) For example, with the
default settings, guess
should contain values of J_at_25
,
RL_at_25
, tau
, Tp
, and Vcmax_at_25
(in that
order).
Examples
# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
PhotoGEA_example_file_path('c3_aci_1.xlsx')
)
# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )
# Organize the data
licor_file <- organize_response_curve_data(
licor_file,
'species_plot',
c(9, 10, 16),
'CO2_r_sp'
)
# Calculate the total pressure in the Licor chamber
licor_file <- calculate_total_pressure(licor_file)
# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_arrhenius(licor_file, c3_arrhenius_bernacchi)
# Define an error function for one curve from the set
error_fcn <- error_function_c3_variable_j(
licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE]
)
# Evaluate the error for J_at_25 = 200, RL_at_25 = 1, tau = 0.5, Tp = 15,
# Vcmax_at_25 = 100
error_fcn(c(200, 1, 0.5, 15, 100))
#> Warning: number of items to replace is not a multiple of replacement length
#> [1] 10580.61
# Make a plot of error vs. Tp when the other parameters are fixed to the values
# above. As Tp increases, it eventually stops limiting the assimilation rate
# and its value stops influencing the error.
tpu_error_fcn <- function(Tp) {error_fcn(c(200, 1, 0.5, Tp, 100))}
tpu_seq <- seq(5, 25)
lattice::xyplot(
sapply(tpu_seq, tpu_error_fcn) ~ tpu_seq,
type = 'b',
xlab = 'Tp (micromol / m^2 / s)',
ylab = 'Negative log likelihood (dimensionless)'
)
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length
#> Warning: number of items to replace is not a multiple of replacement length