---
title: "Linear Response Plateau (LRP)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Linear Response Plateau (LRP)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.5, dpi = 120)
```

```{r setup}
library(trialSizing)
```

## Theory

In a uniformity trial, adjacent basic experimental units (BEU) are grouped into
plots of increasing size `X`, and the coefficient of variation `CV(X)` is
computed among plots of each size. The CV falls as plots get larger, but the
gain shrinks: past a certain size, extra area buys almost no precision. The
Linear Response Plateau model formalizes that as two joined straight lines,

$$
CV_{(X)} =
\begin{cases}
a + bX + \varepsilon, & \text{if } X \le X_o \\
p + \varepsilon,      & \text{if } X > X_o
\end{cases}
$$

The first segment falls with slope $b < 0$; the second is flat at the plateau
$p$. The join between them is the estimate of interest: $X_o$ is the optimal
plot size and $p = CV_{Xo}$ is the CV expected at that size. From the model,

$$
X_o = \frac{p - a}{b}, \qquad CV_{Xo} = a + b X_o
$$

### How the breakpoint is estimated

Fitting a broken line is awkward for general-purpose nonlinear least squares.
The model has a kink, and the residual sum of squares (SSE) is a *stepped*
function of the breakpoint: moving $X_o$ between two observed plot sizes does
not change which points fall in each segment, so the SSE stays flat and its
derivative is zero. Gradient-based fitters such as `nls()` and `nlsLM()` stall
on those flat steps and return whichever basin their starting value happened to
land in.

`fit_lrp()` avoids this by **profiling** the breakpoint. It walks a fine grid of
candidate values for $X_o$; at each candidate the model becomes linear, so the
coefficients come from ordinary least squares. The candidate with the smallest
SSE wins. There are no starting values and no convergence failures, and the
result is the global least-squares optimum by construction.

That does not make the breakpoint unambiguous, though: the profile may have
several basins that fit almost equally well. See
[Competing breakpoints](#competing-breakpoints) below.

### References

Original method: Paranaíba, P. F., Ferreira, D. F. & Morais, A. R. (2009).
Tamanho ótimo de parcelas experimentais: proposição de métodos de estimação.
*Revista Brasileira de Biometria*, 27(2), 255-268.

The implementation is validated against published results in the package tests;
see Cargnelutti Filho, A., Loro, M. V., Ortiz, V. M. & Andretta, J. A. (2025).
Determinação do tamanho de parcela para avaliar a massa de parte aérea de
grão-de-bico. *Revista Vivências*, 21(43), 499-513.

## Data

The examples use the bundled simulated uniformity trial (`?uniformity_trial`):
an 8 × 12 grid of 1 m² BEU per trial, grouped by `calc_cv_shapes()` into plot
shapes. `X` is the plot size in BEU and `CV1`--`CV3` are the CVs of trials 1
to 3.

```{r data, message = FALSE}
grid_mat <- function(t)
  as.matrix(uniformity_trial[uniformity_trial$trial == t,
                             grep("^col", names(uniformity_trial))])

