---
title: "Quadratic Response Plateau (QRP)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quadratic Response Plateau (QRP)}
  %\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

The Quadratic Response Plateau model answers the same question as the LRP, but
replaces the straight descent with a curve:

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

The plateau starts where the parabola reaches its vertex, so the two pieces meet
smoothly (no kink, unlike the LRP). That gives closed forms for both quantities
of interest:

$$
X_o = -\frac{b}{2c}, \qquad CV_{Xo} = p = a - \frac{b^2}{4c}
$$

For the curve to descend and then flatten, $b < 0$ and $c > 0$.

Because the descent is curved rather than straight, the QRP stays above a
straight line for longer and reaches its plateau later. In practice this means
**QRP almost always gives a larger optimum than LRP**, which in turn gives a
larger one than the MCM. That ordering (MCM < LRP < QRP) is reported across many
crops and is not a defect of any of the methods: they answer the same question
with different assumptions about the shape of the decay.

### How the breakpoint is estimated

The published implementations fit this model with `nls()` or `nlsLM()` and fixed
starting values. That is fragile: the starting values are not derived from the
data, and a failure to converge can pass silently.

`fit_qrp()` uses the same grid-search strategy as `fit_lrp()`. Once $X_o$ is
fixed, the model can be written as

$$
CV_{(X)} = A + c\,(X - X_o)^2 \ \text{ for } X \le X_o, \qquad
CV_{(X)} = A \ \text{ for } X > X_o
$$

which is **linear** in the plateau level $A$ and the curvature $c$. So each
candidate breakpoint is fitted by ordinary least squares, and the one with the
smallest residual sum of squares is returned. The reported `a`, `b`, `c` are
converted back to the usual parametrization via $a = A + cX_o^2$ and
$b = -2cX_o$. The result is numerically identical to a converged `nls()` fit,
without the starting values.

### References

Original method: Peixoto, A. P. B., Faria, G. A. & Morais, A. R. (2011).
Modelos de regressão com platô na estimativa do tamanho de parcelas em
experimento de conservação in vitro de maracujazeiro. *Ciência Rural*, 41(11),
1907-1913.

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 bundled simulated uniformity trial (`?uniformity_trial`), the same data
used throughout the package.

```{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
```

## 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_qrp(x = X, cv = CV1, step = 0.01)
fit
```

$X_o \approx 11.9$ m² for trial 1, noticeably larger than the LRP estimate of
about 9.2 m² on the same data: the smooth join always pushes the optimum out.

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

The QRP reports more fit statistics than the LRP: alongside `R2` and `RMSE` it
returns `R2_adj` (adjusted for the extra parameter), `MAE`, `SSE` and `MSE`.

The closed forms can be checked directly against the coefficients:

```{r check}
cf <- fit$coefficients
c(vertex = unname(-cf["b"] / (2 * cf["c"])),
  reported = unname(fit$parameters["Breakpoint"]))
```

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

## Plot

The figure follows the same layout as the LRP, with the quadratic term in the
annotated equation:

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

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

The curvature coefficient `c` is small, so it gets its own decimal control,
`digits_c` (4 by default), separate from `digits_coef` for `a` and `b`.

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

## Several trials at once

The data-frame interface is identical to `fit_lrp()`:

```{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_qrp(trials, x = "x", cv = "cv", trial = "trial", step = 0.01)
res
```

The summary table carries the extra `c` column, since the QRP has three
coefficients. The three breakpoints average 10.25 m², the article's QRP figure.

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

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

## Fine-tuning

`search_range` and `step` behave exactly as in `fit_lrp()`:

```{r tuning}
fit_qrp(X, CV1, search_range = c(6, 15), step = 0.01)$parameters["Breakpoint"]
fit_qrp(X, CV1, step = 0.01)$parameters["Breakpoint"]
```

There is no `method` argument here. The LRP has one because the published
procedure fits the descending line using only the pre-breakpoint points; the QRP
has no such variant.

## Warnings worth heeding

- **Non-positive curvature** (`c <= 0`): the parabola opens downward or is flat,
  so the "descend then plateau" shape 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 optimum.

## Comparing with the other methods

Fitting all three CV-based methods on the same trial shows the usual ordering:

```{r compare}
data.frame(
  method = c("MCM", "LRP", "QRP"),
  Xo = c(fit_mcm(X, CV1)$parameters["Breakpoint"],
         fit_lrp(X, CV1, step = 0.01)$parameters["Breakpoint"],
         fit_qrp(X, CV1, step = 0.01)$parameters["Breakpoint"]),
  row.names = NULL
)
```

Which one to report is a judgement call. The larger optimum is the conservative
choice: it buys more precision at the cost of more field area. The validation
article recommends the LRP value as its overall answer, while noting that LRP
and QRP delivered statistically indistinguishable precision at the optimum.

See `vignette("lrp")` and `vignette("mcm")` for the other two methods, and
`vignette("replicates")` for turning $CV_{Xo}$ into a number of replications.
