Single-Stage to Multi-Stage Trial Analysis

This document illustrates how Mandala connects a single-stage MET analysis with a two-stage workflow. The same breeding-analysis grammar is retained while the analysis moves from plot-level observations to environment-level genotype BLUEs and their stage-1 covariance matrix.

library(mandala)

path <- system.file("extdata", "fullrep_MET_n1000.csv", package = "mandala", mustWork = TRUE)
met <- read.csv(path, stringsAsFactors = FALSE, check.names = FALSE)
for (v in c("geno", "env", "rep", "block", "row", "col")) met[[v]] <- factor(met[[v]])
met$yld <- as.numeric(met$yld)

1. Single-Stage MET

fit_single <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + geno:env + env:rep:block,
  data   = met,
  verbose = FALSE
)

summary(fit_single)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + geno:env + env:rep:block, 
#>     data = met, verbose = FALSE)
#> 
#> Variance Components:
#>      component   estimate  std.error   z.ratio bound %ch
#>           geno 0.43593219 0.09472589  4.602038     P  NA
#>       geno:env 0.07591478 0.02730378  2.780376     P  NA
#>  env:rep:block 0.08022915 0.02319411  3.459032     P  NA
#>       R.sigma2 0.48914215 0.03214592 15.216304     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0907140 0.1699720 41.716954
#>     envY1_L2 -0.2670984 0.2008860 -1.329602
#>     envY1_L3 -0.7784500 0.2008736 -3.875323
#>     envY1_L4 -1.1050189 0.2008708 -5.501142
#>     envY1_L5 -0.4542703 0.2008850 -2.261345
#> 
#> Converged: TRUE  |  Iterations: 10
#> 
#> Model Notes:
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>  random level   estimate std.error    z.ratio
#>    geno    G1 -0.8617860 0.1967470 -4.3801743
#>    geno   G10  0.1274719 0.1967575  0.6478631
#>    geno   G11  1.3460707 0.1965576  6.8482268
#>    geno   G12  0.2550975 0.1971240  1.2940964
#>    geno   G13  1.0727427 0.1966296  5.4556516
#> 
#> logLik: -1235.493   AIC: 2478.987   BIC: 2498.577   logLik_Trunc: -325.744
pred_single <- mandala_predict(fit_single, "geno", verbose = FALSE)
head(pred_single)
#>   geno predicted_value std_error
#> 1   G1        5.773522 0.1590994
#> 2  G10        6.780007 0.1591189
#> 3  G11        8.019827 0.1590331
#> 4  G12        6.909855 0.1592230
#> 5  G13        7.741739 0.1590373
#> 6  G14        6.515366 0.1590225

2. Two-Stage Analysis

Two-stage analysis is useful when a large MET is naturally analyzed in pieces. In stage 1, each environment is analyzed separately to produce adjusted genotype BLUEs and their covariance matrix. In stage 2, those adjusted means are analyzed across environments. This can be efficient, transparent, and practical for large breeding trial networks.

2.1 First-Stage Analysis

First define the per-environment model. Here genotype is fixed in stage 1 so that Mandala can extract adjusted genotype BLUEs for each environment.

s1_prep <- stage1_prep(
  df = met,
  env_col = "env",
  s1_fixed = yld ~ geno,
  s1_random = ~ rep:block,
  s1_classify_term = "geno",
  response_var = "yld"
)

stage1_audit(met, s1_prep)
#> 
#> Stage-1 audit
#> -------------
#> Environments: 10 
#> Total rows: 1000 
#> Response: yld 
#> Classify term: geno 
#> Stage-1 fixed: yld ~ geno 
#> Stage-1 random: ~rep:block 
#> 
#>      env n_rows n_nonmissing_y n_classify_levels
#> 1  Y1_L1    100            100                50
#> 2  Y1_L2    100            100                50
#> 3  Y1_L3    100            100                50
#> 4  Y1_L4    100            100                50
#> 5  Y1_L5    100            100                50
#> 6  Y2_L1    100            100                50
#> 7  Y2_L2    100            100                50
#> 8  Y2_L3    100            100                50
#> 9  Y2_L4    100            100                50
#> 10 Y2_L5    100            100                50

Now fit the first-stage models quietly across environments. Each environment is fit separately, so quiet output is important when many environments are present.

s1 <- mandala_stage1(
  df = met,
  prep = s1_prep,
  verbose = FALSE
)

