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").
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.
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:
min_per_cell observations,
andMixed-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.
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.
Given the ESS, ofemt() runs n_s independent
permutational ANOVAs. Each run:
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().
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.