Pooled analysis of multi-environment trials with CANE

The problem

A single experiment at one site in one season answers a narrow question. Crop improvement and agronomy therefore repeat the same experiment across locations, seasons or years, and the quantity of interest becomes the treatment-by-environment interaction: do treatments rank consistently, or does the best treatment depend on where and when it is grown?

Pooling such a series is not simply a matter of fitting one larger model. The classical workflow has four prerequisites, and each is easy to get wrong:

each environment must first be analysed on its own, with every effect tested against the error term that the design assigns to it;
the error variances of the environments must be compared, because heterogeneous errors invalidate the pooled F-tests;
if they are heterogeneous, the data must be rescaled before pooling; and
the pooled treatment means must be separated using the correct error stratum, and interpreted in the light of the interaction.

CANE performs all four in one design-specific function call.

The five functions

data.frame(
  Design = c("Completely randomised", "Randomised block", "Latin square",
             "Split plot", "Strip plot"),
  Function = c("PooledCRD()", "PooledRBD()", "PooledLSD()", "PooledSPD()",
               "PooledStPD()"),
  `Error strata` = c(1, 1, 1, 2, 3),
  `Design arguments` = c("-", "Replication", "Row, Column",
                         "Replication, MainPlot, SubPlot",
                         "Replication, RowFactor, ColumnFactor"),
  check.names = FALSE)
#>                  Design     Function Error strata                     Design arguments
#> 1 Completely randomised  PooledCRD()            1                                    -
#> 2      Randomised block  PooledRBD()            1                          Replication
#> 3          Latin square  PooledLSD()            1                          Row, Column
#> 4            Split plot  PooledSPD()            2       Replication, MainPlot, SubPlot
#> 5            Strip plot PooledStPD()            3 Replication, RowFactor, ColumnFactor

All five share the same interface: data, Response, Location, the design-specific factors, alpha, Mult_Comp_Test, and optionally transform and env_effect. There is no omnibus function that handles every design; the user chooses the function matching the design that was actually conducted.

A worked randomised block series

The sugarcane dataset holds commercial cane sugar yield (t/ha) of 18 genotypes grown in a randomised block design with two replications at two locations over three crop cycles, giving six environments.

data(sugarcane)
str(sugarcane)
#> 'data.frame':    216 obs. of  6 variables:
#>  $ env : Factor w/ 6 levels "SGA-IP","SGA-IIP",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ loc : Factor w/ 2 levels "SGA","TCSM": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ crop: Factor w/ 3 levels "IP","IIP","R": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ rep : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ gen : Factor w/ 18 levels "Co 12008","Co 12009",..: 1 2 3 4 5 6 7 8 9 10 ...
#>  $ ccsy: num  16.19 13.55 19.12 6.27 20.83 ...
with(sugarcane, table(env, rep))
#>           rep
#> env         1  2
#>   SGA-IP   18 18
#>   SGA-IIP  18 18
#>   SGA-R    18 18
#>   TCSM-IP  18 18
#>   TCSM-IIP 18 18
#>   TCSM-R   18 18

fit <- PooledRBD(sugarcane, Response = "ccsy", Location = "env",
                 Treatment = "gen", Replication = "rep",
                 alpha = 0.05, Mult_Comp_Test = 1)

Step 1: the environments on their own

fit$Individual_ANOVA[["SGA-IP"]]
#>  Source      Df Sum Sq Mean Sq F value Pr(>F)  Sig Error
#>  Replication  1 89.164 89.16   8.314   0.01032 *   Error
#>  Treatment   17 682.24 40.13   3.742   0.00475 **  Error
#>  Error       17 182.32 10.72                            
#>  Total       35 953.72                                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 'Error' names the mean square used as the denominator of each F-test.

The error mean squares of the six environments are what the next step compares:

fit$Environment_Key
#>      Level  n   Error_MS Error_df
#> 1   SGA-IP 36 10.7245624       17
#> 2  SGA-IIP 36 20.5435136       17
#> 3    SGA-R 36 10.6075908       17
#> 4  TCSM-IP 36  0.1605052       17
#> 5 TCSM-IIP 36  1.3576430       17
#> 6   TCSM-R 36  3.3007993       17