tab1 <- calc_cv_shapes(grid_mat("T1"))
X   <- tab1$x
CV1 <- tab1$cv
CV2 <- calc_cv_shapes(grid_mat("T2"))$cv
CV3 <- calc_cv_shapes(grid_mat("T3"))$cv
```

Note that `X` repeats: several shapes share the same size (four different
shapes of 6 BEU, for instance). Supply one CV per shape, never the mean per
size, or the fit loses the within-size variability it needs.

## Basic use

Every fit below passes `step = 0.01` rather than the default `0.001`. That is
only to keep this vignette quick to build, and to keep every number on this page
consistent with the others: the coarser grid already resolves $X_o$ to two
decimals, and only the third decimal differs. See
[Fine-tuning](#fine-tuning) for what `step` does.

```{r fit}
fit <- fit_lrp(x = X, cv = CV1, step = 0.01)
fit
```

$X_o \approx 9.2$ m² with $CV_{Xo} \approx 6.9\%$ for trial 1. The competing
minima listed underneath are explained
[further down](#competing-breakpoints).

`summary()` splits the model coefficients from the fit statistics:

```{r summary}
summary(fit)
```

Everything is also available directly:

```{r access}
fit$coefficients
fit$parameters["Breakpoint"]
round(fit$residuals[1:5], 3)
```

`predict()` evaluates the fitted model at any plot size, which is handy to check
the plateau behaviour:

```{r predict}
predict(fit, newx = c(1, 5, 12, 15))
```

The last two values are identical: past the breakpoint the model is flat.

## Plot

The `plot()` method draws the article-style figure. Title and styling go here,
not in `fit_lrp()`:

```{r plot}
plot(fit, title = "Trial 1")
```

For Portuguese-language figures, switch the decimal mark and the conditional
word:

```{r plot-ptbr}
plot(fit, title = "Ensaio 1", decimal_mark = ",", cond_word = "se")
```

Other useful arguments: `annotate_model = FALSE` drops the equation block,
`base_size` / `label_size` / `title_size` scale the text, `theme` accepts any
ggplot2 theme, and `xlab` / `ylab` relabel the axes.

To export, set `save = TRUE`. TIFF is written with LZW compression; PDF and EPS
are vector formats, usually the better choice for line art like this.

```{r save, eval = FALSE}
plot(fit, title = "Trial 1",
     save = TRUE, file = "trial1.tiff", format = "tiff",
     dpi = 300, width = 18, height = 12, units = "cm")
```

## Several trials at once

Pass a data frame plus the column names. One model is fitted per trial:

```{r multi}
trials <- rbind(
  data.frame(x = X, cv = CV1, trial = "Trial 1"),
  data.frame(x = X, cv = CV2, trial = "Trial 2"),
  data.frame(x = X, cv = CV3, trial = "Trial 3")
)

res <- fit_lrp(trials, x = "x", cv = "cv", trial = "trial", step = 0.01)
res$summary
```

The three breakpoints differ from trial to trial, as real trials do; their mean
is the natural overall recommendation. The `n_local` column counts competing
basins per trial.

The individual fits stay accessible, so any single trial can be inspected or
plotted as usual:

```{r multi-access}
res$fits[["Trial 2"]]
```

`plot()` on the multi-trial object arranges the panels in a grid of at most six
(3 rows × 2 columns), paginating into separate figures beyond that. It needs
the **patchwork** package.

```{r multi-plot, eval = FALSE}
plot(res, decimal_mark = ",", cond_word = "se", label_size = 3)
```

Column names default to `"x"`, `"cv"` and `"trial"`, so `fit_lrp(trials)` alone
works when the data frame already uses those names. Missing columns produce a
clear error rather than a cryptic one:

```{r error, error = TRUE}
fit_lrp(trials, x = "x", cv = "CV", trial = "trial")
```

## Competing breakpoints {#competing-breakpoints}

Uniformity trials rarely have many distinct plot sizes, so the SSE profile is
stepped and often has several local minima. The grid search always returns the
global optimum, but a second basin may fit nearly as well. When that happens the
breakpoint is **not sharply identified**, and different software will
legitimately report different values.

`fit_lrp()` makes this visible instead of hiding it. Every fit carries
`$local_minima`, listing each basin with its SSE and its excess over the
optimum. Basins fitting within `local_min_tol` (10% by default) are flagged in
the `competing` column and starred by `print()`. Nothing is signalled by a
warning: stepped profiles almost always have several minima, so a warning would
fire on nearly every fit and stop meaning anything.

```{r competing}
fit2 <- fit_lrp(X, CV2, step = 0.01)
```

```{r competing-table}
lm2 <- fit2$local_minima
lm2[c("breakpoint", "SSE", "SSE_excess")] <-
  round(lm2[c("breakpoint", "SSE", "SSE_excess")], 3)
lm2
```

Read the `SSE_excess` column as "how much worse than the optimum". When a
runner-up basin fits only a few percent worse than the reported optimum, those
two plot sizes are effectively tied and the breakpoint carries more uncertainty
than its decimals suggest. A competitor 30% worse, in contrast, can be safely
dismissed.

The whole profile is available for inspection or plotting:

```{r profile, fig.height = 3.5}
prof <- fit2$sse_profile
plot(prof$breakpoint, prof$SSE, type = "l",
     xlab = "Breakpoint", ylab = "SSE", las = 1)
