Topics / Chapter 7

Comparing several treatments

To compare three or more fertilizer rates, varieties or practices in a single test and find out which ones differ.

Why not just run t tests over and over?

You want to compare four nitrogen rates. Pairwise t tests would need 6 comparisons. If each carries a 5% risk of error, the probability of wrongly declaring “a difference” in at least one of them is not 5% but roughly 26%:

1−0.956≈0.261 - 0.95^6 \approx 0.26

This risk grows quickly with the number of comparisons. Analysis of variance (ANOVA) compares all the groups in a single test and keeps the risk at 5%. That is why almost every field trial is evaluated with analysis of variance.

The idea: variability between and within groups

The total variability in a trial has two sources:

  • Treatment (between groups): different nitrogen rates giving different yields. This is the effect we want to find.
  • Error (within groups): even plots receiving the same rate yield differently. Soil, moisture, measurement differences, everything we can’t control; in other words, chance.

The F value is the ratio of the two:

F=Treatment mean squareError mean squareF = \frac{\text{Treatment mean square}}{\text{Error mean square}}

If nitrogen has no effect at all, differences between groups are also just chance and F comes out around 1. The larger F is than 1, the more the treatment effect dominates chance. The F table tells you how large it needs to be.

How should the trial be set up?

The model calculated here is for the completely randomized design. In this design treatments are assigned to plots completely at random and each treatment is repeated on several plots (replicates). The setup has two requirements:

  • Replication: each treatment should be tested on at least 2, preferably 3–4 plots. Error variability can’t be measured with a single plot.
  • Randomization: treatments must be assigned to plots at random. Otherwise the rate that lands in the fertile corner of the field gets an unfair advantage.

This design is suitable only if the field is uniform, that is, if there are no systematic differences between plots other than the treatments. If the field is not uniform (slope, direction of the irrigation line, an old field boundary), the randomized complete block design is used.

Question: Is there a difference among three or more treatments (fertilizer rates, varieties, irrigation intervals)? If so, which ones differ?
Example: Four nitrogen rates (kg/da) were tested with four replicates each in a completely randomized design. Yield (kg/da) was measured. Do the rates differ?

n = 4
n = 4
n = 4
n = 4
H0:μ1=μ2=⋯=μ4H_0: \mu_1 = \mu_2 = \cdots = \mu_{4}H1:At least one mean differs from the othersH_1: \text{At least one mean differs from the others}
Analysis of variance table
Source of variationdfSSMSF
Treatment323053.197684.465.67 **
Error121404.25117.02
Total1524457.44

df: degrees of freedom, SS: sum of squares, MS: mean square. * significant at 5%, ** at 1%, ns: not significant. CV = 2.97%.

Treatment means
TreatmentMeanChart
0 kg N308 c
5 kg N351.25 b
10 kg N396.25 a
15 kg N401.75 a

Duncan (5%) critical ranges: R2 = 16.67, R3 = 17.44, R4 = 17.92. Means sharing a letter are not significantly different.

Pairwise comparisons: Duncan’s multiple range test

Means are ranked; a larger critical difference is used for means far apart in rank and a smaller one for neighbours. Widely used in agricultural field trials.

  1. Rank the means from largest to smallest. The critical range depends on how many means (p) the pair spans.

    Rp=qαp; p; 12MSE4,αp=1−(1−α)p−1R_p = q_{\alpha_p;\,p;\,12} \sqrt{\dfrac{MS_E}{4}}, \qquad \alpha_p = 1 - (1 - \alpha)^{p-1}
  2. Critical ranges:

    R2=16.67,R3=17.44,R4=17.92R_{2} = 16.67,\quad R_{3} = 17.44,\quad R_{4} = 17.92
  3. If a wider range is not significant, the differences between means inside it are also declared not significant.

  4. Means whose difference is below the critical value share a letter. 5 of 6 pairs differ significantly.

ComparisonDifferencepCritical differenceResult
0 kg N – 5 kg N43.25216.67Significant
0 kg N – 10 kg N88.25317.44Significant
0 kg N – 15 kg N93.75417.92Significant
5 kg N – 10 kg N45216.67Significant
5 kg N – 15 kg N50.5317.44Significant
10 kg N – 15 kg N5.5216.67Not significant
012345673.49calculated 65.67 →F (3; 12)

The red area is the rejection region (α = 0.05). The green line is the value calculated from your data.

H₀ rejected: the difference is significant

The difference between treatments is statistically significant: at least one group mean differs from the others. See the comparison table for which groups differ.

p < 0.0001 (significant at 1%, **)

How was it calculated?