Step 2: are the error variances homogeneous?

fit$Bartlett_Test
#> 
#>  Bartlett's test for homogeneity of error variances (on error mean squares)
#> 
#> data:  environment-wise error mean squares
#> Bartlett's K-squared = 81.17, df = 5, p-value = 4.77e-16
#> ratio of largest to smallest error mean square = 128
fit$Levene_Test
#> 
#>  Levene's test (median-centred, Brown-Forsythe) on environment-wise residuals,
#>  denominator df = total error df
#> 
#> data:  environment-wise residuals
#> F = 12.152, num.df = 5, denom.df = 102, p-value = 3.033e-09

The largest error mean square is more than a hundred times the smallest, so the error variances are clearly heterogeneous. Bartlett’s test is sensitive to non-normality, which is why CANE reports the median-centred Brown-Forsythe form of Levene’s test beside it, and the environment-wise Shapiro-Wilk tests as well:

fit$Normality_Test
#>   Environment  n         W     p.value
#> 1      SGA-IP 36 0.9664854 0.337459671
#> 2     SGA-IIP 36 0.9618970 0.245817020
#> 3       SGA-R 36 0.9734697 0.527607435
#> 4     TCSM-IP 36 0.9810379 0.779801819
#> 5    TCSM-IIP 36 0.8850167 0.001362167
#> 6      TCSM-R 36 0.9939775 0.999237292

Both tests agree here, so the decision to rescale is secure. When they disagree, prefer Levene, and use transform = "always" or transform = "never" to take the decision yourself rather than leaving it to a threshold.

Step 3: the transformation

fit$Transformation$applied
#> [1] TRUE
cat(fit$Transformation$reason)
#> transform = "auto": Bartlett's test gave P = 4.77e-16 against alpha = 0.05, so the error variances were judged heterogeneous and Aitken's transformation was applied. Levene's test (robust to non-normality) gave P = 3.033e-09.

Aitken’s transformation divides each observation by the square root of the error mean square of its own environment. That yields a useful self-check: after the transformation each environment contributes an error sum of squares equal to its own error degrees of freedom, so the pooled error mean square must come out at exactly 1.

cat(fit$Transformation$diagnostic$note)
#> On the Aitken scale the pooled error mean square of a single-error design equals 1 by construction, because each environment's error sum of squares is divided by its own error mean square. Observed: 1.

Step 4: the pooled analysis

fit$Pooled_ANOVA
#>  Source          Df  Sum Sq Mean Sq F value Pr(>F)  Sig Error         
#>  Environment       5 30753. 6151.   1654.   2.5e-09 *** Rep within Env
#>  Rep within Env    6 22.312 3.719   3.719   0.0022  **  Pooled error  
#>  Treatment        17 503.72 29.63   29.63   < 1e-16 *** Pooled error  
#>  Treatment x Env  85 1233.5 14.51   14.51   < 1e-16 *** Pooled error  
#>  Pooled error    102 102.00 1.000                                     
#>  Total           215 32614.                                           
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 'Error' names the mean square used as the denominator of each F-test.

The Error column names the mean square used as the denominator of each F-test. In a randomised block series the environment main effect is tested against the replication-within-environment mean square, because replications are the units within which environments are compared; the treatment effect and the interaction are tested against the pooled error.

Note the caveat attached to the environment main effect on a transformed scale: each environment has been divided by its own error standard deviation, so that term is no longer expressed in t/ha and is not interpretable as a yield difference. Interpretation belongs to the treatment and interaction terms.

Step 5: separating the treatment means

