Introduction to Mandala

Mandala is an R package for fitting mixed models commonly used in plant breeding trials. A typical workflow is:

  1. prepare the data;
  2. fit a model with mandala();
  3. inspect summary();
  4. test fixed effects;
  5. predict genotype or treatment means;
  6. check diagnostics and runtime.

This introduction uses two small MET datasets shipped with Mandala:

Set MANDALA_RUN_VIGNETTES=true before knitting to run all examples.

library(mandala)

The mandala() Function

The main fitting function is mandala().

fit <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + env:rep,
  data   = df
)

The most commonly used arguments are:

Argument Purpose Common use
fixed Fixed-effect formula yld ~ env
random Random-effect formula ~ geno + env:rep
data Analysis data frame one row per plot/observation
matrix_list Optional relationship matrices genomic relationship matrix
R_formula Residual structure default is independent residuals
engine Computational engine usually "auto"
method Matrix method usually "sparse"
em_iters Optional EM warmup iterations default is 0
mme_trace_mode Trace strategy for MME fits "auto", "exact", or "adaptive"
return_mme_inv Store full MME inverse default "auto" is memory-safe
verbose Print iteration log TRUE for learning, FALSE for routine use

By default, Mandala starts directly with AI-REML. To run an EM warmup, use em_iters = 5 or another positive integer.

Here are a few model patterns:

# Genotype random effect
mandala(yld ~ 1, ~ geno, data = df)

# Multi-environment trial
mandala(yld ~ env, ~ geno + env:rep, data = df)

# Genotype-by-environment model
mandala(yld ~ env, ~ geno + geno:env + env:rep, data = df)

# Genomic model
mandala(yld ~ env, ~ GM(geno, GRM), data = df, matrix_list = G_list)

Some common random-effect terms in field-trial models are:

Term Typical meaning
geno Genotype main effect; commonly used for BLUPs.
env:rep Replication nested within environment.
env:rep:block Block nested within replication and environment.
geno:env Genotype-by-environment deviation.
env:row Field row effect within environment.
env:col Field column effect within environment.
GM(geno, GRM) Genomic relationship model for genotype.

Post-Model Fitting Functions

After mandala() returns a fitted model object, the usual next steps are summary, inference, diagnostics, and prediction.

Task Function or object Default behavior and common options
Model summary summary(fit) Prints the model statement, variance components, fixed effects, convergence status, and first random effects.
Fixed-effect tests mandala_fixed_tests(fit) Default is an incremental GLS test with Satterthwaite denominator df for small/medium models. Use type = "selected" for a scalable selected-covariance test. Common denDF choices are "satterthwaite", "residual", "containment", and "stratum".
Diagnostic plots mandala_diagnostic_plots(fit) Produces base-R fitted/residual diagnostics. Use response = "yld" when the response name should be shown.
BLUE extraction fit$BLUEs Data frame of fixed-effect estimates, standard errors, and z-ratios when available.
BLUP extraction fit$BLUPs Data frame of random-effect predictions, standard errors, and z-ratios when available.
Prediction mandala_predict(fit, "geno") Predicts marginal means for a classification term. Useful options include present = TRUE for observed combinations and verbose = FALSE for quiet output.
Heritability h2_estimates(fit, genotype = "geno") Computes genotype-centered heritability summaries when the model contains a genotype term.

The most common post-fit workflow is:

summary(fit)
mandala_fixed_tests(fit)
head(fit$BLUEs)
head(fit$BLUPs)
head(mandala_predict(fit, "geno", verbose = FALSE))
mandala_diagnostic_plots(fit, response = "yld")

Model Evaluation

Mandala stores common model evaluation quantities in the fitted object.

Quantity Where to find it Interpretation
Log-likelihood fit$logLik Larger values indicate better fit to the data, conditional on the model framework.
AIC fit$AIC Penalizes model complexity; smaller values are preferred when comparing candidate models fitted to the same response and data.
BIC fit$BIC Similar to AIC, but with a stronger penalty for model complexity; smaller values are preferred.

Use these values as quick model-comparison summaries, and use fixed-effect tests, prediction behavior, diagnostics, and biological sense before making a final model choice.

Data Preparation

Read the data, convert design variables to factors, and make sure the response is numeric.

path <- system.file("extdata", "fullrep_MET_n1000.csv",
                    package = "mandala", mustWork = TRUE)

df <- read.csv(path, stringsAsFactors = FALSE, check.names = FALSE)