FcalcF_{calc} = 65.667

  1. Correction factor (CF): the square of the grand total divided by the number of observations.

    CF=(∑X)2N=5829216=2123577.56CF = \dfrac{(\sum X)^2}{N} = \dfrac{5829^2}{16} = 2123577.56
  2. Total sum of squares: the overall variability of all observations.

    SST=∑X2−CF=2148035−2123577.56=24457.44SS_T = \sum X^2 - CF = 2148035 - 2123577.56 = 24457.44
  3. Treatment sum of squares: variability coming from differences between groups.

    SSTr=12322+14052+15852+160724−CF=23053.19SS_{Tr} = \dfrac{1232^2 + 1405^2 + 1585^2 + 1607^2}{4} - CF = 23053.19
  4. Error sum of squares: variability within groups (due to chance).

    SSE=SST−SSTr=24457.44−23053.19=1404.25SS_E = SS_T - SS_{Tr} = 24457.44 - 23053.19 = 1404.25
  5. Divide the sums of squares by their degrees of freedom to get the mean squares.

    MSTr=SSTrk−1=23053.193=7684.396MS_{Tr} = \dfrac{SS_{Tr}}{k-1} = \dfrac{23053.19}{3} = 7684.396
  6. MSE=SSEN−k=1404.2512=117.021MS_E = \dfrac{SS_E}{N-k} = \dfrac{1404.25}{12} = 117.021
  7. F: the ratio of between-group to within-group variability. If the treatment has no effect it comes out around 1.

    F=MSTrMSE=7684.396117.021=65.667F = \dfrac{MS_{Tr}}{MS_E} = \dfrac{7684.396}{117.021} = 65.667
  8. Compare the calculated value with the table (critical) value.

    Fcalc=65.667>Ftab=3.49F_{calc} = 65.667 > F_{tab} = 3.49
  9. The calculated value is in the rejection region: H₀ is rejected. The p-value (< 0.0001) is smaller than α = 0.05.

  10. The coefficient of variation of the experiment shows how precisely it was carried out.

    CV=MSEXˉ⋅100=10.818364.313⋅100=2.97%CV = \dfrac{\sqrt{MS_E}}{\bar{X}} \cdot 100 = \dfrac{10.818}{364.313} \cdot 100 = 2.97\%

When the field isn’t uniform: the randomized complete block design

Real fields are rarely uniform. On a sloping field the plots at the foot are wetter and more productive than those at the top. Plots at the head of an irrigation line get more water than those at the end. Assigning treatments to such a field completely at random causes two problems:

  • By chance, most plots of one variety might land at the foot of the slope, and that variety looks better than it is.
  • The large difference between top and foot gets mixed into the error. As the error swells, F shrinks and real treatment differences disappear from view.

The randomized complete block design solves this:

  1. The field is divided into strips (blocks) perpendicular to the direction of variability. On a slope, each block is a strip at the same height, so the plots within a block are similar.
  2. Every treatment appears once in each block. The number of blocks equals the number of replicates.
  3. Treatments are randomized within each block.

The analysis now splits the total variability into three: block, treatment and error. The difference between top and foot goes into the block sum of squares and leaves the error. With a smaller error, treatment differences show up much more clearly:

SST=SSBl+SSTr+SSESS_T = SS_{Bl} + SS_{Tr} + SS_E

In the example below, 5 wheat varieties were tested in 4 blocks on a sloping field. Once blocks are separated, the difference between varieties is significant at the 1% level. The box under the calculator also shows what happens if the same data are analysed ignoring the blocks: the difference is not significant (p ≈ 0.30). Same data, two different decisions. The only difference is whether the design of the trial was taken into account in the analysis.

Question: When the field is not uniform, do the treatments differ? Once the block effect is removed, which treatments differ?
Example: Five wheat varieties were tested in 4 blocks on a sloping field. Blocks run across the slope: block 1 at the top, block 4 at the bottom. Yield (kg/da) was measured. Do the varieties differ?

TreatmentTotal
2051
2140
1992
2233
2044
Block total236225822677283910460

If a plot was lost (animals, flooding, a measurement error), leave its cell empty; its value will be estimated.

H0:μ1=μ2=⋯=μ5H_0: \mu_1 = \mu_2 = \cdots = \mu_{5}H1:At least one treatment mean differsH_1: \text{At least one treatment mean differs}
Analysis of variance table
Source of variationdfSSMSF
Block323823.67941.264.05 **
Treatment49042.52260.6318.23 **
Error121487.9123.99
Total1934354

df: degrees of freedom, SS: sum of squares, MS: mean square. * significant at 5%, ** at 1%, ns: not significant. CV = 2.13%.

