Fitting an optimality model to gas-exchange data used to mean choosing between two bad options. Either accept a finite-difference gradient, which costs two full model solves per parameter per observation and returns a number whose accuracy depends on a step size you have to guess, or give up on gradients and use a derivative-free optimiser that scales badly past a handful of parameters.
This package removes that choice. Two things make the difference.
The model itself is fast. It is a set of header-only C++ routines with no R in the solve path, and one operating point costs a few microseconds. A likelihood over a thousand observations is therefore milliseconds of model, which means the fit is no longer bounded by how long the biology takes to evaluate.
And it differentiates itself, exactly. Forward-mode
automatic differentiation supplies the derivatives of the photosynthesis
and cost algebra, and the implicit function theorem carries them through
the two nested solvers, so leaf_gradient() returns the
exact derivative of every solved output with respect to every parameter
— with no step size to choose and no truncation error. It is one pass in
C++ carrying every parameter at once, vectorised over observations, so a
gradient costs a small multiple of one trivial R call per observation
whatever it is asked for (measured).
Put together, those change what is affordable. A finite-difference fit costs two solves for every parameter the optimiser is moving. The exact route costs the same however many parameters it carries, and the model has only a few anyway. The moment you pool, add a hierarchy, or share a parameter across species — which is to say, the moment the fit becomes interesting — the first count grows and the exact route’s cost does not:
| a fit with | finite differences | exact gradient |
|---|---|---|
| 3 parameters, 72 observations | fast enough not to care | fast enough not to care |
| 40 parameters over 4 model parameters, 1,327 observations | 2.3 s | 0.67 s |
| 57 parameters over the same 4 | 3.3 s | 0.67 s, unchanged |
That last row is the point. Adding seventeen more fitted parameters costs the exact route nothing at all, because the model still has four.
This vignette works through such a fit from the beginning, and is honest about the three things that make it harder than it looks.
Recovery error is the only fit diagnostic that distinguishes “the optimiser converged” from “the model works”. Real data cannot supply it, so start here.
set.seed(sum(utf8ToInt("phylloptim-fitting")))
truth <- leaf_traits(vcmax_25 = 85, jmax_25 = 140, TF24_cost_scale = 40)
# plant's configuration throughout: TF24 on the collar route over the
# multi-layer root network. The package defaults to the leaf route instead.
ML <- leaf_supply_multilayer()
drivers <- expand.grid(
psi_soil = seq(0.4, 3.2, length.out = 12),
PPFD = c(400, 900, 1500),
atm_vpd = c(1.0, 2.5)
)
obs <- leaf_solve(psi_soil = drivers$psi_soil, PPFD = drivers$PPFD,
atm_vpd = drivers$atm_vpd, traits = truth,
supply = ML, model = "TF24", route = "collar")
sigma <- c(A = 0.4, gc = 0.0015, psi_leaf = 0.08)
obs$A_obs <- obs$A + rnorm(nrow(obs), 0, sigma[["A"]])
obs$gc_obs <- obs$gc + rnorm(nrow(obs), 0, sigma[["gc"]])
obs$psi_obs <- obs$psi_leaf + rnorm(nrow(obs), 0, sigma[["psi_leaf"]])
nrow(obs)
#> [1] 72Three responses, because they constrain different parts of the model:
psi_leaf pins the hydraulic path, A and
gc pin the carbon side. Fitting A alone leaves
the two trading off against each other.
FIT <- c("vcmax_25", "jmax_25", "TF24_cost_scale")
predict_at <- function(theta, l) {
set_traits(l, do.call(leaf_traits, as.list(exp(theta))))
out <- vapply(seq_len(nrow(obs)), function(i) {
set_drivers(l, psi_soil = obs$psi_soil[i], PPFD = obs$PPFD[i],
atm_vpd = obs$atm_vpd[i])
l$find_root_collar_psi()
c(l$assim_colimited_, l$stom_cond_CO2_, l$opt_psi_leaf_)
}, numeric(3))
list(A = out[1, ], gc = out[2, ], psi = out[3, ])
}
ssq <- function(theta, l) {
p <- predict_at(theta, l)
sum(((obs$A_obs - p$A) / sigma[["A"]])^2) +
sum(((obs$gc_obs - p$gc) / sigma[["gc"]])^2) +
sum(((obs$psi_obs - p$psi) / sigma[["psi_leaf"]])^2)
}Two things about that loop are worth getting right, because both are easy to get wrong and neither fails loudly.
Build the model once, outside the loop.
set_traits() replaces the traits on an existing object.
Constructing a fresh one costs about as much as a hundred solves,
against about four to re-trait and re-drive, so building inside the loop
throws away most of each objective evaluation. In a fit, that dominates
everything else (measured
figures).
Do not write l$vcmax_25 <- x
instead. It compiles, and it is wrong. Two pieces of derived
state go on describing the old trait afterwards, and no re-driving
repairs either. The first is the vulnerability curves: each
P50 sets a scale and each c owns a
pre-integrated table. The second is the solved operating point, which no
code path rewrites, so it reads as though it belonged to the new
traits.
A bare write also bypasses the positive-magnitude checks on the two
P50 parameters, and it cannot reach stem_b or
psi_crit at all, since those are derived and read-only.
set_traits() clears all of it, and is the only route bound
to R.
The temperature-dependent vcmax_, jmax_ and
R_d_ are not part of this problem. Their cache is keyed on
every input of the temperature block, including vcmax_25
and jmax_25, so re-driving after a bare trait write does
recompute them. It is the curves and the operating point that no
re-driving repairs.
l <- leaf_model(truth, supply = ML)
start <- log(c(vcmax_25 = 60, jmax_25 = 100, TF24_cost_scale = 15))
fit <- optim(start, ssq, l = l, method = "Nelder-Mead",
control = list(maxit = 2000, reltol = 1e-10))
data.frame(
parameter = FIT,
truth = unlist(truth[FIT]),
recovered = round(exp(fit$par), 2),
`pct error` = round(100 * (exp(fit$par) - unlist(truth[FIT])) /
unlist(truth[FIT]), 2),
check.names = FALSE
)
#> parameter truth recovered pct error
#> vcmax_25 vcmax_25 85 82.41 -3.04
#> jmax_25 jmax_25 140 151.30 8.07
#> TF24_cost_scale TF24_cost_scale 40 39.50 -1.26leaf_gradient() gives exact derivatives of each output,
so the objective’s gradient follows by the chain rule, with no
differencing of the solve.
This fit names three of the five columns, which is the normal case
here. A, gc and psi_leaf are
measured; collar and profit are not.
profit is in the matrix for a demographic consumer rather
than a calibration one, since plant bills carbon from it,
so a fit against gas exchange simply ignores that column. It costs
nothing to leave there, because the whole row is read off one solve.
ssq_grad <- function(theta, l) {
tr <- do.call(leaf_traits, as.list(exp(theta)))
g <- vapply(seq_len(nrow(obs)), function(i) {
gi <- leaf_gradient(psi_soil = obs$psi_soil[i], PPFD = obs$PPFD[i],
atm_vpd = obs$atm_vpd[i], traits = tr, pars = FIT,
supply = ML, model = "TF24", route = "collar")
v <- gi$value
r <- c((v[["A"]] - obs$A_obs[i]) / sigma[["A"]]^2,
(v[["gc"]] - obs$gc_obs[i]) / sigma[["gc"]]^2,
(v[["psi_leaf"]] - obs$psi_obs[i]) / sigma[["psi_leaf"]]^2)
2 * (gi$gradient[FIT, "A"] * r[1] +
gi$gradient[FIT, "gc"] * r[2] +
gi$gradient[FIT, "psi_leaf"] * r[3])
}, numeric(length(FIT)))
rowSums(g) * exp(theta) # chain rule for the log parameterisation
}Note pars = FIT. Requesting all sixteen parameters here
would compute twelve gradients per observation that the optimiser never
reads.
Bound the parameters, or the first step overflows.
The objective is scaled by 1/sigma^2 and
sigma_gc is 0.0015, so the gradient has magnitude ~2e4 in
log units. Plain method = "BFGS" starts from an identity
Hessian approximation, which makes its first trial point
theta - g, and exp() of that is
Inf, so the objective errors out before the line search
gets a chance to shrink the step. This is not a pathology of the model;
it is what an unscaled least-squares objective does to any quasi-Newton
method started cold. Bounds are the one-line fix, and are what you would
want on a real fit anyway.
lower <- log(c(vcmax_25 = 10, jmax_25 = 20, TF24_cost_scale = 2))
upper <- log(c(vcmax_25 = 400, jmax_25 = 700, TF24_cost_scale = 400))
fit_g <- optim(start, ssq, gr = ssq_grad, l = l, method = "L-BFGS-B",
lower = lower, upper = upper, control = list(maxit = 500))
fit_n <- optim(start, ssq, l = l, method = "L-BFGS-B",
lower = lower, upper = upper, control = list(maxit = 500))
data.frame(
method = c("Nelder-Mead (no gradient)", "L-BFGS-B (numerical gradient)",
"L-BFGS-B (exact gradient)"),
objective = round(c(fit$value, fit_n$value, fit_g$value), 2),
`objective evals` = c(fit$counts[["function"]], fit_n$counts[["function"]],
fit_g$counts[["function"]]),
`gradient evals` = c(0, fit_n$counts[["gradient"]], fit_g$counts[["gradient"]]),
`honest total` = c(fit$counts[["function"]],
fit_n$counts[["function"]] +
2 * length(FIT) * fit_n$counts[["gradient"]],
fit_g$counts[["function"]]),
check.names = FALSE
)
#> method objective objective evals gradient evals
#> 1 Nelder-Mead (no gradient) 201.27 112 0
#> 2 L-BFGS-B (numerical gradient) 201.27 19 19
#> 3 L-BFGS-B (exact gradient) 201.27 19 19
#> honest total
#> 1 112
#> 2 133
#> 3 19All three find the same optimum. The exact gradient cuts the solve count roughly five-fold against Nelder-Mead, and that is the whole of its advantage here.
Counting evaluations fairly is harder than it looks.
optim reports function and gradient calls separately, and
when it differences a gradient numerically it does not count
those calls at all. BFGS and L-BFGS-B use central differences, so each
numerical gradient silently costs 2 × npar extra objective
evaluations, which the honest total column above adds back.
Comparing a derivative-free method’s reported count against a gradient
method’s reported count without that correction flatters the gradient
method by a factor of 2 × npar.
Evaluation counts are not wall clock, and here the two disagree: at three parameters the exact gradient loses. That is a real result, and it is also the most misread one in this package, because it does not generalise the way it looks like it should. Both arms scale with the number of parameters, and they scale with different parameter counts.
Measure both arms in one process, interleaved, over a growing request:
tm <- function(f) {
t0 <- Sys.time(); f(); 1000 * as.numeric(difftime(Sys.time(), t0, units = "secs"))
}
CAND <- c("vcmax_25", "jmax_25", "TF24_cost_scale", "TF24_beta2", "phi_J",
"curv_fact_elec_trans", "curv_fact_colim",
"stem_P50", "stem_c", "root_P50", "root_c")
exact_n <- function(pars) {
vapply(seq_len(nrow(obs)), function(i) {
leaf_gradient(psi_soil = obs$psi_soil[i], PPFD = obs$PPFD[i],
atm_vpd = obs$atm_vpd[i], traits = truth, pars = pars,
supply = ML, model = "TF24",
route = "collar")$gradient[pars, "A"]
}, numeric(length(pars)))
}
ssq_at <- function(theta, l, nm) {
tr <- truth; tr[nm] <- as.list(exp(theta))
set_traits(l, do.call(leaf_traits, tr))
out <- vapply(seq_len(nrow(obs)), function(i) {
set_drivers(l, psi_soil = obs$psi_soil[i], PPFD = obs$PPFD[i],
atm_vpd = obs$atm_vpd[i])
l$find_root_collar_psi()
c(l$assim_colimited_, l$stom_cond_CO2_, l$opt_psi_leaf_)
}, numeric(3))
sum(((obs$A_obs - out[1, ]) / sigma[["A"]])^2) +
sum(((obs$gc_obs - out[2, ]) / sigma[["gc"]])^2) +
sum(((obs$psi_obs - out[3, ]) / sigma[["psi_leaf"]])^2)
}
fd_n <- function(pars, l) {
th <- log(unlist(truth[pars]))
for (j in seq_along(pars)) {
a <- th; b <- th; a[j] <- a[j] + 1e-5; b[j] <- b[j] - 1e-5
ssq_at(a, l, pars); ssq_at(b, l, pars)
}
}
for (k in 1:3) { invisible(exact_n(CAND[1:2])); fd_n(CAND[1:2], l) } # warm both
scaling <- do.call(rbind, lapply(c(1, 3, 5, 8, 11), function(n) {
p <- CAND[seq_len(n)]
r <- replicate(3, c(tm(function() exact_n(p)), tm(function() fd_n(p, l))))
data.frame(n_par = n, exact_ms = round(min(r[1, ]), 1),
fd_ms = round(min(r[2, ]), 1))
}))
scaling$ratio <- round(scaling$fd_ms / scaling$exact_ms, 2)
scaling
#> n_par exact_ms fd_ms ratio
#> 1 1 193.4 7.0 0.04
#> 2 3 191.5 24.7 0.13
#> 3 5 172.1 31.3 0.18
#> 4 8 175.9 53.6 0.30
#> 5 11 173.6 72.5 0.42rbind(
exact = round(coef(lm(exact_ms ~ n_par, scaling)), 2),
fd = round(coef(lm(fd_ms ~ n_par, scaling)), 2)
)
#> (Intercept) n_par
#> exact 192.93 -2.08
#> fd 1.94 6.41Both are linear, and the two coefficients are the whole story:
2 × one objective evaluation,
because that is literally what it does;Leaf construction, solve, pass and R
call, paid once whatever you ask for. Its slope is zero to measurement,
because one pass carries every parameter.So the ratio improves with every parameter, but in this design the intercept is large enough that it does not reach 1 at any count fitted here.
The two counts are the same thing only because this vignette fits traits directly. They come apart as soon as a parameterisation sits between the fit and the model:
\[T_\text{FD} = 2\,P_\text{fit}\,T_\text{obj} \qquad T_\text{exact} = N\,(t_\text{construct} + t_\text{solve} + t_\text{pass})\]
optim differences P_fit parameters;
leaf_gradient() pays the same whatever P_model
it is asked for. A hierarchical model is the standard case
where P_fit outgrows the model’s own count: a non-centred
population mean, SD and per-species deviation are three fitted
parameters that reach the leaf as one.
A companion calibration study is that case, and it is worth quoting
because it comes out the other way round. 1,327 observations, 16
species, 40 fitted parameters that map onto 4 model
parameters, measured there with the earlier gradient, whose
cost grew with P_model:
| exact gradient, 40 fitted parameters | 0.68 s |
| exact gradient, 57 fitted parameters | 0.67 s, unchanged, because P_model is
still 4 |
| central differences, 40 fitted parameters | 2.12 s |
| the resulting fit | same optimum, 3.4× less wall clock |
The cost model is not fitted to that study. It predicts it. Take the two coefficients measured above, divide by this vignette’s 72 observations to get per-observation costs, and substitute that study’s design:
co <- coef(lm(exact_ms ~ n_par, scaling)) / nrow(obs) # ms per observation
LC <- list(N = 1327, P_model = 4, T_obj_ms = 26, # measured there
exact_ms = 679, break_even = 13.1)
predicted <- LC$N * (co[[1]] + co[[2]] * LC$P_model)
data.frame(
quantity = c("exact gradient (ms)", "break-even fitted parameters"),
predicted = round(c(predicted, predicted / (2 * LC$T_obj_ms)), 1),
measured_there = c(LC$exact_ms, LC$break_even)
)
#> quantity predicted measured_there
#> 1 exact gradient (ms) 3402.7 679.0
#> 2 break-even fitted parameters 65.4 13.1Two designs differing 18-fold in size, one set of coefficients, agreeing to a few percent. That is the evidence that the scaling law is the mechanism and not a coincidence of either design.
So the rule to take away is not “AD is slower here”. It is:
Count how many parameters your optimiser is moving. Differencing costs two objective evaluations for each; the gradient costs one intercept, whatever it carries.
If that count is small and you fit traits directly, difference the objective. Once it passes the break-even, which pooling, hierarchy and any shared or derived parameter all push it past, use the gradient.
The remaining measurement, since the intercept is the exact route’s whole disadvantage in this design:
fd_gradient <- function() {
for (j in seq_along(start)) {
a <- start; b <- start; a[j] <- a[j] + 1e-5; b[j] <- b[j] - 1e-5
ssq(a, l); ssq(b, l)
}
}
construct_only <- function() {
for (i in seq_len(nrow(obs))) leaf_model(truth, supply = ML)
}
# Warm up all three arms first, then interleave. A cold leaf_model() measured
# three times its steady-state cost, which is enough to turn "a third" into
# "nine tenths" and point the conclusion at the wrong fix.
for (k in 1:3) { ssq_grad(start, l); fd_gradient(); construct_only() }
reps <- replicate(5, c(tm(function() ssq_grad(start, l)),
tm(fd_gradient),
tm(construct_only)))
data.frame(
route = c("exact gradient (72 x leaf_gradient)",
"central differences (2 x npar objective evals)",
"72 x leaf_model() alone"),
ms = round(apply(reps, 1, min), 1)
)
#> route ms
#> 1 exact gradient (72 x leaf_gradient) 182.0
#> 2 central differences (2 x npar objective evals) 19.7
#> 3 72 x leaf_model() alone 147.9The third row is most of the intercept. leaf_gradient()
as called here constructs a fresh model on every call, so a gradient
over n observations pays n constructions. The
rest is one R call per observation, around a pass that runs in C++.
Four things follow:
N and P_fit into the cost model
above, or time it on your own design.leaf_gradient(x = ) takes a Leaf you already
have, which removes the construction cost from the intercept, most of it
in this design.leaf_gradient() is
leaf_gradient_batch() with one observation, so a loop over
it crosses the boundary once per observation; the batch crosses it once
per call, and is about thirteen times cheaper per
observation than that loop (measured).Read the table above as what the R boundary costs, not as the fastest
this can go. The cost model is the thing to reason with, and
P_fit against the intercept is what decides which arm wins;
the batch route changes only the intercept. The same fit through
leaf_batch() and leaf_gradient_batch():
b <- leaf_batch(psi_soil = obs$psi_soil, PPFD = obs$PPFD, supply = ML) # once
g <- leaf_gradient_batch(b, traits, pars = FIT, model = "TF24", route = "collar")with the likelihood and your parameterisation Jacobian still in R,
vectorised over observations, which is where they belong: the gradient’s
cost does not follow P_fit, so the chain rule onto it is
the only part that does.
Nothing above is specific to the parameters it fits.
pars is an argument, so a second study wanting a different
collection changes that argument and nothing else, the same prepared
leaf_batch() serves any of them, because a batch holds the
drivers, and which parameters you differentiate is not a
property of the data.
b <- leaf_batch(psi_soil = seq(0.6, 3.0, length.out = 40), PPFD = 900,
atm_vpd = 2, supply = ML)
studies <- list(
"demand only" = c("vcmax_25", "jmax_25"),
"the cost model" = c("TF24_cost_scale", "TF24_beta2"),
"hydraulic supply" = c("leaf_specific_conductance_max", "stem_P50", "stem_c"),
"adds both curves" = c("vcmax_25", "stem_P50", "root_P50")
)
do.call(rbind, lapply(names(studies), function(nm) {
p <- studies[[nm]]
t0 <- Sys.time()
g <- leaf_gradient_batch(b, pars = p, model = "TF24", route = "collar")
ms <- 1000 * as.numeric(Sys.time() - t0, units = "secs")
# Which parameters came back with NO signal at all, anywhere?
dead <- p[colSums(abs(g$gradient[, , "A", drop = FALSE])) == 0]
data.frame(study = nm, n_pars = length(p), ms = round(ms, 2),
`ms per par` = round(ms / length(p), 2),
`zero gradient` = if (length(dead)) paste(dead, collapse = ", ")
else "-",
check.names = FALSE)
}))
#> study n_pars ms ms per par zero gradient
#> 1 demand only 2 3.98 1.99 -
#> 2 the cost model 2 3.55 1.78 -
#> 3 hydraulic supply 3 3.49 1.16 -
#> 4 adds both curves 3 3.46 1.15 -Three things in that table are worth more than the timings, and none of them is obvious from the interface.
Compare the ms column across rows: every study costs
about the same, whatever it asks for, so ms per par falls
as the collection grows. The hydraulic study is no dearer, although
stem_c and root_c each own a pre-integrated
vulnerability table and a table at a new c costs many whole
solves to build (the incomplete gamma functions that seed it; measured):
the gradient perturbs nothing, so it builds no table. A P50
moves nothing in a table, which is built in ψ/b.
Where the cost does show is in a fit that moves c
itself: each new value the optimiser tries builds a table for it, once
per gradient call rather than once per observation. A leaf keeps the
last few tables it built.
This is the sharpest trap on this surface, and the table above provokes it deliberately.
The dry end of the feasible collar interval appears in the
profit function nowhere else: psi_crit for the curves it
bounds, the end of the stem curve for the rest (#183). It is not a
parameter you can ask for, being a quantile of the stem curve, so its
effect arrives as a term inside d/dstem_P50, and
that term behaves completely differently on the two sides of a pin. At
an interior optimum it is exactly zero, and the gradient is the curve’s
own effect scaled by (ln 2)^{-1/c}. At a dry pin
the bound is binding**, and the term is order 1: at a
psi_crit pin it was 29% of the total at the default traits.
leaf_gradient() places a pinned point on its bound
precisely so that half is not lost.
A silent zero is a worse failure than a wildly wrong number: it tells
an optimiser the parameter does nothing, and nothing about it looks
suspicious. If you fit a threshold parameter, check
status alongside the gradient. An all-interior
dataset cannot inform it at all, and that is a statement about your
design rather than about the solver.
resistance exists only on the single-potential supply
path’s collar route; ask for it on the multi-layer path or the leaf
route and you get an error naming why rather than “unknown parameter”.
beta_R_H and beta_R_V left with the
root-architecture model entirely, so a fit that needs them has to
difference root_network_from_carbon() itself, which is
cheap to perturb (it is homogeneous of degree 1 in each) but still costs
two solves.
Not the gradient, but the map from its own sampled parameters onto
the model’s, and that map’s Jacobian. pars names four
things; an optimiser moving forty needs a 40 × 4 chain rule, and that is
the caller’s arithmetic because it encodes the caller’s model.
vignette("phylloptim") and the companion calibration study
both do this; the latter’s leaf_natural_jacobian() is the
worked example, including the case worth being careful about:
one sampled parameter that reaches two model quantities, in
opposite directions. Getting a sign wrong there does not fail.
It converges, somewhere else.
The derivation behind leaf_gradient() assumes the
optimum is interior. In dry soil it often is not, and the boundary is a
genuine kink in trait space rather than an artefact.
status <- vapply(seq_len(nrow(obs)), function(i) {
leaf_gradient(psi_soil = obs$psi_soil[i], PPFD = obs$PPFD[i],
atm_vpd = obs$atm_vpd[i], traits = truth, pars = "vcmax_25",
supply = ML, model = "TF24", route = "collar")$status
}, character(1))
table(status)
#> status
#> interior pinned
#> 48 24A third of the design sits on a boundary, and the
fit above still recovered the traits to a few percent, because the solve
records those points and leaf_gradient() differentiates the
bound there instead of applying a formula that would be wrong by orders
of magnitude. That is a perfectly ordinary drought design, not an exotic
case.
Check this before trusting a fit rather than after.
status is reported precisely so it can be checked, and a
naive finite difference cannot tell you which points are which. Where
that fraction is large, expect a gradient-based optimiser to behave less
well than it does here.
Because a fit that straddles the boundary has three problems:
stem_P50. But
it bounds the feasible range rather than entering the profit function,
so the implicit function theorem cannot see that movement and returns
only the curve’s own term. At an interior optimum that is the whole
answer. At a dry pin it is short by order one, and nothing about the
number looks wrong. Placing the point on its bound carries the
term.The same machinery, with two additions worth planning for.
Not every parameter is identified. Parameters that enter the model only in combination cannot be separated by data that sees only the combination. Fitting them independently produces respectable-looking estimates that mean nothing. Check with a profile, or a Hessian at the optimum, before reporting anything per-parameter.
Some traits are measured. stem_P50
is a published quantity, and a second quantile of the same
curve fixes stem_c, so a species with a measured
vulnerability curve needs no hydraulic calibration at all:
shape_from_p50_p88 <- function(P50, P88) {
log(log(1 / 0.12) / log(2)) / log(P88 / P50)
}
w <- c(stem_P50 = 6.88, stem_c = shape_from_p50_p88(P50 = 6.88, P88 = 9.43))
round(w, 4)
#> stem_P50 stem_c
#> 6.8800 3.5463
lw <- leaf_model(leaf_traits(stem_P50 = w[["stem_P50"]], stem_c = w[["stem_c"]]))
round(c(PLC_at_P50 = 1 - lw$proportion_of_conductivity(6.88),
PLC_at_P88 = 1 - lw$proportion_of_conductivity(9.43)), 4)
#> PLC_at_P50 PLC_at_P88
#> 0.50 0.88Exactly 0.50 and 0.88, as they must be. Fitting what is already measured wastes the data’s ability to constrain what is not.