df$geno  <- factor(df$geno)
df$env   <- factor(df$env)
df$rep   <- factor(df$rep)
df$block <- factor(df$block)
df$row   <- factor(df$row)
df$col   <- factor(df$col)
df$yld   <- as.numeric(df$yld)
df$disease <- as.numeric(df$disease)

head(df)
#>   geno   env loc year rep block row col variety design_class target_n plot_id    yield  disease   height   heading      yld family sire dam
#> 1  G37 Y1_L1  L1   Y1   1    B1   3   2     G37  fullrep_MET     1000       1 6.659879 22.71252 26.94052 112.92483 6.659879   F004   P1 P71
#> 2  G50 Y1_L1  L1   Y1   1    B3   7   1     G50  fullrep_MET     1000       2 6.544204 36.19650 25.16138 110.67029 6.544204   F005   P1 P46
#> 3  G30 Y1_L1  L1   Y1   1    B3   6   7     G30  fullrep_MET     1000       3 7.189111 24.51450 26.56731 107.27796 7.189111   F003   P1 P24
#> 4  G25 Y1_L1  L1   Y1   1    B2   5   5     G25  fullrep_MET     1000       4 9.500000 34.69108 32.11104  90.33920 9.500000   F003   P1 P24
#> 5  G36 Y1_L1  L1   Y1   1    B2   5   1     G36  fullrep_MET     1000       5 7.693915 39.46150 27.25189  95.58215 7.693915   F004   P1 P71
#> 6   G4 Y1_L1  L1   Y1   1    B1   1   8      G4  fullrep_MET     1000       6 6.290341 25.56863 27.99881  99.76546 6.290341   F001   P1 P32

For later examples, use a small helper to keep the code focused on the model.

read_example_data <- function(file, factors) {
  dat <- read.csv(
    system.file("extdata", file, package = "mandala", mustWork = TRUE),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  dat[factors] <- lapply(dat[factors], factor)
  dat$yld <- as.numeric(dat$yld)
  dat
}

trial_summary <- function(dat) {
  data.frame(
    rows = nrow(dat),
    genotypes = length(unique(dat$geno)),
    environments = if ("env" %in% names(dat)) length(unique(dat$env)) else NA_integer_,
    reps = if ("rep" %in% names(dat)) length(unique(dat$rep)) else NA_integer_,
    blocks = if ("block" %in% names(dat)) length(unique(dat$block)) else NA_integer_,
    rows_in_field = if ("row" %in% names(dat)) length(unique(dat$row)) else NA_integer_,
    cols_in_field = if ("col" %in% names(dat)) length(unique(dat$col)) else NA_integer_
  )
}

To keep the first example quick, use a small subset of the shipped data.

df_intro <- subset(df, env %in% levels(df$env)[1:5] &
                     geno %in% levels(df$geno)[1:50])
df_intro <- droplevels(df_intro)
trial_summary(df_intro)
#>   rows genotypes environments reps blocks rows_in_field cols_in_field
#> 1  500        50            5    2      3             7             8

For a quick first analysis, the smallest complete workflow is:

fit <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + env:rep,
  data   = df_intro
)

summary(fit)
mandala_fixed_tests(fit)
mandala_predict(fit, "geno", verbose = FALSE)

Model Run 1: Simple MET with Genotype Random

For a simple multi-environment trial, fit environment as fixed and genotype plus replication-within-environment as random.

This first run uses verbose = TRUE so that you can see the AI-REML iteration log. In routine analyses, use verbose = FALSE.