stage1_report_check(s1)
#> 
#> Stage-1 report check
#> --------------------
#> Successful environments: 10 
#> Failed environments    : 0 
#>      env n_rows n_BLUEs has_vcov vcov_dim
#> 1  Y1_L1    100      50     TRUE    50x50
#> 2  Y1_L2    100      50     TRUE    50x50
#> 3  Y1_L3    100      50     TRUE    50x50
#> 4  Y1_L4    100      50     TRUE    50x50
#> 5  Y1_L5    100      50     TRUE    50x50
#> 6  Y2_L1    100      50     TRUE    50x50
#> 7  Y2_L2    100      50     TRUE    50x50
#> 8  Y2_L3    100      50     TRUE    50x50
#> 9  Y2_L4    100      50     TRUE    50x50
#> 10 Y2_L5    100      50     TRUE    50x50

The stage-1 object stores one fitted model per environment. For example, inspect the first few genotype BLUEs carried forward from stage 1.

stage1_blues <- do.call(rbind, lapply(s1$results, function(x) {
  out <- x$BLUEs[, c("geno", "predicted_value", "std_error")]
  out$env <- x$env
  out[, c("env", "geno", "predicted_value", "std_error")]
}))

head(stage1_blues)
#>           env geno predicted_value std_error
#> Y1_L1.1 Y1_L1   G1        5.314903 0.5489310
#> Y1_L1.2 Y1_L1  G10        6.320838 0.5489249
#> Y1_L1.3 Y1_L1  G11        8.869945 0.5489375
#> Y1_L1.4 Y1_L1  G12        7.033992 0.5489310
#> Y1_L1.5 Y1_L1  G13        8.770005 0.5489310
#> Y1_L1.6 Y1_L1  G14        7.436263 0.5489310

Stage-1 heritability is a useful trial-quality summary. Mandala estimates it from separate per-environment models where genotype is fitted as random. Poor or unstable environments can then be reviewed before stage 2.

s1_h2 <- mandala.h2(
  df = met,
  out = s1,
  env_col = "env",
  genotype_col = "geno",
  response_col = "yld",
  h2_fixed = yld ~ 1,
  h2_random = ~ geno + rep:block,
  verbose = FALSE
)

head(s1_h2)
#>     env n_obs n_genotype  sigma_g2  sigma_e2 n_rep_arithmetic n_rep_harmonic  mean_PEV mean_PEVD   H2_plot H2_standard H2_entrymean_PEV H2_Cullis status
#> 1 Y1_L1   100         50 0.5449739 0.5996125                2              2 0.3013242 0.6016202 0.4761318   0.6451074        0.4470852 0.4480284     ok
#> 2 Y1_L2   100         50 0.6659502 0.4879891                2              2 0.2590033 0.5108521 0.5771102   0.7318578        0.6110771 0.6164488     ok
#> 3 Y1_L3   100         50 0.6405839 0.4140495                2              2 0.2588452 0.4507066 0.6073996   0.7557543        0.5959230 0.6482064     ok
#> 4 Y1_L4   100         50 0.4442142 0.4921793                2              2 0.2591018 0.4967718 0.4743884   0.6435053        0.4167188 0.4408420     ok
#> 5 Y1_L5   100         50 0.6177487 0.4315737                2              2 0.2285999 0.4481608 0.5887120   0.7411186        0.6299468 0.6372629     ok
#> 6 Y2_L1   100         50 0.4628072 0.5312939                2              2 0.2799074 0.5484042 0.4655534   0.6353278        0.3951964 0.4075241     ok

2.2 Second-Stage Preparation

mandala_stage2_prep() returns a data frame of stage-1 BLUEs and an RMAT covariance matrix aligned to that data frame.

s2 <- mandala_stage2_prep(s1)
stage2_df <- s2$data
RMAT <- s2$RMAT

head(stage2_df)
#>     env geno      yld SE_stage1
#> 1 Y1_L1   G1 5.314903 0.5489310
#> 2 Y1_L1  G10 6.320838 0.5489249
#> 3 Y1_L1  G11 8.869945 0.5489375
#> 4 Y1_L1  G12 7.033992 0.5489310
#> 5 Y1_L1  G13 8.770005 0.5489310
#> 6 Y1_L1  G14 7.436263 0.5489310
dim(RMAT)
#> [1] 500 500

The first rows and columns of RMAT show the covariance among the first-stage BLUEs passed into stage 2.

round(as.matrix(RMAT[1:5, 1:5]), 4)
#>           Y1_L1|G1 Y1_L1|G10 Y1_L1|G11 Y1_L1|G12 Y1_L1|G13
#> Y1_L1|G1    0.3013    0.0007    0.0000    0.0014    0.0014
#> Y1_L1|G10   0.0007    0.3013    0.0007    0.0007    0.0007
#> Y1_L1|G11   0.0000    0.0007    0.3013    0.0000    0.0000
#> Y1_L1|G12   0.0014    0.0007    0.0000    0.3013    0.0014
#> Y1_L1|G13   0.0014    0.0007    0.0000    0.0014    0.3013

2.3 Second-Stage Fit