fit$Treatments_Comparison
#> Multiple comparison of means
#>   Procedure     : Tukey's HSD
#>   Error stratum : Pooled error  (MSerror = 1 on 102 df)
#>   Scale         : Aitken-transformed
#>   alpha = 0.05 | treatments = 18 | replication = 12 | critical difference = 1.462
#> 
#>  Treatment    Mean Groups Mean_original_scale
#>   Co 12009 13.9320      a             15.8937
#>   Co 16018 13.6682      a             15.5943
#>   Co 17012 12.8414     ab             15.1987
#>   Co 14027 12.6390     ab             15.2795
#>   Co 12008 11.8366     bc             12.8636
#>   Co 86032 11.8043    bcd             14.2367
#>   Co 16010 11.6770    bcd             14.3450
#>   Co 17001 11.3938   bcde             15.8328
#>   Co 17003 11.3920   bcde             13.5122
#>   Co 16009 11.0494   cdef             12.2808
#>   Co 18024 10.5471   cdef             12.8915
#>   Co 18009 10.3948  cdefg             13.0860
#>   Co 14005 10.3435   defg             12.9514
#>   Co 14002 10.0568    efg             14.7262
#>   Co 17004 10.0188    efg             13.5124
#>   Co 15018  9.7537     fg             11.1913
#>   Co 15015  9.0423     gh             10.7336
#>   Co 17013  7.8221      h              8.0148
#> 
#> Means sharing a letter do not differ significantly at alpha = 0.05.

Because the analysis was carried out on the transformed scale, the means in the Mean column and the grouping letters both refer to that scale. The arithmetic means on the original scale are given alongside, for reporting in t/ha. They are not an inverse transformation: no single multiplier inverts Aitken’s rescaling, because each environment was divided by a different constant.

The interaction matters

fit$GEI_Note
#>   The term Treatment x Env (P < 2.2e-16) is significant at alpha = 0.05.
#>   Treatment rankings are therefore not consistent across environments, and
#>   the pooled treatment means in 'Treatments_Comparison' average over that
#>   inconsistency. Base any recommendation on the environment-wise comparisons
#>   in 'Location_wise' as well, and consider a stability or
#>   multiplicative-model analysis (for example AMMI or GGE) before advancing a
#>   treatment on its pooled mean alone.

The interaction is significant, so the pooled ranking conceals genotypes that excel in some environments and not others. The environment-wise comparisons are already in the object:

fit$Location_wise[["TCSM-IP"]]$groups[1:5, ]
#>   Treatment    Mean Groups
#> 1  Co 12009 19.8965      a
#> 2  Co 16018 19.0625     ab
#> 3  Co 12008 17.8655     bc
#> 4  Co 17012 17.7115     bc
#> 5  Co 14027 16.3740     cd

Designs with more than one error stratum

A split plot design has two error terms and a strip plot design three. This is where the choice of denominator is not a detail: the main plot effect of a split plot tested against the subplot error is simply the wrong test.

data(spd_example)
sp <- PooledSPD(spd_example, "Yield", "Location", "Replication", "MainPlot",
                "SubPlot", 0.05, 1)
sp$Pooled_ANOVA
#>  Source          Df Sum Sq   Mean Sq F value  Pr(>F)   Sig Error          
#>  Environment      1 433280.  433300. 17.66    0.005674 **  Rep within Env 
#>  Rep within Env   6 147240.  24540.  1.929    0.222068     Main plot error
#>  Main plot        1 502420.  502400. 39.48    0.000756 *** Main plot error
#>  Main plot x Env  1 46.021   46.02   0.003617 0.953998     Main plot error
#>  Main plot error  6 76349.   12720.  8.645    4.58e-05 *** Subplot error  
#>  Subplot          2 177930.  88970.  60.44    4.27e-10 *** Subplot error  
#>  Subplot x Env    2 2213.1   1107.   0.7518   0.482313     Subplot error  
#>  Main plot x SP   2 960.17   480.1   0.3262   0.724836     Subplot error  
#>  MP x SP x Env    2 1128.1   564.0   0.3832   0.685776     Subplot error  
#>  Subplot error   24 35326.   1472.                                        
#>  Total           47 1376900.                                              
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 'Error' names the mean square used as the denominator of each F-test.

Mean separation follows the same logic. Main plot means are compared with the main plot error, subplot means with the subplot error:

sp$Treatments_Comparison$MainPlot$statistics
#>    MSerror Df     Mean       CV Critical_difference
#> 1 12724.77  6 5144.129 2.192871             79.6806
sp$Treatments_Comparison$SubPlot$statistics
#>    MSerror Df     Mean        CV Critical_difference
#> 1 1471.905 24 5144.129 0.7458096            33.87375

The strip plot design returns three sets of comparisons, one per stratum:

