Fits a hyperbolic C4 assimilation model to an experimental curve

Fits a model to an experimentally measured C4 CO2 response curve using the data in the exdf object along with a few other user-supplied parameters. This function can accomodate alternative column names for the variables taken from the Licor file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

fit_c4_aci_hyperbola(
    replicate_exdf,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    sd_A = 'RMSE',
    OPTIM_FUN = optimizer_nmkb(1e-7),
    lower = list(),
    upper = list(),
    fit_options = list(),
    error_threshold_factor = 0.147,
    hard_constraints = 0,
    calculate_confidence_intervals = TRUE
  )

Arguments

replicate_exdf: An exdf object representing one CO2 response curve.
a_column_name: The name of the column in replicate_exdf that contains the net assimilation in micromol m^(-2) s^(-1).
ci_column_name: The name of the column in replicate_exdf that contains the intercellular CO2 concentration in micromol mol^(-1).
sd_A: A value of the standard deviation of measured A values, or the name of a method for determining the deviation; currently, the only supported option is 'RMSE'.
OPTIM_FUN: An optimization function that accepts the following input arguments: an initial guess, an error function, lower bounds, and upper bounds. It should return a list with the following elements: par, convergence, value, and (optionally) message. See optimizers for a list of available options.
lower: A list of named numeric elements representing lower bounds to use when fitting. Values supplied here override the default values (see details below). For example, lower = list(Vmax = 10) sets the lower limit for Vmax to 10 micromol / m^2 / s.
upper: A list of named numeric elements representing upper bounds to use when fitting. Values supplied here override the default values (see details below). For example, upper = list(Vmax = 200) sets the upper limit for Vmax to 200 micromol / m^2 / s.
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 in exdf_obj of the same name; and a numeric value means that the parameter will be set to that value. For example, fit_options = list(rL = 0, Vmax = 'fit', c4_curvature = 'column') means that rL will be set to 0, Vmax will be fit, and c4_curvature will be set to the values in the c4_curvature column of replicate_exdf.
error_threshold_factor: To be passed to confidence_intervals_c4_aci_hyperbola when calculate_confidence_intervals is TRUE.
hard_constraints: To be passed to calculate_c4_assimilation_hyperbola; see that function for more details.
calculate_confidence_intervals: A logical value indicating whether or not to estimate confidence intervals for the fitting parameters using confidence_intervals_c4_aci_hyperbola.

Details

This function calls calculate_c4_assimilation_hyperbola to calculate values of net assimilation. The user-supplied optimization function is used to vary the values of c4_curvature, c4_slope, rL, and Vmax to find ones that best reproduce the experimentally measured values of net assimilation. By default, the following options are used for the fits:

c4_curvature: lower = -10, upper = 10, fit_option = 'fit'
c4_slope: lower = -50, upper = 1000, fit_option = 'fit'
rL: lower = -10, upper = 100, fit_option = 'fit'
Vmax: lower = -50, upper = 1000, fit_option = 'fit'

With these settings, all of the parameters are fit during the process (see fit_options above). The bounds are chosen liberally to avoid any bias.

An initial guess for the parameters is generated by calling initial_guess_c4_aci_hyperbola. Note that any fixed values specified in the fit options will override the values returned by the guessing function.

The fit is made by creating an error function using error_function_c4_aci_hyperbola and minimizing its value using OPTIM_FUN, starting from the initial guess described above. The optimizer_nmkb optimizer is used by default since it has been found to reliably return great fits. However, it is a fast optimizer that can get stuck in local minima. If it seems to be returning bad fits, consider using the optimizer_deoptim optimizer instead, but be aware that the fits will take more time to complete.

Unlike the model represented by calculate_c4_assimilation, the model in calculate_c4_assimilation_hyperbola is smooth in the sense that small changes in the input parameters cause small changes in its outputs. Because of this, it is a fairly easy model to fit.

This function assumes that replicate_exdf represents a single C4 A-Ci curve. To fit multiple curves at once, this function is often used along with by.exdf and consolidate.

Value

A list with two elements:

fits: An exdf object including the original contents of replicate_exdf along with several new columns:
- The fitted values of net assimilation will be stored in a column whose name is determined by appending '_fit' to the end of a_column_name; typically, this will be 'A_fit'.
- Residuals (measured - fitted) will be stored in a column whose name is determined by appending '_residuals' to the end of a_column_name; typically, this will be 'A_residuals'.
- Values of fitting parameters will be stored in the c4_curvature, c4_slope, rL, and Vmax columns.
- The other outputs from calculate_c4_assimilation_hyperbola will be stored in columns with the usual names: Ag, Ainitial, Amax, An, c4_curvature, c4_slope, rL, Vinitial, Vmax, and c4_assimilation_hyperbola_msg.
fits_interpolated: An exdf object including the calculated assimilation rates at a fine spacing of Ci values (step size of 1 micromol mol^(-1)).
parameters: An exdf object including the identifiers, fitting parameters, and convergence information for the A-Ci curve:
- The best-fit values are stored in the c4_curvature, c4_slope, rL, and Vmax. If calculate_confidence_intervals is TRUE, upper and lower limits for each of these parameters will also be included.
- The convergence column indicates whether the fit was successful (==0) or if the optimizer encountered a problem (!=0).
- The feval column indicates how many cost function evaluations were required while finding the optimal parameter values.
- The residual stats as returned by residual_stats are included as columns with the default names: dof, RSS, RMSE, etc.
- The Akaike information criterion is included in the AIC column.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c4_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'
)

# Fit just one curve from the data set (it is rare to do this).
one_result <- fit_c4_aci_hyperbola(
  licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE]
)

# Fit all curves in the data set (it is more common to do this)
aci_results <- consolidate(by(
  licor_file,
  licor_file[, 'species_plot'],
  fit_c4_aci_hyperbola
))

# View the fitting parameters for each species / plot
col_to_keep <- c(
  'species', 'plot',                                       # identifiers
  'c4_curvature', 'c4_slope', 'rL', 'Vmax',                # best estimates for parameter values
  'dof', 'RSS', 'MSE', 'RMSE', 'RSE',                      # residual stats
  'convergence', 'convergence_msg', 'feval', 'optimum_val' # convergence info
)

aci_results$parameters[ , col_to_keep, TRUE]
#>   species [UserDefCon] (NA) plot [UserDefCon] (NA)
#> 1                     maize                      5
#> 2                   sorghum                      2
#> 3                   sorghum                      3
#>   c4_curvature [fit_c4_aci_hyperbola] (dimensionless)
#> 1                                           0.6976615
#> 2                                           0.5614384
#> 3                                           0.5236451
#>   c4_slope [fit_c4_aci_hyperbola] (mol m^(-2) s^(-1))
#> 1                                           1.0100960
#> 2                                           0.9970135
#> 3                                           0.9282023
#>   rL [fit_c4_aci_hyperbola] (micromol m^(-2) s^(-1))
#> 1                                          1.3225714
#> 2                                          2.7695450
#> 3                                          0.7192808
#>   Vmax [fit_c4_aci_hyperbola] (micromol m^(-2) s^(-1))
#> 1                                             65.12738
#> 2                                             72.43272
#> 3                                             71.18932
#>   dof [residual_stats] (NA) RSS [residual_stats] ((micromol m^(-2) s^(-1))^2)
#> 1                         9                                         104.63148
#> 2                         9                                         220.78920
#> 3                         9                                          27.39311
#>   MSE [residual_stats] ((micromol m^(-2) s^(-1))^2)
#> 1                                          8.048575
#> 2                                         16.983784
#> 3                                          2.107162
#>   RMSE [residual_stats] (micromol m^(-2) s^(-1))
#> 1                                       2.837001
#> 2                                       4.121139
#> 3                                       1.451607
#>   RSE [residual_stats] (micromol m^(-2) s^(-1))
#> 1                                      3.409651
#> 2                                      4.952992
#> 3                                      1.744614
#>   convergence [fit_c4_aci_hyperbola] ()
#> 1                                     0
#> 2                                     0
#> 3                                     0
#>   convergence_msg [fit_c4_aci_hyperbola] () feval [fit_c4_aci_hyperbola] ()
#> 1                    Successful convergence                             318
#> 2                    Successful convergence                             263
#> 3                    Successful convergence                             359
#>   optimum_val [fit_c4_aci_hyperbola] ()
#> 1                              32.00192
#> 2                              36.85588
#> 3                              23.29093

# View the fits for each species / plot
plot_c4_aci_hyperbola_fit(aci_results, 'species_plot', ylim = c(0, 100))


# View the residuals for each species / plot
lattice::xyplot(
  A_residuals ~ Ci | species_plot,
  data = aci_results$fits$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', aci_results$fits$units$Ci, ']'),
  ylab = paste('Assimilation rate residuals [', aci_results$fits$units$A_residuals, ']')
)