Treatment means
TreatmentMeanChart
Variety 1512.75 c
Variety 2535 b
Variety 3498 c
Variety 4558.25 a
Variety 5511 c

Duncan (5%) critical ranges: R2 = 17.16, R3 = 17.96, R4 = 18.44, R5 = 18.76. Means sharing a letter are not significantly different.

Pairwise comparisons: Duncan’s multiple range test

Means are ranked; a larger critical difference is used for means far apart in rank and a smaller one for neighbours. Widely used in agricultural field trials.

  1. Rank the means from largest to smallest. The critical range depends on how many means (p) the pair spans.

    Rp=qαp; p; 12MSE4,αp=1−(1−α)p−1R_p = q_{\alpha_p;\,p;\,12} \sqrt{\dfrac{MS_E}{4}}, \qquad \alpha_p = 1 - (1 - \alpha)^{p-1}
  2. Critical ranges:

    R2=17.16,R3=17.96,R4=18.44,R5=18.76R_{2} = 17.16,\quad R_{3} = 17.96,\quad R_{4} = 18.44,\quad R_{5} = 18.76
  3. If a wider range is not significant, the differences between means inside it are also declared not significant.

  4. Means whose difference is below the critical value share a letter. 7 of 10 pairs differ significantly.

ComparisonDifferencepCritical differenceResult
Variety 1 – Variety 222.25217.16Significant
Variety 1 – Variety 314.75317.96Not significant
Variety 1 – Variety 445.5317.96Significant
Variety 1 – Variety 51.75217.16Not significant
Variety 2 – Variety 337418.44Significant
Variety 2 – Variety 423.25217.16Significant
Variety 2 – Variety 524317.96Significant
Variety 3 – Variety 460.25518.76Significant
Variety 3 – Variety 513217.16Not significant
Variety 4 – Variety 547.25418.44Significant

Did blocking help?

The blocks differ significantly (F = 64.05, p < 0.0001): the field really had a gradient, and blocking removed that variation from the error.

Relative efficiency: 10.95. A completely randomized design would need about 11 times as many replicates for the same precision.

Had the same data been analysed ignoring blocks (as a completely randomized design): F = 1.34, p = 0.3009, so the treatment difference would not be significant. Block variation would have inflated the error and hidden the real difference.

01234563.259calculated 18.23 →F (4; 12)

The red area is the rejection region (α = 0.05). The green line is the value calculated from your data.

H₀ rejected: the difference is significant

After removing the block effect, the difference between treatments is statistically significant. Letters show which treatments differ.

p < 0.0001 (significant at 1%, **)

How was it calculated?

FcalcF_{calc} = 18.232

  1. Correction factor (CF): the square of the grand total divided by the number of observations.

    CF=(∑X)2t⋅r=1046025⋅4=5470580CF = \dfrac{(\sum X)^2}{t \cdot r} = \dfrac{10460^2}{5 \cdot 4} = 5470580
  2. Total sum of squares: the overall variability of all plots.

    SST=∑X2−CF=34354SS_T = \sum X^2 - CF = 34354
  3. Block sum of squares: squares of the block totals divided by the number of treatments. Variability due to slope, soil or moisture differences in the field.

    SSBl=23622+25822+26772+⋯5−CF=23823.6SS_{Bl} = \dfrac{2362^2 + 2582^2 + 2677^2 + \cdots}{5} - CF = 23823.6
  4. Treatment sum of squares: squares of the treatment totals divided by the number of blocks.

    SSTr=20512+21402+19922+⋯4−CF=9042.5SS_{Tr} = \dfrac{2051^2 + 2140^2 + 1992^2 + \cdots}{4} - CF = 9042.5
  5. Error sum of squares: what remains after removing blocks and treatments, the truly random part.

    SSE=SST−SSBl−SSTrSS_E = SS_T - SS_{Bl} - SS_{Tr}
  6. SSE=34354−23823.6−9042.5=1487.9SS_E = 34354 - 23823.6 - 9042.5 = 1487.9
  7. Degrees of freedom: blocks r − 1, treatments t − 1, error (t − 1)(r − 1).

    dfB=3,dfM=4,dfH=4⋅3=12df_B = 3, \quad df_M = 4, \quad df_H = 4 \cdot 3 = 12
  8. Mean squares:

    MSTr=9042.54=2260.625MS_{Tr} = \dfrac{9042.5}{4} = 2260.625
  9. MSE=1487.912=123.992MS_E = \dfrac{1487.9}{12} = 123.992
  10. F value for treatments:

    F=MSTrMSE=2260.625123.992=18.232F = \dfrac{MS_{Tr}}{MS_E} = \dfrac{2260.625}{123.992} = 18.232
  11. Compare the calculated value with the table (critical) value.

    Fcalc=18.232>Ftab=3.259F_{calc} = 18.232 > F_{tab} = 3.259
  12. The calculated value is in the rejection region: H₀ is rejected. The p-value (< 0.0001) is smaller than α = 0.05.

  13. Coefficient of variation of the experiment:

    CV=MSEXˉ⋅100=11.135523⋅100=2.13%CV = \dfrac{\sqrt{MS_E}}{\bar{X}} \cdot 100 = \dfrac{11.135}{523} \cdot 100 = 2.13\%

