Skip to contents

Methodology

Source: vignettes/methodology.Rmd
methodology.Rmd

This vignette describes the statistical protocol implemented by the ofemeantest package, intended for readers who want to understand what the function computes — beyond the user-facing recipe in vignette("getting-started").

1. Problem setting

On-farm experiments (OFE) on commercial fields are typically laid out as side-by-side strips with no replication of treatments and no randomisation. Yield-monitor data over these strips are dense, irregular and strongly spatially autocorrelated, so classical ANOVA on the raw observations grossly inflates the effective sample size and breaks its independence assumption.

The protocol implemented here adapts ANOVA to this setting by (i) aggregating the raw observations into a regular grid, (ii) estimating the autocorrelation-corrected effective sample size, and (iii) inferring treatment effects from a distribution of permutational p-values computed on random subsamples of that ESS.

2. Grid construction

make_ofe_grid() builds a regular grid (square cells, optional rotation and origin shift) over the bounding box of the trial data and assigns each yield-monitor point to a cell. A cell is selected for analysis if and only if:

  • it contains at least min_per_cell observations, and
  • all those observations come from a single treatment.

Mixed-treatment cells (those that straddle the boundary between strips) are dropped. The per-cell response used downstream is the median of the raw observations falling in that cell, which is robust to the typical outliers of yield-monitor data.

Implementation details: the grid origin and cell size are passed through deterministically; the same arguments always produce the same cells, so ofemt(grid = NULL, ...) and ofemt(grid = make_ofe_grid(...)) are interchangeable.

3. Effective sample size (ESS)

A one-way ANOVA on cell medians yields residuals \(e_i\). Their spatial autocorrelation is estimated with the approximate profile likelihood estimator (APLE) of \(\rho\) (Li, Calder & Cressie, 2007), and the effective sample size is computed from Griffith (2005):

\[ n^{*} = n \left(1 - a \cdot \frac{n-1}{n} \cdot \left(1 - \exp(-b\rho + c\sqrt{\rho})\right)\right) \]

with the coefficients given in the published formula. We also report Moran’s I of the residuals as a sanity check on the magnitude of the autocorrelation. \(\rho\) is constrained to \([0, 1]\) before entering the ESS formula.

4. Permutational ANOVA on subsamples

Given the ESS, ofemt() runs n_s independent permutational ANOVAs. Each run:

  1. samples \(\lceil n^{*}/k \rceil\) cells without replacement from each of the \(k\) treatments;
  2. fits a one-way ANOVA on the subsample;
  3. permutes the treatment labels n_p times and records the proportion of permuted F-statistics that exceed the observed one (this is the per-run p-value).

The reported test statistic for the comparison is the median of the n_s per-run p-values. The full distribution is retained in res$perm_runs and can be inspected with plot_pvalue_hist().

5. Multiple comparisons

When the trial contains more than two treatments, every pairwise comparison is tested as above and the results are adjusted for multiplicity via stats::p.adjust() (Bonferroni, Holm or Benjamini–Hochberg, selected by the p_adjust_method argument).

The adjustment is applied within each run, across the \(\binom{k}{2}\) comparisons produced by that run, and the reported p_adj is the median of the adjusted values. Adjusting inside the run keeps the two summaries on the same footing: p_value and p_adj are both medians of a genuine empirical distribution over the n_s runs, and the median line drawn by plot_pvalue_hist() is exactly the value reported in the table. (Adjusting the median instead — the reverse order — would mix a summary statistic with a correction computed from a single vector of medians, and the histogram could then disagree with the table.) Both per-run quantities are kept in res$perm_runs as p_value and p_adj.

A compact letter display is then computed from the adjusted p-values with multcompView::multcompLetters(). Treatments are ordered by decreasing median response before the letters are assigned, so "a" marks the highest-yielding group and the letters read monotonically down the Means comparison table. Treatment labels are swapped for internal placeholders before the letters are computed and restored afterwards, so labels containing spaces, +, -, parentheses or accents are reported verbatim.

6. Assumptions and limitations

  • The protocol assumes the trial layout produces enough cells per treatment for the subsampling step to be meaningful.
  • The ESS formula is an approximation valid for approximately normally distributed data with non-negative autocorrelation.
  • The current implementation supports a single categorical treatment. Multi-factor designs, factor interactions and continuous covariates are not yet supported.

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.
  • Li, H., Calder, C.A. & Cressie, N. (2007). Beyond Moran’s I: Testing for spatial dependence based on the spatial autoregressive model. Geographical Analysis 39(4): 357–375.
  • Vega A., Córdoba M., Balzarini M. (2019). Protocol for automating error removal from yield maps. Precision Agriculture 20: 1030–1044.