--- title: "Modified Maximum Curvature (MCM)" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Modified Maximum Curvature (MCM)} %\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 two plateau models cut the CV curve into pieces. The Modified Maximum Curvature method takes a different route: it fits a single smooth curve and then asks where that curve *bends the most*. The relationship between CV and plot size is described by a power function, $$ CV_{(X)} = \frac{a}{X^{b}} + \varepsilon = a X^{-b} + \varepsilon $$ with $a > 0$ and $b > 0$, so the CV decreases as plots grow. This curve never becomes exactly flat, which is why the method cannot simply read off a plateau. Instead it uses the curvature of the line at each point, $$ K = \frac{y''}{\left(1 + y'^2\right)^{3/2}} $$ and takes the optimal plot size as the point of maximum curvature, the "elbow" where the steep initial drop turns into the slow tail. Solving $\mathrm{d}K/\mathrm{d}X = 0$ for $y = aX^{-b}$ gives a closed form: $$ X_o = \left[\frac{a^{2} b^{2} (2b + 1)}{b + 2}\right]^{\frac{1}{2b + 2}}, \qquad CV_{Xo} = a X_o^{-b} $$ Because it locates an elbow rather than a plateau, the MCM systematically returns a **smaller** optimum than the LRP and QRP, and correspondingly a higher CV at the optimum. This is expected, not an error. A note on the word "modified": in the original paper it refers to correcting $X_c$ by the cost factors $K_1$ and $K_2$ (labour per plot and per unit area). Modern Brazilian plot-size studies drop that cost step, and `fit_mcm()` follows them: $X_o$ is returned in the units of `x`. ### Two ways to estimate a and b This is the one place in the package where the two reference articles genuinely disagree, so both are available. **`method = "loglinear"`** is the classic route. Taking logs linearizes the model, $$ \log CV_{(X)} = \log a - b \log X + \varepsilon $$ so `a` and `b` come from a straight-line regression. Meier & Lessman weight that regression by the degrees of freedom of each point, following Federer (1955): CVs from large plots are computed from fewer plots and are less reliable, so they should count less. Pass `df` to enable that weighting. **`method = "nls"`** (the default) fits $CV = aX^{-b}$ directly on the original scale by nonlinear least squares, seeded from the log-log fit. This is what the modern articles do, and it is consistent with `fit_lrp()` and `fit_qrp()`. If `nls()` fails to converge, the function falls back to the log-linear estimate with a warning rather than erroring out. The two give visibly different answers, so choose deliberately: match whichever article you are comparing against. ### References Original method: Meier, V. D. & Lessman, K. J. (1971). Estimation of optimum field plot shape and size for testing yield in *Crambe abyssinica* Hochst. *Crop Science*, 11(5), 648-650. (After Lessman & Atkins, 1963.) Weighting: Federer, W. T. (1955). *Experimental Design*. Macmillan, New York. 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`) again. `n` is the number of plots of each size, returned by `calc_cv_shapes()`, and supplies the degrees of freedom for the optional weighting. ```{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 n <- tab1$n CV1 <- tab1$cv CV2 <- calc_cv_shapes(grid_mat("T2"))$cv CV3 <- calc_cv_shapes(grid_mat("T3"))$cv ``` Both `x` and `cv` must be strictly positive: the log-log fit and the starting values are undefined otherwise. ## Basic use ```{r fit} fit <- fit_mcm(x = X, cv = CV1) fit ``` $X_o \approx 4.5$ m² with $CV_{Xo} \approx 12.7\%$ for trial 1 -- as expected, well below the LRP's ~9.2 m². The MCM is the most conservative of the three methods. ```{r summary} summary(fit) ``` The closed form can be verified straight from the coefficients: ```{r check} a <- unname(fit$coefficients["a"]); b <- unname(fit$coefficients["b"]) c(formula = ((a^2 * b^2 * (2 * b + 1)) / (b + 2))^(1 / (2 * b + 2)), reported = unname(fit$parameters["Breakpoint"])) ``` ```{r predict} predict(fit, newx = c(1, 6, 18)) ``` Note there is no AIC or BIC here. The MCM is not nested with the plateau models and, under `"loglinear"`, is not even fitted on the same scale, so an information criterion would not be comparable. Compare methods on `R2` and `RMSE`, both reported on the original scale. ## Choosing the estimator ```{r methods} rbind( nls = fit_mcm(X, CV1, method = "nls")$parameters[c("Breakpoint", "R2")], loglinear = fit_mcm(X, CV1, method = "loglinear")$parameters[c("Breakpoint", "R2")], loglin_df = fit_mcm(X, CV1, method = "loglinear", df = n - 1)$parameters[c("Breakpoint", "R2")] ) ``` The differences are modest here but real, and they grow when the CVs of the large plots are noisy, which is exactly the case the weighting is meant to handle. The `df` weighting only applies with `method = "loglinear"`; `df = NULL` (the default) is the unweighted fit. ## Plot The figure has no flat segment: the power curve keeps descending, and the dotted guides mark the maximum-curvature point. ```{r plot} plot(fit, title = "Trial 1") ``` ```{r plot-ptbr} plot(fit, title = "Ensaio 1", decimal_mark = ",") ``` There is no `cond_word` argument, since the annotation has no conditional clause: a single power equation and $R^2$, plus $X_o$ and $CV_{Xo}$ by the breakpoint. ```{r save, eval = FALSE} plot(fit, title = "Trial 1", save = TRUE, file = "trial1_mcm.tiff", format = "tiff", dpi = 300) ``` ## Several trials at once ```{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_mcm(trials, x = "x", cv = "cv", trial = "trial") res ``` The three MCM breakpoints cluster below 5 m²; their mean is the overall recommendation, and all three sit well under the LRP and QRP optima. The degrees of freedom can be supplied as a column, applied to every trial: ```{r multi-df, eval = FALSE} trials$df <- rep(n - 1, 3) fit_mcm(trials, x = "x", cv = "cv", df = "df", trial = "trial", method = "loglinear") ``` ```{r multi-plot, eval = FALSE} plot(res, label_size = 3) ``` ## Warnings worth heeding - **`b <= 0`**: the CV is not decreasing with plot size, so a maximum-curvature point is meaningless for these data. - **nls fallback**: if the message says the fit fell back to the log-linear estimate, the reported values come from `"loglinear"`, not `"nls"`. The `method` field of the object records what was actually used. ## Where to go next `vignette("lrp")` and `vignette("qrp")` cover the plateau models, which use the same inputs and give larger optima. `vignette("paranaiba")` estimates the plot size straight from the raw grid instead of from CV values, and `vignette("replicates")` turns $CV_{Xo}$ into a number of replications.