data(stpd_example)
st <- PooledStPD(stpd_example, "Yield", "Location", "Replication", "RowFactor",
                 "ColumnFactor", 0.05, 1)
st$Pooled_ANOVA
#>  Source             Df Sum Sq  Mean Sq F value Pr(>F)   Sig Error             
#>  Environment         1 54493.  54490.  1.845   0.22323      Rep within Env    
#>  Rep within Env      6 177230. 29540.  15.28   5.46e-05 *** Intersection error
#>  Row factor          1 104430. 104400. 14.98   0.00826  **  Row strip error   
#>  Row x Env           1 8318.7  8319.   1.194   0.31653      Row strip error   
#>  Row strip error     6 41819.  6970.   3.605   0.02795  *   Intersection error
#>  Column factor       2 113960. 56980.  25.46   4.81e-05 *** Column strip error
#>  Column x Env        2 262.65  131.3   0.05868 0.94327      Column strip error
#>  Column strip error 12 26854.  2238.   1.158   0.40199      Intersection error
#>  Row x Column        2 429.41  214.7   0.1111  0.89579      Intersection error
#>  Row x Col x Env     2 1230.9  615.4   0.3184  0.73330      Intersection error
#>  Intersection error 12 23197.  1933.                                          
#>  Total              47 552220.                                                
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 'Error' names the mean square used as the denominator of each F-test.
names(st$Treatments_Comparison)
#> [1] "RowFactor"    "ColumnFactor" "Interaction"

Choosing a multiple comparison procedure

Mult_Comp_Test = 1 gives Tukey’s HSD, which controls the family-wise error rate over all pairwise comparisons and is the sensible default when all pairwise differences are of interest. 2 gives Duncan’s multiple range test, whose graduated protection levels buy power at the cost of weaker error control. 3 gives Fisher’s LSD, which makes no adjustment at all; with means it conducts independent tests at level , so the chance of at least one false declaration grows quickly. CANE warns when LSD is requested for more than three treatments. Reserve it for a small number of comparisons planned before the data were seen.

Fixed or random environments

By default environments are treated as fixed, which is appropriate when they are a specific set of target sites or seasons. If the environments are better regarded as a sample from a wider population to which the conclusion should generalise, the treatment main effect should be tested against the treatment-by-environment interaction:

fr <- PooledRBD(sugarcane, "ccsy", "env", "gen", "rep", 0.05, 1,
                env_effect = "random")
r <- fr$Pooled_ANOVA
r[r$Source %in% c("Treatment", "Treatment x Env"), ]
#>  Source          Df Sum Sq Mean Sq F value Pr(>F) Sig Error          
#>  Treatment       17 503.72 29.63   2.042   0.0171 *   Treatment x Env
#>  Treatment x Env 85 1233.5 14.51   14.51   <1e-16 *** Pooled error   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 'Error' names the mean square used as the denominator of each F-test.

The pooled comparison of means then uses the same mean square, so the grouping letters and the F-test rest on the same denominator:

fr$Treatments_Comparison$parameters[, c("test", "error_stratum", "scale")]
#>          test   error_stratum              scale
#> 1 Tukey's HSD Treatment x Env Aitken-transformed

What CANE is not for

CANE implements the fixed-effects pooled analysis of a balanced series of experiments. It is deliberately narrow. For

unbalanced series, or treatment sets that only partly overlap between environments, use a mixed-model package such as nlme or lme4, which handle missing cells properly;
augmented designs, in which a few checks are replicated and the test entries are not, use augmentedRCBD (archived on CRAN in April 2026 for an administrative reason; it remains installable from its own repository);
the single-trial stage of a mixed-model pipeline, use statgenSTA;
explicit modelling of heterogeneous variances or of a spatial covariance structure, use nlme, lme4, sommer or asreml;
multiplicative decomposition of the interaction and stability analysis (AMMI, GGE biplots, Shukla’s variance, Eberhart-Russell regression), use metan or ammistability;
best linear unbiased prediction of genotype effects for selection, use a mixed-model package.

The pooled ANOVA is the first step of a multi-environment analysis, not the last one. Its value is that it tests the main effects and the interaction exactly and transparently, which tells you whether the more elaborate analyses are worth running.