fit <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + env:rep,
  data   = df_intro,
  verbose = TRUE
)
#> Initial data rows: 500 
#> Final data rows after NA handling: 500 
#> EM run skipped. To run "n" EM iterations, use option [em_iters=5]
#> Starting AI-REML logLik = -661.4512
#>  Iter     LogLik   Sigma2  DF     wall    Step  Ridge       Restrained
#>     1  -660.9977     0.636  495  00:37:42  1.00  2.10e-06  ( 0 restrained)
#>     2  -660.4752     0.636  495  00:37:42  1.00  2.13e-06  ( 0 restrained)
#>     3  -659.8747     0.636  495  00:37:42  1.00  2.16e-06  ( 0 restrained)
#>     4  -659.3100     0.636  495  00:37:43  1.00  2.22e-06  ( 0 restrained)
#>     5  -658.7902     0.636  495  00:37:43  1.00  2.31e-06  ( 0 restrained)
#>     6  -658.3260     0.636  495  00:37:43  1.00  2.46e-06  ( 0 restrained)
#>     7  -657.9287     0.636  495  00:37:43  1.00  2.70e-06  ( 0 restrained)
#>     8  -657.6100     0.636  495  00:37:44  1.00  3.10e-06  ( 0 restrained)
#>     9  -657.3796     0.636  495  00:37:44  1.00  3.72e-06  ( 0 restrained)
#>    10  -657.2443     0.636  495  00:37:44  1.00  4.68e-06  ( 0 restrained)
#>    11  -657.2049     0.636  495  00:37:44  1.00  6.07e-06  ( 0 restrained)
#>    12  -657.2049     0.636  495  00:37:45    NA  7.83e-06  ( 0 restrained)
#> No further logLik improvement at iter 12, but logLik/step has plateaued. Declaring convergence.
#> Main AI-REML Loop finished [total elapsed = 3.04 sec]. Combining results.
summary(fit)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + env:rep, data = df_intro, 
#>     verbose = TRUE)
#> 
#> Variance Components:
#>  component   estimate  std.error   z.ratio bound %ch
#>       geno 0.48298533 0.11050286  4.370795     P  NA
#>    env:rep 0.02543052 0.02413931  1.053490     P  NA
#>   R.sigma2 0.63565365 0.04280714 14.849242     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0849820 0.1695037 41.798392
#>     envY1_L2 -0.2941118 0.1953038 -1.505919
#>     envY1_L3 -0.8354016 0.1953038 -4.277446
#>     envY1_L4 -1.1108598 0.1953038 -5.687854
#>     envY1_L5 -0.4772864 0.1953038 -2.443815
#> 
#> Converged: TRUE  |  Iterations: 12
#> 
#> 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.9013381 0.2543792 -3.5432849
#>    geno   G10  0.1133713 0.2543792  0.4456783
#>    geno   G11  1.4471950 0.2543792  5.6891242
#>    geno   G12  0.2758597 0.2543792  1.0844427
#>    geno   G13  0.9672133 0.2543792  3.8022493
#> 
#> logLik: -657.205   AIC: 1320.410   BIC: 1333.023   logLik_Trunc: -202.330

Use mandala_fixed_tests() for term-level tests of fixed effects. The default uses a Satterthwaite-style test for small and medium models, and Mandala uses safer selected-covariance paths for large memory-safe fits.

mandala_fixed_tests(fit)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df     denDF      F.inc       Wald           Pr status
#> (Intercept)  1 19.181286 3175.77253 3175.77253 8.903907e-23     ok
#> env          4  4.987437   10.08131   40.32523 1.311952e-02     ok
#> 
#> Note:
#> - Sequential GLS Wald tests using fitted full-model variance components.
#> - Denominator df use a Satterthwaite approximation.
#> - For classical split-plot or known-stratum designs, consider denDF = 'containment' or denDF = 'stratum'.

When genotype is random, genotype effects are extracted as BLUPs.

geno_blups <- subset(fit$BLUPs, random == "geno")
head(geno_blups)
#>   random level   estimate std.error    z.ratio
#> 1   geno    G1 -0.9013381 0.2543792 -3.5432849
#> 2   geno   G10  0.1133713 0.2543792  0.4456783
#> 3   geno   G11  1.4471950 0.2543792  5.6891242
#> 4   geno   G12  0.2758597 0.2543792  1.0844427
#> 5   geno   G13  0.9672133 0.2543792  3.8022493
#> 6   geno   G14  0.1304795 0.2543792  0.5129329

Prediction

Predict genotype means:

pred_geno <- mandala_predict(fit, "geno", verbose = FALSE)
head(pred_geno)
#>   geno predicted_value std_error
#> 1   G1        5.640112 0.2426179
#> 2  G10        6.654821 0.2426179
#> 3  G11        7.988645 0.2426179
#> 4  G12        6.817310 0.2426179
#> 5  G13        7.508663 0.2426179
#> 6  G14        6.671930 0.2426179

Predict genotype-by-environment means:

pred_ge <- mandala_predict(fit, "geno:env", present = TRUE, verbose = FALSE)
head(pred_ge)
#>   geno   env predicted_value std_error
#> 1   G1 Y1_L1        6.183644 0.2722515
#> 2   G1 Y1_L2        5.889532 0.2722515
#> 3   G1 Y1_L3        5.348242 0.2722515
#> 4   G1 Y1_L4        5.072784 0.2722515
#> 5   G1 Y1_L5        5.706357 0.2722515
#> 6  G10 Y1_L1        7.198353 0.2722515

Heritability and Diagnostics