The default and recommended stage-2 covariance syntax is `R_formula =
vcov(RMAT)`. 
fit_two <- mandala_stage2(
  data = stage2_df,
  RMAT = RMAT,
  fixed = yld ~ env,
  random = ~ geno + geno:env,
  R_formula = ~ vcov(RMAT),
  env = "env",
  id = "geno",
  verbose = FALSE
)

summary(fit_two)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + geno:env, data = data, 
#>     matrix_list = matrix_list, R_formula = ~vcov(RMAT), verbose = FALSE, 
#>     method = "sparse", engine = "mme")
#> 
#> Variance Components:
#>  component   estimate    std.error      z.ratio bound %ch
#>       geno 0.42911835 9.355504e-02 4.586801e+00     P  NA
#>   geno:env 0.04026488 2.097152e+06 1.919979e-08     P  NA
#>   R.sigma2 0.05086729 2.097152e+06 2.425541e-08     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0853885 0.1301334 54.447126
#>     envY1_L2 -0.2795267 0.1372235 -2.037018
#>     envY1_L3 -0.7529652 0.2190914 -3.436763
#>     envY1_L4 -1.1002689 0.1605382 -6.853625
#>     envY1_L5 -0.4612663 0.1383428 -3.334227
#> 
#> Converged: TRUE  |  Iterations: 11
#> 
#> Model Notes:
#> - Model size class: complex (n=500, q=550, MME dim=560).
#> 
#> Random Effects (BLUPs) [first 5]:
#>  random level   estimate std.error    z.ratio
#>    geno    G1 -0.8355748 0.1984764 -4.2099466
#>    geno   G10  0.1335397 0.1984965  0.6727559
#>    geno   G11  1.3643650 0.1981645  6.8850122
#>    geno   G12  0.2506781 0.1989430  1.2600499
#>    geno   G13  1.0502786 0.1983903  5.2940020
#> 
#> logLik: -511.109   AIC: 1028.217   BIC: 1040.801   logLik_Trunc: -60.829
pred_two <- mandala_predict(fit_two, "geno", verbose = FALSE)
head(pred_two)
#>   geno predicted_value std_error
#> 1   G1        5.805882 0.1717629
#> 2  G10        6.784090 0.1717769
#> 3  G11        8.026464 0.1716367
#> 4  G12        6.902327 0.1719026
#> 5  G13        7.709430 0.1717100
#> 6  G14        6.512140 0.1716185

2.4 Compare Single-Stage and Two-Stage Results

cmp <- merge(
  pred_single[, c("geno", "predicted_value")],
  pred_two[, c("geno", "predicted_value")],
  by = "geno",
  suffixes = c("_single", "_two")
)

r_single_two <- cor(
  cmp$predicted_value_single,
  cmp$predicted_value_two,
  use = "complete.obs"
)

r_single_two
#> [1] 0.9994991

plot(
  cmp$predicted_value_single,
  cmp$predicted_value_two,
  pch = 16,
  xlab = "Single-stage genotype mean",
  ylab = "Two-stage genotype mean",
  main = "Single-stage vs two-stage genotype means"
)
abline(0, 1, col = "gray40", lty = 2)
legend(
  "topleft",
  legend = paste0("r = ", round(r_single_two, 3)),
  bty = "n"
)

2.5 Reliability

rel_two <- stage2_reliability(fit_two, classify_term = "geno", preds = pred_two)
head(rel_two)
#>    geno predicted_value std_error        PEV reliability  accuracy weight_invPEV
#> 29  G35        6.147453 0.1716135 0.02945118   0.9313682 0.9650742      33.95450
#> 33  G39        6.835278 0.1716142 0.02945142   0.9313676 0.9650739      33.95422
#> 6   G14        6.512140 0.1716185 0.02945291   0.9313641 0.9650721      33.95250
#> 15  G22        6.966521 0.1716205 0.02945360   0.9313625 0.9650713      33.95171
#> 26  G32        6.159648 0.1716294 0.02945667   0.9313554 0.9650676      33.94817
#> 46  G50        6.150468 0.1716299 0.02945682   0.9313550 0.9650674      33.94800

stage2_variance_summary(fit_two, stage2_df, RMAT)
#>         Component   Variance        PVE
#> 1             env         NA         NA
#> 2  genetic_effect 0.42911835 0.57526136
#> 3 homogeneous_GxE 0.05086729 0.06819095
#> 4    Stage1_error 0.26596807 0.35654769

Selected Literature

Single-stage and two-stage analyses are both standard mixed-model strategies for multi-environment trial analysis. Single-stage models analyze plot-level data directly, while two-stage methods first estimate environment-specific genotype effects and then combine those estimates with appropriate weighting or covariance information. The literature below provides useful background on REML mixed models, MET covariance modeling, and efficient stage-wise analysis.