abline(v = fit2$parameters["Breakpoint"], col = "forestgreen", lwd = 2)
abline(v = fit2$local_minima$breakpoint, col = "red", lty = 2)
```

### Reproducing a fit from another package

When a colleague's `nls()`-based script or another package reports a different
breakpoint, it is usually because their starting value landed in another basin.
The `start` argument reproduces that solution *without changing the reported
fit*, so you can see exactly what they got and what it cost in fit quality:

```{r compat}
soy_x  <- c(1, 2, 4, 8, 2, 4, 8, 16, 4, 8, 16, 32, 5, 10, 20, 40)
soy_cv <- c(18.699092, 14.130115, 10.321934, 7.990773, 12.690754, 9.995547,
            7.916291, 7.588785, 9.276139, 7.130636, 4.777755, 4.279987,
            8.411412, 5.818440, 3.431264, 3.335962)

seeded <- fit_lrp(soy_x, soy_cv, start = 10.17, step = 0.01)
seeded$compat
```

The reported fit is still the optimum at 5.79; `$compat` shows that a fitter
seeded near 10 would stop at 10.17 with an SSE 15% higher. That is a large
enough gap to prefer the optimum, but the number is there for you to judge, and
to cite when reconciling results across software.

### Deciding between basins

A statistically better fit is not automatically the better agronomic answer. A
very early breakpoint can be a "false plateau": the trial may lack the range of
plot sizes needed for a genuine plateau to emerge, so the model latches onto a
flat stretch that is really just noise. When two basins are close in SSE, weigh
the number of points supporting each segment, whether the plateau looks
plausible on the plot, and what plot sizes are feasible in the field.

Once you have decided, `search_range` documents the choice explicitly:

```{r search-range-basin}
fit_lrp(soy_x, soy_cv, search_range = c(8, 20), step = 0.01)$parameters["Breakpoint"]
```

## Fine-tuning

**`search_range`** restricts where the breakpoint is searched, which helps when
an outlier in the transition zone pulls the estimate somewhere implausible.
It must fall inside the data range, and a warning is issued if the breakpoint
lands on the edge of the range you gave.

```{r search-range}
fit_lrp(X, CV1, search_range = c(5, 12), step = 0.01)$parameters["Breakpoint"]
```

**`step`** is the grid resolution, 0.001 by default. Coarser grids are faster
and slightly less precise:

```{r step}
fit_lrp(X, CV1, step = 0.01)$parameters["Breakpoint"]
```

**`method`** chooses how the linear coefficients are estimated at each candidate
breakpoint. `"segment"` (default) uses only the observations with `x <= X0`,
reproducing the published procedure. `"ramp"` uses every observation on the
basis `pmin(x, X0)`, the standard least-squares LRP; the plateau points then
also inform the descending line, which can shift the optimum.

```{r method}
fit_lrp(X, CV1, method = "ramp", step = 0.01)$parameters["Breakpoint"]
```

**`local_min_tol`** sets how close a competitor must fit to be flagged as
competing. It changes labelling only, never the fit: lower it (say `0.02`) to
mark only near-ties, raise it to mark any competitor.

```{r local-min-tol}
fit_lrp(X, CV2, local_min_tol = 0.02, step = 0.01)$local_minima$competing
```

## Warnings worth heeding

Two situations trigger a warning, and both mean the estimate deserves a second
look rather than direct use:

- **Non-negative slope**: the CV is not decreasing with plot size, so the
  "descent then plateau" reading does not hold for these data.
- **Breakpoint at the edge of the search range**: the trial may not span enough
  plot sizes to bracket the true optimum, so the CV may still be falling at the
  largest size evaluated.

A third situation is reported without a warning, because it is too common for
one to be informative: a **competing local minimum**, where another breakpoint
fits nearly as well and the estimate is not sharply identified. Check the
`competing` column of `$local_minima` before reporting a single value.

## Where to go next

`vignette("qrp")` and `vignette("mcm")` cover the other two CV-based methods,
which take the same inputs and usually give a larger and a smaller optimum
respectively. The QRP joins its plateau smoothly, so its SSE profile is a single
smooth basin and competing minima are rare. `vignette("replicates")` shows how
$CV_{Xo}$ feeds into the number of replications.