h2_estimates(fit, genotype = "geno")
#>               variable             value
#> 1        source_random              TRUE
#> 2         source_fixed             FALSE
#> 3             sigma_g2 0.482985333960737
#> 4             sigma_e2 0.635653653489121
#> 5                n_rep                10
#> 6        mean_PEV_BLUP  0.06470879263804
#> 7           avsed_BLUE                NA
#> 8           vdBLUE_avg                NA
#> 9            H2_Cullis 0.866023276302422
#> 10           H2_Piepho                NA
#> 11             H2_Plot 0.431761577577222
#> 12         H2_Standard 0.883697220716757
#> 13   H2_Cullis_clipped 0.866023276302422
#> 14   H2_Piepho_clipped                NA
#> 15     H2_Plot_clipped 0.431761577577222
#> 16 H2_Standard_clipped 0.883697220716757
mandala_diagnostic_plots(fit, response = "yld")

Model Run 2: Add a Covariate

Covariates are added to the fixed-effect formula. Here disease severity is used as a simple numeric covariate.

fit_cov <- mandala(
  fixed  = yld ~ env + disease,
  random = ~ geno + env:rep,
  data   = df_intro,
  verbose = FALSE
)

summary(fit_cov)
#> Model statement:
#> mandala(fixed = yld ~ env + disease, random = ~geno + env:rep, 
#>     data = df_intro, verbose = FALSE)
#> 
#> Variance Components:
#>  component  estimate  std.error   z.ratio bound %ch
#>       geno 0.4910459 0.11246442  4.366234     P  NA
#>    env:rep 0.0257118 0.02432445  1.057035     P  NA
#>   R.sigma2 0.6323724 0.04263275 14.833020     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  6.6750082 0.3069750 21.744465
#>     envY1_L2 -0.3007941 0.1958994 -1.535452
#>     envY1_L3 -0.8180585 0.1961530 -4.170513
#>     envY1_L4 -1.1919477 0.2022651 -5.892997
#>     envY1_L5 -0.6358293 0.2193507 -2.898688
#> 
#> Converged: TRUE  |  Iterations: 12
#> 
#> 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.86833563 0.2553539 -3.4005179
#>    geno   G10  0.09945221 0.2545611  0.3906810
#>    geno   G11  1.43636965 0.2545666  5.6424128
#>    geno   G12  0.19707192 0.2591809  0.7603643
#>    geno   G13  1.01015779 0.2556567  3.9512281
#> 
#> logLik: -659.805   AIC: 1325.611   BIC: 1338.219   logLik_Trunc: -205.850
mandala_fixed_tests(fit_cov)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df      denDF       F.inc        Wald           Pr status
#> (Intercept)  1  19.181286 3133.269382 3133.269382 1.012462e-22     ok
#> env          4   4.987437   10.024631   40.098524 1.327895e-02     ok
#> disease      1 432.314715    2.576582    2.576582 1.091858e-01     ok
#> 
#> Note:
#> - Sequential GLS Wald tests using fitted full-model variance components.
#> - Denominator df use a Satterthwaite approximation.
#> - For classical split-plot or known-stratum designs, consider denDF = 'containment' or denDF = 'stratum'.

Compare the first two models with log-likelihood, AIC, and BIC.

Note: these two models have different fixed-effect terms, so this table is mainly for illustrating where to find model evaluation values. Interpret model comparisons carefully and only compare models that answer the same biological question.

data.frame(
  model  = c("env", "env + disease"),
  logLik = c(fit$logLik, fit_cov$logLik),
  AIC    = c(fit$AIC, fit_cov$AIC),
  BIC    = c(fit$BIC, fit_cov$BIC)
)
#>           model    logLik      AIC      BIC
#> 1           env -657.2049 1320.410 1333.023
#> 2 env + disease -659.8055 1325.611 1338.219

Model Run 3: Add Genotype-by-Environment Random Effects

The term geno:env models genotype-by-environment deviation.

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

