--- title: "Getting started with ofemeantest" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with ofemeantest} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = "center", fig.width = 6, fig.height = 4 ) ``` `ofemeantest` implements a cell-based, permutation-based protocol to compare treatments in **unreplicated on-farm experiments** (OFE) where strips of management are laid out side-by-side over a field and yield is recorded densely (e.g., from a yield monitor). The package handles three things: 1. **Build a regular grid** over the trial and aggregate the dense point data into cell-level medians, keeping only cells that fall inside a single treatment strip and contain enough observations. 2. **Estimate the effective sample size (ESS)** from the spatial autocorrelation of the residuals of a one-way ANOVA on cell medians, to avoid pseudo-replication from spatially correlated data. 3. **Test treatment differences** with a permutational ANOVA, repeated over many random subsamples of size ESS, and report the median *p*-value as the field-specific test result. ## Installation ```{r install, eval = FALSE} # install.packages("pak") pak::pkg_install("PPaccioretti/ofemeantest") ``` ## A worked example The bundled `ofe_f2` dataset comes from a corn (Zea mays L.) trial where a single fertilized strip (~2.2 ha) was compared against an adjacent control strip on the same field. Raw yield-monitor data were cleaned following Vega et al. (2019). ```{r load} library(ofemeantest) data("ofe_f2") head(ofe_f2) ``` ### Quick run Call `ofemt()` and let it build the grid internally: ```{r quick, message = FALSE, warning = FALSE} res <- ofemt( data = ofe_f2, y = "Yield_tn", x = "Treatment", cellsize = 9, min_per_cell = 4, n_p = 2000, n_s = 200, alpha = 0.05 ) res ``` Key fields in the printed result: - **Selected cells**: cells retained after filtering by single treatment and `min_per_cell`. - **n / ESS**: nominal sample size and effective sample size given the spatial autocorrelation of the residuals. - **Rho / Moran's I**: spatial autocorrelation indicators. - **Means comparison**: per-treatment median yield (in cells) and compact letter display from the multiple-comparison procedure. - **ANOVA permutation test**: median *p*-value per pairwise comparison. ### Step-by-step For finer control — inspecting the grid, tweaking the cell size, or reusing the same selection in several downstream analyses — build the grid explicitly with `make_ofe_grid()` and pass it to `ofemt()`. ```{r stepwise, message = FALSE, warning = FALSE} g <- make_ofe_grid( data = ofe_f2, x = "Treatment", cellsize = 9, min_per_cell = 4 ) names(g) plot_grid_selection(g, data = ofe_f2) res2 <- ofemt( data = ofe_f2, y = "Yield_tn", x = "Treatment", grid = g, n_p = 2000, n_s = 200 ) ``` Both calls produce the same numeric output when the grid is built with matching arguments — `ofemt(grid = NULL, ...)` is equivalent to `ofemt(grid = make_ofe_grid(...))` under the hood. ### Tuning the grid `plot_grid_selection()` draws the three layers on one set of axes: the full grid in grey, the cells that survived the filters shaded in blue, and the observations as points. Overlaying the points on the cell boundaries is what makes the geometric arguments legible — a cell is dropped either because it straddles two treatments or because too few points landed inside it, and both are visible at a glance. By default the plot is built with **ggplot2**, which keeps the legend outside the panel — it can never land on top of the data, and the result does not change with the size of the graphics device. Pass `engine = "base"` for base graphics instead; if ggplot2 is not installed, that is what you get anyway, with a message. Because the ggplot2 engine returns a `ggplot` object, you can keep customising it: ```{r engine, eval = FALSE} plot_grid_selection(g, data = ofe_f2) + ggplot2::labs(subtitle = "Lote 2, campaña 21/22") plot_grid_selection(g, data = ofe_f2, engine = "base") ``` With dense yield-monitor data the observations can swamp the cell boundaries; lower `point_size` until the grid shows through. The title of each plot repeats the parameters used, so successive calls can be compared directly: ```{r tuning, eval = FALSE} op <- par(mfrow = c(1, 2)) plot_grid_selection( make_ofe_grid(ofe_f2, x = "Treatment", cellsize = 9, min_per_cell = 4), data = ofe_f2 ) # Shift the origin by half a cell and rotate to follow the strips plot_grid_selection( make_ofe_grid( ofe_f2, x = "Treatment", cellsize = 9, min_per_cell = 4, shift = c(4.5, 4.5), angle_deg = 10, buffer = 5 ), data = ofe_f2 ) par(op) ``` ### Keeping the geometries in the result By default `ofemt()` returns tables only. Set `keep_components` to embed the spatial objects in the result, which lets you plot the analysis that actually ran rather than rebuilding the grid by hand: | `keep_components` | Adds to the result | `plot()` shows | |---|---|---| | `"none"` (default) | nothing | *(errors — no geometries)* | | `"light"` | `grid`, the full `ofe_grid` (`grid_all`, `grid_sel`, `cell_stats`, `params`) | grid + selected cells | | `"full"` | `"light"` plus `cell_medians` (one point per selected cell, with its median response and ANOVA residual) and `points_joined` (every observation with its `CellID`) | grid + selected cells + points | Cost scales accordingly: `"light"` stores one polygon per grid cell, `"full"` adds one row per observation. ```{r keep, eval = FALSE} res_full <- ofemt( ofe_f2, y = "Yield_tn", x = "Treatment", cellsize = 9, min_per_cell = 4, keep_components = "full" ) # No further arguments needed — every layer comes from the object itself plot(res_full) # The per-cell medians and residuals that fed the spatial diagnostics head(res_full$cell_medians) ``` ### Reproducibility `ofemt()` exposes a `seed` argument that controls the random subsampling inside the permutation runs. The default is `seed = 7L`, so two calls with the same inputs return identical p-values: ```{r seed, eval = FALSE} identical( ofemt( ofe_f2, y = "Yield_tn", x = "Treatment", cellsize = 9, min_per_cell = 4, seed = 7L ), ofemt( ofe_f2, y = "Yield_tn", x = "Treatment", cellsize = 9, min_per_cell = 4, seed = 7L ) ) #> TRUE ``` Pass `seed = NULL` to let results vary across runs (e.g., when exploring sensitivity to the random draws). ### Inspecting the permutation distribution Each call retains the per-run *p*-values in `res$perm_runs` so you can sanity-check that the median *p*-value is not an artefact of a long tail. `perm_runs` has one row per comparison per sampling run, with columns `Comparison`, `p_value` (that run's permutation *p*-value), `p_adj` (the same value after adjusting for multiplicity *within* the run) and `run`. `plot_pvalue_hist()` draws that distribution with two reference lines per panel: a solid line at the median — the value reported in the `ANOVA permutation test` table — and a dashed line at `alpha`. ```{r hist, eval = FALSE} # Requires the optional 'ggplot2' package. plot_pvalue_hist(res) ``` When the analysis was run with a multiplicity adjustment, the histogram shows the *adjusted* values by default, so the median line and the reported `p_adj` always refer to the same quantity. Use `which = "raw"` to see the unadjusted distribution instead: ```{r hist_adj, eval = FALSE} res_bonf <- ofemt( ofe_f2, y = "Yield_tn", x = "Treatment", cellsize = 9, min_per_cell = 4, p_adjust_method = "bonferroni" ) plot_pvalue_hist(res_bonf) # adjusted plot_pvalue_hist(res_bonf, which = "raw") # unadjusted ``` ## Where to go next - `vignette("methodology")` — technical description of the protocol (effective sample size, permutational ANOVA, multiplicity adjustment). - `?ofemt` — full argument reference for the main function. - `?make_ofe_grid` — details on grid construction and selection. ## References - Córdoba, M., Paccioretti, P., & Balzarini, M. (2025). A new method to compare treatments in unreplicated on-farm experimentation. Precision Agriculture, 26(1), 4. https://doi.org/10.1007/s11119-024-10206-0 - Griffith, D.A. (2005). Effective geographic sample size in the presence of spatial autocorrelation. *Annals of the Association of American Geographers* 95(4): 740–760. - Vega A., Córdoba M., Balzarini M. (2019). Protocol for automating error removal from yield maps. *Precision Agriculture* 20: 1030–1044.