Missing plot

A plot’s data can be lost to animal damage, flooding or a measurement error. In a randomized block design this is a problem, because the calculation assumes every treatment is present in every block. The solution is to estimate the missing value from the other plots of that treatment and that block (Yates’ formula):

X=t⋅T+r⋅B−G(t−1)(r−1)X = \frac{t \cdot T + r \cdot B - G}{(t-1)(r-1)}

Here TT is the sum of the known values for the missing plot’s treatment, BB the sum for its block, and GG the sum of all known plots. The estimate is “the value expected for this treatment and this block” and adds nothing to the error sum of squares. The analysis then proceeds as usual, with two corrections:

  • Because the estimate is not a real measurement, 1 is subtracted from the error degrees of freedom.
  • The estimate makes the treatment sum of squares slightly too large. This bias is calculated and subtracted from the treatment SS.

In the calculator, just leave the cell empty; the estimated value appears faintly in the cell. If several plots are missing, the formula is repeated for each in turn until the estimates stop changing. But the more plots are missing, the less reliable the trial becomes.

Relative efficiency

How much blocking gained is measured by the relative efficiency. In the example it is about 11: matching the same precision with a completely randomized design would take about 11 times as many replicates. A relative efficiency close to 1 means the field was uniform already and blocking gained nothing. Blocking has a cost too: the error degrees of freedom fall, so the table value rises slightly.

Don’t lay out blocks the wrong way: blocks must run perpendicular to the direction of variability. A block running along the slope contains everything from top to foot, so variability within the block doesn’t drop and blocking is useless. In this design every treatment must also be present in every block. If a plot is lost the analysis can still be done, but the missing value has to be estimated first, as explained above.

F is significant: what next?

The F test only says “at least one group differs”; it doesn’t say which. For that, a multiple comparison is made. The logic is the same in every method: if the difference between two means is larger than a critical difference, the two treatments differ. What separates the methods is how large that critical difference is chosen.

MethodCritical differenceWhen?
LSDSmallest. One value for all pairs: tα/22⋅MSE/rt_{\alpha/2} \sqrt{2 \cdot MS_E / r}Few treatments, comparisons planned in advance
DuncanGrows with how many means lie between the ranked pair (R2<R3<⋯R_2 < R_3 < \cdots)Very common in agricultural field trials
Tukey HSDLargest. One value for all pairs: qα; kMSE/rq_{\alpha;\,k} \sqrt{MS_E / r}Many treatments; when a false “difference” must especially be avoided

The number of pairwise comparisons rises fast with the number of treatments (10 pairs for 5 treatments, 45 for 10). Because LSD allows a 5% risk for every pair, the overall chance of at least one false “difference” grows. Tukey keeps that overall risk at 5% and therefore finds fewer differences. Duncan sits in between.

Results are shown with letters. Means sharing a letter are not significantly different. In the completely randomized example, the 10 and 15 kg N rates both get “a” whichever method you choose. Since 15 kg brings no extra benefit over 10 kg, the economical rate to recommend to farmers is 10 kg. In the randomized block example you can see the letters change when you switch method: Tukey groups together some varieties that LSD and Duncan call “different”.

Order matters: a multiple comparison is made only after the F test is significant. Moving on to pairwise comparisons when F is not significant brings back the error risk that analysis of variance was meant to avoid. Choose the method before looking at the results, too: switching methods until you like the letters makes the result misleading. Use whichever method your course or the journal expects.

The coefficient of variation of the trial

Field trial reports always quote a CV. It is the square root of the error mean square as a percentage of the grand mean, and it shows how carefully the trial was run:

CV=MSEXˉ×100CV = \frac{\sqrt{MS_E}}{\bar{X}} \times 100

In field trials a CV below 10–15% is generally considered good. A very high CV means there are large differences between plots other than the treatments, and treatment effects may be lost in that noise.

Assumptions

Analysis of variance assumes three things to give correct results: the data in each group are roughly normal, the group variances are similar, and the observations are independent. If one group’s values are much more scattered than the others, or the data are clearly skewed (as count and percentage data often are), the result can mislead. The data may then need to be transformed.