summary(fit_gxe)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + geno:env + env:rep:block, 
#>     data = df_intro, verbose = FALSE)
#> 
#> Variance Components:
#>      component   estimate  std.error   z.ratio bound %ch
#>           geno 0.46644783 0.10938081  4.264439     P  NA
#>       geno:env 0.11537671 0.04366255  2.642464     P  NA
#>  env:rep:block 0.06791418 0.03022383  2.247041     P  NA
#>       R.sigma2 0.49054218 0.04542669 10.798545     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0916433 0.1680490 42.199859
#>     envY1_L2 -0.2720598 0.1945190 -1.398628
#>     envY1_L3 -0.7833986 0.1945045 -4.027663
#>     envY1_L4 -1.1027840 0.1945004 -5.669830
#>     envY1_L5 -0.4599871 0.1945197 -2.364732
#> 
#> 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.9547980 0.2678439 -3.5647550
#>    geno   G10  0.1172227 0.2678493  0.4376444
#>    geno   G11  1.4057557 0.2678604  5.2480906
#>    geno   G12  0.2333505 0.2680692  0.8704860
#>    geno   G13  1.0058777 0.2673602  3.7622561
#> 
#> logLik: -647.053   AIC: 1302.107   BIC: 1318.925   logLik_Trunc: -192.179
mandala_fixed_tests(fit_gxe)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df    denDF       F.inc       Wald           Pr status
#> (Intercept)  1 65.85139 3289.424282 3289.42428 6.104505e-58     ok
#> env          4 26.39067    9.836288   39.34515 5.239533e-05     ok
#> 
#> Note:
#> - Sequential GLS Wald tests using fitted full-model variance components.
#> - Denominator df use a Satterthwaite approximation.
#> - For classical split-plot or known-stratum designs, consider denDF = 'containment' or denDF = 'stratum'.

head(subset(fit_gxe$BLUPs, random == "geno"))
#>   random level   estimate std.error    z.ratio
#> 1   geno    G1 -0.9547980 0.2678439 -3.5647550
#> 2   geno   G10  0.1172227 0.2678493  0.4376444
#> 3   geno   G11  1.4057557 0.2678604  5.2480906
#> 4   geno   G12  0.2333505 0.2680692  0.8704860
#> 5   geno   G13  1.0058777 0.2673602  3.7622561
#> 6   geno   G14  0.1242625 0.2678037  0.4640061
head(subset(fit_gxe$BLUPs, random == "geno:env"))
#>      random     level    estimate std.error    z.ratio
#> 51 geno:env  G1.Y1_L1 -0.23676425 0.2932719 -0.8073198
#> 52 geno:env G10.Y1_L1 -0.29301479 0.2932568 -0.9991749
#> 53 geno:env G11.Y1_L1  0.09039697 0.2937413  0.3077435
#> 54 geno:env G12.Y1_L1 -0.06690822 0.2932999 -0.2281222
#> 55 geno:env G13.Y1_L1  0.24132560 0.2932461  0.8229455
#> 56 geno:env G14.Y1_L1  0.09668341 0.2932745  0.3296687

Model Run 4: Genotype as Fixed

If genotype is fitted as a fixed effect, genotype estimates are extracted from the BLUE table.

fit_fixed <- mandala(
  fixed  = yld ~ env + geno,
  random = ~ env:rep,
  data   = df_intro,
  verbose = FALSE
)

summary(fit_fixed)
#> Model statement:
#> mandala(fixed = yld ~ env + geno, random = ~env:rep, data = df_intro, 
#>     verbose = FALSE)
#> 
#> Variance Components:
#>  component   estimate  std.error  z.ratio bound %ch
#>    env:rep 0.02543052 0.02413931  1.05349     P  NA
#>   R.sigma2 0.63565365 0.04280714 14.84924     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  6.0650194 0.2852470 21.262342
#>     envY1_L2 -0.2941118 0.1953038 -1.505919
#>     envY1_L3 -0.8354016 0.1953038 -4.277446
#>     envY1_L4 -1.1108598 0.1953038 -5.687854
#>     envY1_L5 -0.4772864 0.1953038 -2.443815
#> 
#> Converged: TRUE  |  Iterations: 12
#> 
#> 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
#>  env:rep Y1_L1.1 0.039356958 0.1302044 0.3022705
#>  env:rep Y1_L2.1 0.079135566 0.1302044 0.6077794
#>  env:rep Y1_L3.1 0.162604016 0.1302044 1.2488363
#>  env:rep Y1_L4.1 0.005918547 0.1302044 0.0454558
#>  env:rep Y1_L5.1 0.090000712 0.1302044 0.6912262
#> 
#> logLik: -600.522   AIC: 1205.044   BIC: 1213.245   logLik_Trunc: -190.675
mandala_fixed_tests(fit_fixed)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df      denDF        F.inc        Wald           Pr status
#> (Intercept)  1   8.995518 11218.284577 11218.28458 3.067931e-15     ok
#> env          4   4.965140    10.081308    40.32523 1.327726e-02     ok
#> geno        49 441.000000     8.598247   421.31409 1.103116e-39     ok
#> 
#> Note:
#> - Sequential GLS Wald tests using fitted full-model variance components.
#> - Denominator df use a Satterthwaite approximation.
#> - For classical split-plot or known-stratum designs, consider denDF = 'containment' or denDF = 'stratum'.

