Guide
How to Analyse an RCBD Field Trial in R: ANOVA, Assumption Checks, LSD and DMRT with agricolae
In short
To analyse an RCBD trial in R, fit aov(response ~ block + treatment) with both columns as factors, check residual normality (Shapiro-Wilk) and equal variances (Bartlett or Levene) plus residual plots, then separate means with agricolae's LSD.test() or duncan.test(). Their $groups output gives the compact letters you plot on a ggplot2 bar chart with standard-error bars.
- Design
- Randomised complete block (one plot per treatment per block)
- Model
- aov(yield ~ block + treatment)
- Error df
- (t - 1)(r - 1)
- Assumption checks
- Shapiro-Wilk on residuals, Bartlett, Levene, residual plots
- Mean separation
- agricolae LSD.test() and duncan.test()
- Code tested in
- R 4.5.2
What is an RCBD, and why does the block term matter?
A randomised complete block design (RCBD) is the workhorse of agronomic field trials. The field is divided into blocks that are as uniform as possible inside, and every treatment appears once in every block, in a random order. Blocks usually follow a known gradient: slope, distance from an irrigation channel, a previous crop, or shade from a tree line. Because each block holds a full set of treatments, differences between blocks can be removed from the error term instead of inflating it.
That is the whole point of the analysis below. If you analyse an RCBD as if it were a completely randomised design, the block-to-block variation stays in the residual, the error mean square grows, and real treatment differences can disappear. The statistical model for one observation is:
y_ij = mu + tau_i + beta_j + e_ij
where mu is the overall mean, tau_i the effect of treatment i, beta_j the effect of block j, and e_ij the random error. There is one observation per treatment per block, so there is no treatment-by-block interaction to estimate; that variation is the error.
Step 1: lay the data out correctly
Keep one row per plot and one column per variable. Block and treatment must be factors, not numbers, or R will fit them as continuous slopes. A hormone concentration such as 100 ppm is still a treatment level in an RCBD analysis, so store it as a factor too.
| Column | Type in R | Example values |
|---|---|---|
| block | factor | B1, B2, B3, B4 |
| treatment | factor | Control, GA3_50, GA3_100, NAA_50, NAA_100 |
| yield (or any response) | numeric | 4.6 (t ha-1) |
The example below uses a small, made-up dataset (five treatments, four blocks) so that the code runs as written. The numbers are illustrative only and are not results from any published trial. With your own data, replace the data.frame() call with read.csv("your_file.csv", stringsAsFactors = TRUE).
install.packages(c("agricolae", "car", "ggplot2")) # once
library(agricolae)
library(car)
library(ggplot2)
trial <- data.frame(
block = factor(rep(c("B1", "B2", "B3", "B4"), each = 5)),
treatment = factor(rep(c("Control", "GA3_50", "GA3_100", "NAA_50", "NAA_100"), times = 4)),
yield = c(3.4, 4.1, 5.2, 4.0, 4.6,
3.1, 4.4, 5.6, 3.8, 4.9,
3.6, 4.0, 5.0, 4.3, 4.4,
3.3, 4.6, 5.4, 4.1, 4.8)
)
str(trial)
Step 2: fit the two-way ANOVA with aov()
The model formula lists the block and the treatment as additive terms. Putting block first is a convention; with balanced data and one plot per treatment per block, the order does not change the sums of squares.
model <- aov(yield ~ block + treatment, data = trial)
summary(model)
cv.model(model) # coefficient of variation (%), from agricolae
With t treatments and r blocks, the degrees of freedom in the ANOVA table should always add up as follows. Checking them is the fastest way to catch a block or treatment column that was read as numeric.
| Source | Degrees of freedom | In the example (t = 5, r = 4) |
|---|---|---|
| Block | r - 1 | 3 |
| Treatment | t - 1 | 4 |
| Error | (t - 1)(r - 1) | 12 |
| Total | tr - 1 | 19 |
Read the treatment line first: the F value is the treatment mean square divided by the error mean square, and its p-value tells you whether at least one treatment mean differs. The block line tells you whether blocking was worth it; a large block F means the blocks absorbed real field variation. Report the coefficient of variation alongside the ANOVA, because reviewers use it to judge how well the trial was managed.
Step 3: check the assumptions before you trust the F-test
ANOVA assumes the residuals are roughly normal, that variances are similar across treatments, and that errors are independent (randomisation takes care of the last one). Test and plot; do not rely on one p-value.
shapiro.test(residuals(model)) # normality of residuals
bartlett.test(yield ~ treatment, data = trial) # equal variances (sensitive to non-normality)
leveneTest(yield ~ treatment, data = trial) # equal variances (more robust), from car
par(mfrow = c(1, 2))
plot(model, which = 1:2) # residuals vs fitted, normal Q-Q
- Shapiro-Wilk on the residuals, not on the raw yield column. Raw yields are a mixture of treatment groups and are often non-normal even when the model is fine.
- Bartlett's test is powerful but reacts to non-normality as well as unequal variance. Levene's test (median-centred by default in
car) is the safer second opinion. - Residuals vs fitted should look like a shapeless band. A funnel that widens with the mean is the classic sign that a transformation is needed.
- Normal Q-Q points should sit close to the line, with only mild departures at the tails.
If the checks fail, the usual agronomic fixes are a square-root transformation for counts (insects per plant, pods per plant), an arcsine square-root or logit transformation for percentages, and a log transformation when the spread grows with the mean. Re-fit the model on the transformed variable, re-check, and report means on the original scale.
Step 4: separate the means with LSD and DMRT
A significant F-test does not say which treatments differ. The agricolae package, maintained by Felipe de Mendiburu, provides the post-hoc tests most often requested by agricultural journals in South Asia. Both functions accept the fitted model directly and pull the error mean square and error degrees of freedom from it.
lsd <- LSD.test(model, "treatment", p.adj = "none", console = TRUE)
dmrt <- duncan.test(model, "treatment", console = TRUE)
lsd$groups # means with compact letters
dmrt$groups
| Test | Function | What it controls | When to use it |
|---|---|---|---|
| Fisher's LSD | LSD.test() | Per-comparison error only (with p.adj = "none") | After a significant F-test, with a small number of planned treatments |
| Duncan's multiple range test (DMRT) | duncan.test() | Critical range grows with the number of means spanned | Common in Nepali and Indian agronomy journals; widely expected by reviewers |
| Tukey's HSD | HSD.test() | Family-wise error across all pairs | Many treatments, or when a reviewer asks for a conservative test |
The $groups table is the compact-letter display: treatments that share a letter are not significantly different at the chosen alpha (0.05 by default). In the example data, GA3_100 gets "a", NAA_100 "b", GA3_50 and NAA_50 share "c", and the control is "d". Apply LSD or DMRT only after a significant F-test; if the F-test is not significant, report that and stop, rather than hunting for letters.
Step 5: plot the means with letters in ggplot2
A bar chart with standard-error bars and letters is the figure most agronomy papers expect. The $means table from LSD.test() already contains a standard error for each mean, calculated from the pooled error mean square, so the bars match the test.
plot_df <- data.frame(treatment = rownames(lsd$means),
mean = lsd$means$yield,
se = lsd$means$se)
letters_df <- data.frame(treatment = rownames(lsd$groups),
groups = lsd$groups$groups)
plot_df <- merge(plot_df, letters_df, by = "treatment")
p <- ggplot(plot_df, aes(x = reorder(treatment, -mean), y = mean)) +
geom_col(fill = "grey70", width = 0.6) +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.15) +
geom_text(aes(y = mean + se + 0.2, label = groups)) +
labs(x = "Treatment", y = expression(Yield~(t~ha^-1))) +
theme_classic()
p
ggsave("yield_letters.png", p, width = 6, height = 4, dpi = 300)
The column called yield in lsd$means takes the name of your response variable, so change lsd$means$yield if your column is pod_yield or plant_height. The whole script above was run end to end in R 4.5.2 before publication.
Common mistakes in RCBD analysis
- Dropping the block term. It moves block variation into the error and weakens every test.
- Numeric factors. Blocks coded 1-4 and doses coded 50, 100, 150 must be converted with
factor(), or the degrees of freedom will be wrong. - Averaging subsamples first, or not at all. If you measured five plants per plot, average them to one plot value before the ANOVA. Treating plants as replicates is pseudo-replication.
- Letters without a significant F. Mean separation is only meaningful after the overall test.
- Testing normality on raw data instead of on model residuals.
How this workflow connects to real field trials
This is the same analysis structure that sits behind the RCBD field experiments in Bishal Neupane's publication record. The okra growth-regulator trial in Sunsari (first author Manwee Gupta, with Neupane as second and corresponding author) compared GA3 and NAA foliar sprays in randomised blocks, and the rice genotype evaluation in eastern Terai (first author Laxmi Katuwal) tested nine spring rice genotypes against Hardinath-1 in the same design. During his PMAMP Morang internship in 2025, Neupane ran his own RCBD okra experiment across three hormone concentrations and analysed more than 200 morphological and yield observations in R with ANOVA and LSD, the steps shown above. The full list of papers is on the publications page, and the wider toolkit, including Google Earth Engine and Python, is on the methods page.
For a satellite-side companion to this field-statistics guide, see vegetation indices from Sentinel-2.
What to report in your results section
- The design: RCBD, number of treatments, number of replications (blocks), and plot size.
- The software: R version, and the agricolae and car packages.
- The assumption checks you ran and any transformation applied.
- The F-test result for treatments, the coefficient of variation, and the mean separation test with its alpha level.
- Means with letters in a table or figure, on the original scale.
Supervisors and research groups interested in field experimentation combined with geospatial analysis are welcome to get in touch through the contact page.
Frequently asked questions
How do I run ANOVA for an RCBD in R?
Store block and treatment as factors, then fit model <- aov(yield ~ block + treatment, data = trial) and run summary(model). The block term removes field variation from the error. Check that the error degrees of freedom equal (t - 1)(r - 1); if they do not, a factor was probably read as a number. Add cv.model(model) from agricolae to report the coefficient of variation.
Should I use LSD or Duncan's test for my agriculture thesis?
Both are acceptable after a significant F-test. Fisher's LSD is simple and suits a few planned treatments; Duncan's multiple range test adjusts the critical range by how many means are compared and is widely used in Nepali and Indian agronomy journals. If you have many treatments or a reviewer wants stricter control, use Tukey's HSD with HSD.test().
Do I test normality on raw yield data or residuals?
Test the residuals of the fitted model with shapiro.test(residuals(model)). Raw yields combine several treatment groups with different means, so they can look non-normal even when the ANOVA assumptions hold. Always pair the test with a normal Q-Q plot, because with small trials the test has little power and with large ones it flags trivial departures.
How do I add significance letters to a ggplot2 bar chart?
Run LSD.test() or duncan.test() from agricolae, take the treatment names and letters from the $groups table, and merge them with the means and standard errors from the $means table. Then add geom_text(aes(label = groups)) positioned just above each error bar. Treatments sharing a letter are not significantly different at alpha 0.05.
What if my RCBD data fail the assumption checks?
Transform and re-fit. Square-root transformation suits counts such as insects or pods per plant, arcsine square-root or logit suits percentages, and log suits data whose spread grows with the mean. Re-check the residuals on the transformed scale and present back-transformed means. If no transformation works, consider a non-parametric Friedman test, which also respects the block structure.