geno_blues <- subset(fit_fixed$BLUEs, grepl("^geno", effect))
head(geno_blues)
#>          effect estimate std.error  z.ratio
#> genoG10 genoG10 1.148255  0.356554 3.220423
#> genoG11 genoG11 2.657622  0.356554 7.453632
#> genoG12 genoG12 1.332128  0.356554 3.736119
#> genoG13 genoG13 2.114470  0.356554 5.930295
#> genoG14 genoG14 1.167614  0.356554 3.274720
#> genoG15 genoG15 1.296277  0.356554 3.635569

Model Run 5: Sparse p-rep MET with Field Design Terms

In a sparse MET, not every genotype is tested in every environment. A smaller set of shared genotypes connects environments, while many entries are tested in only one or a few environments. In a p-rep sparse MET, some observed genotype-environment cells are replicated, while others are represented once. This is common when many new breeding lines must be screened with limited field space.

For sparse field trials, it is often useful to include field-design terms such as row and column within environment.

sparse_met <- read_example_data(
  "sparse_prep_MET_n1000.csv",
  factors = c("geno", "env", "rep", "block", "row", "col")
)
sparse_met <- subset(sparse_met, env %in% levels(sparse_met$env)[1:5] &
                       geno %in% levels(sparse_met$geno)[1:160])
sparse_met <- droplevels(sparse_met)

trial_summary(sparse_met)
#>   rows genotypes environments reps blocks rows_in_field cols_in_field
#> 1  160        88            5    2      4             8             9

Check how many genotypes are shared between each pair of environments.

env_geno <- table(sparse_met$geno, sparse_met$env) > 0
t(env_geno) %*% env_geno
#>        
#>         Y1_L1 Y1_L2 Y1_L3 Y1_L4 Y1_L5
#>   Y1_L1    21     5     5     5     5
#>   Y1_L2     5    20     5     5     5
#>   Y1_L3     5     5    21     6     5
#>   Y1_L4     5     5     6    25     5
#>   Y1_L5     5     5     5     5    22

First fit a sparse MET model with genotype and field-design terms.

fit_sparse <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + env:row + env:col,
  data   = sparse_met,
  verbose = TRUE
)
#> Initial data rows: 160 
#> Final data rows after NA handling: 160 
#> EM run skipped. To run "n" EM iterations, use option [em_iters=5]
#> Starting AI-REML logLik = -229.8085
#>  Iter     LogLik   Sigma2  DF     wall    Step  Ridge       Restrained
#>     1  -225.5826     0.480  155  00:38:02  1.00  4.40e-07  ( 0 restrained)
#>     2  -221.8156     0.407  155  00:38:02  1.00  5.72e-07  ( 0 restrained)
#>     3  -219.0474     0.392  155  00:38:02  1.00  7.68e-07  ( 0 restrained)
#>     4  -217.1245     0.393  155  00:38:02  1.00  9.61e-07  ( 0 restrained)
#>     5  -215.8019     0.395  155  00:38:02  1.00  1.17e-06  ( 0 restrained)
#>     6  -214.9646     0.398  155  00:38:02  1.00  1.43e-06  ( 0 restrained)
#>     7  -214.4960     0.400  155  00:38:02  1.00  1.72e-06  ( 0 restrained)
#>     8  -214.2872     0.402  155  00:38:02  1.00  2.04e-06  ( 0 restrained)
#>     9  -214.2321     0.404  155  00:38:02  1.00  2.38e-06  ( 0 restrained)
#>    10  -214.2245     0.403  155  00:38:02  1.00  2.56e-06  ( 0 restrained)
#>    11  -214.2242     0.403  155  00:38:02  1.00  2.64e-06  ( 0 restrained)
#>    12  -214.2242     0.403  155  00:38:02  1.00  2.66e-06  ( 0 restrained)
#>    13  -214.2242     0.403  155  00:38:02  1.00  2.66e-06  ( 0 restrained)
#> Converged at iter 13 by logLik plateau: delta_loglik=8.67e-07
#> Main AI-REML Loop finished [total elapsed = 0.17 sec]. Combining results.

summary(fit_sparse)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + env:row + env:col, 
#>     data = sparse_met, verbose = TRUE)
#> 
#> Variance Components:
#>  component   estimate  std.error   z.ratio bound %ch
#>       geno 0.56448895 0.14658704 3.8508789     P  NA
#>    env:row 0.05011499 0.05497765 0.9115522     P  NA
#>    env:col 0.03400427 0.05366276 0.6336662     P  NA
#>   R.sigma2 0.40298179 0.08914755 4.5203909     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect    estimate std.error    z.ratio
#>  (Intercept)  7.69696188 0.2086182 36.8949721
#>     envY1_L2  0.05901633 0.2774721  0.2126928
#>     envY1_L3 -0.63520767 0.2757905 -2.3032257
#>     envY1_L4 -0.80770644 0.2737804 -2.9501987
#>     envY1_L5 -0.11439023 0.2730971 -0.4188628
#> 
#> Converged: TRUE  |  Iterations: 13
#> 
#> 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 -1.3456035 0.2464350 -5.460277
#>    geno   G10 -0.5990870 0.2370225 -2.527553
#>    geno  G103  0.9278005 0.4239502  2.188466
#>    geno  G105 -0.4902806 0.5129145 -0.955872
#>    geno  G106 -0.9480363 0.5136796 -1.845579
#> 
#> logLik: -214.224   AIC: 436.448   BIC: 448.622   logLik_Trunc: -71.789
mandala_fixed_tests(fit_sparse)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df    denDF       F.inc       Wald           Pr status
#> (Intercept)  1 42.69698 4466.980210 4466.98021 7.624345e-45     ok
#> env          4 19.30548    4.156088   16.62435 1.360264e-02     ok
#> 
#> Note:
#> - Sequential GLS Wald tests using fitted full-model variance components.
#> - Denominator df use a Satterthwaite approximation.
#> - For classical split-plot or known-stratum designs, consider denDF = 'containment' or denDF = 'stratum'.
head(mandala_predict(fit_sparse, "geno", verbose = FALSE))
#>   geno predicted_value std_error
#> 1   G1        6.051701 0.2338620
#> 2  G10        6.798217 0.2225257
#> 3 G103        8.325105 0.4228712
#> 4 G105        6.907024 0.5136291
#> 5 G106        6.449268 0.5146811
#> 6 G108        7.729100 0.5151287

Now add a genotype-by-environment random term. For sparse data, this model is more ambitious because not every genotype is observed in every environment. Here verbose = TRUE is used so the iteration log shows what happens if a variance component reaches the boundary.

fit_sparse_gxe <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + geno:env + env:row + env:col,
  data   = sparse_met,
  verbose = TRUE
)
#> Initial data rows: 160 
#> Final data rows after NA handling: 160 
#> EM run skipped. To run "n" EM iterations, use option [em_iters=5]
#> Starting AI-REML logLik = -237.4903
#>  Iter     LogLik   Sigma2  DF     wall    Step  Ridge       Restrained
#>     1  -232.1732     0.484  155  00:38:03  1.00  1.36e-07  ( 0 restrained)
#>     2  -229.4271     0.490  155  00:38:03  1.00  1.98e-07  ( 0 restrained)
#>     3  -226.7808     0.471  155  00:38:03  1.00  2.27e-07  ( 0 restrained)
#>     4  -221.3603     0.468  155  00:38:03  1.00  2.93e-07  ( 0 restrained)
#>     5  -219.1754     0.431  155  00:38:03  1.00  4.62e-07  ( 0 restrained)
#>     6  -216.6474     0.435  155  00:38:03  1.00  6.58e-07  ( 0 restrained)
#>     7  -215.7104     0.420  155  00:38:03  1.00  9.39e-07  ( 0 restrained)
#>     8  -214.9268     0.420  155  00:38:03  1.00  1.24e-06  ( 0 restrained)
#>     9  -214.6426     0.418  155  00:38:03  1.00  1.61e-06  ( 0 restrained)
#>    10  -214.5001     0.417  155  00:38:03  1.00  1.86e-06  ( 0 restrained)
#>    11  -214.4293     0.417  155  00:38:03  1.00  2.01e-06  ( 0 restrained)
#>    12  -214.3762     0.417  155  00:38:03  1.00  2.06e-06  ( 0 restrained)
#>    13  -214.3414     0.416  155  00:38:03  1.00  2.10e-06  ( 0 restrained)
#>    14  -214.3173     0.416  155  00:38:03  1.00  2.13e-06  ( 0 restrained)
#>    15  -214.3005     0.416  155  00:38:03  1.00  2.15e-06  ( 0 restrained)
#>    16  -214.2886     0.416  155  00:38:03  1.00  2.17e-06  ( 0 restrained)
#>    17  -214.2802     0.416  155  00:38:03  1.00  2.19e-06  ( 0 restrained)
#>    18  -214.2740     0.416  155  00:38:03  1.00  2.20e-06  ( 0 restrained)
#>  Froze geno:env at 2.274000e-03
#>    19  -214.2696     0.416  155  00:38:03  1.00  2.21e-06  ( 1 restrained)
#>    20  -214.2309     0.402  155  00:38:03  1.00  2.37e-06  ( 1 restrained)
#>    21  -214.2309     0.402  155  00:38:03  1.00  2.55e-06  ( 1 restrained)
#>    22  -214.2309     0.402  155  00:38:03  1.00  2.54e-06  ( 1 restrained)
#> Converged at iter 22 by logLik plateau: delta_loglik=1.05e-07
#> Main AI-REML Loop finished [total elapsed = 0.10 sec]. Combining results.
#> Warning: Boundary restraint occurred in 4 AI iterations; max restrained = 1. Check model fit.

summary(fit_sparse_gxe)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + geno:env + env:row + 
#>     env:col, data = sparse_met, verbose = TRUE)
#> 
#> Variance Components:
#>  component   estimate   std.error      z.ratio bound %ch
#>       geno 0.56284409  0.14625926 3.848263e+00     P  NA
#>   geno:env 0.00227400 50.01304688 4.546814e-05     B  NA
#>    env:row 0.05011187  0.05842397 8.577279e-01     P  NA
#>    env:col 0.03424646  0.05482966 6.245973e-01     P  NA
#>   R.sigma2 0.40201088  0.08893421 4.520317e+00     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect    estimate std.error    z.ratio
#>  (Intercept)  7.69746915 0.2088724 36.8524966
#>     envY1_L2  0.05758943 0.2779669  0.2071809
#>     envY1_L3 -0.63663462 0.2762237 -2.3047788
#>     envY1_L4 -0.80741659 0.2742087 -2.9445331
#>     envY1_L5 -0.11506859 0.2735863 -0.4205934
#> 
#> Converged: TRUE  |  Iterations: 22
#> 
#> 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 -1.3440800 0.2470235 -5.441103
#>    geno   G10 -0.5990672 0.2376488 -2.520809
#>    geno  G103  0.9245769 0.4248471  2.176258
#>    geno  G105 -0.4886965 0.5130068 -0.952612
#>    geno  G106 -0.9456476 0.5137532 -1.840665
#> 
#> logLik: -214.231   AIC: 438.462   BIC: 453.679   logLik_Trunc: -71.795
mandala_fixed_tests(
  fit_sparse_gxe,
  type  = "selected",
  denDF = "satterthwaite"
)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: selected 
#> Denominator df: satterthwaite 
#> 
#>             Df    denDF       F.inc       Wald           Pr status
#> (Intercept)  1 21.00825 1358.106507 1358.10651 1.413877e-20     ok
#> env          4 18.09014    4.138436   16.55374 1.491738e-02     ok
#> 
#> Note:
#> - Selected coefficient covariance was used; no dense marginal V or full MME inverse is required.
#> - Boundary/fixed variance components are excluded from the selected Satterthwaite df propagation.

The two models can be compared as nested random-effect structures fitted to the same response and fixed-effect model.

data.frame(
  model  = c("without GxE", "with GxE"),
  logLik = c(fit_sparse$logLik, fit_sparse_gxe$logLik),
  AIC    = c(fit_sparse$AIC, fit_sparse_gxe$AIC),
  BIC    = c(fit_sparse$BIC, fit_sparse_gxe$BIC)
)
#>         model    logLik      AIC      BIC
#> 1 without GxE -214.2242 436.4484 448.6221
#> 2    with GxE -214.2309 438.4618 453.6789

Takeaway

The basic Mandala workflow is:

fit <- mandala(yld ~ env, ~ geno + env:rep, data = df)
summary(fit)
mandala_fixed_tests(fit)
mandala_predict(fit, "geno")
h2_estimates(fit, genotype = "geno")
mandala_diagnostic_plots(fit)

Selected Literature

Mandala’s core workflow follows the linear mixed-model tradition used in plant and animal breeding. Restricted maximum likelihood was introduced by Patterson and Thompson (1971) for variance-component estimation in designed experiments. Henderson’s mixed-model equations provide the classical basis for BLUEs and BLUPs in mixed models. Average-information REML, described by Gilmour, Thompson, and Cullis (1995), is a widely used computational approach for variance-parameter estimation. Spatial and multi-environment breeding-trial analyses are closely related to the mixed-model framework described by Smith, Cullis, and Thompson (2001). Heritability summaries for unbalanced plant breeding trials should be interpreted with attention to prediction error and selection response, as discussed by Piepho and Moehring (2007).