Mandala with Advanced Variance Model Structures

This document shows a small set of representative advanced variance structures in Mandala. The goal is not to list every possible structure. Instead, the examples focus on model families that clarify how covariance structures change the biological interpretation of a mixed model.

The genomic prediction, spatial-analysis, and stagewise-analysis vignettes cover GM(), pspline2D(), and vcov(RMAT) in more detail. The examples focus on:

library(mandala)

G-Side and R-Side Structures

The same statistical idea can often appear in two different places in a mixed model:

Location Example Interpretation
G side random = ~ geno:diag(env) genetic variance differs by environment
R side R_formula = ~ diag(env) residual/error variance differs by environment
G side random = ~ geno:comsym(env) genotype-by-environment effects are equally correlated
R side R_formula = ~ comsym(env):id(geno) residual errors are equally correlated across environments within genotype
G side random = ~ FA(geno, env, k = 1) genotype-by-environment covariance is factor analytic

That distinction is important. A G-side structure describes covariance among random effects. An R-side structure describes covariance or heterogeneity in the residual errors.

Models Shown

The examples use the same fixed-effect part, yld ~ env, so that the variance structures are easy to compare.

Model Syntax Highlight Main Question
Reference MET geno + geno:env What is the ordinary baseline model?
G-side diagonal geno + geno:diag(env) Does genetic variance differ by environment?
R-side diagonal R_formula = ~ diag(env) Does residual noise differ by environment?
Compound symmetry geno + geno:comsym(env) Is one common GxE correlation adequate?
R-side compound symmetry R_formula = ~ comsym(env):id(geno) Are stage-two residuals correlated across environments?
FA without genotype main effect FA(geno, env, k = 1) Can FA alone describe the genotype-environment signal?
FA with genotype main effect geno + FA(geno, env, k = 1) What changes when a genotype main-effect term is included?
FA + GM FA(geno, env, k = 1) + GM(geno, GRM) Can FA GxE be combined with genomic main effects?
Environment-specific spatial residuals R_formula = ~ by(env):ar1(row):ar1(col) Does each environment need its own row-column residual pattern?
Unstructured GxE geno + us(env):geno What does the most flexible GxE covariance look like?

Data

The examples use the small full-rep MET dataset shipped with Mandala. The dataset has genotypes tested across environments with replication and field layout factors.

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)

data.frame(
  n_obs = nrow(met),
  n_geno = nlevels(met$geno),
  n_env = nlevels(met$env),
  n_rep = nlevels(met$rep),
  n_block = nlevels(met$block)
)
#>   n_obs n_geno n_env n_rep n_block
#> 1  1000     50    10     2       3

The FA+GM example uses a genomic relationship matrix.

GRM <- readRDS(system.file("extdata", "sim_GRM_1000.rds",
                           package = "mandala", mustWork = TRUE))

To keep this vignette fast, the examples below use a smaller balanced subset of the same dataset. The syntax is unchanged for the full dataset.

demo_geno <- levels(met$geno)[1:20]
demo_env <- levels(met$env)[1:5]

met_demo <- droplevels(subset(met, geno %in% demo_geno & env %in% demo_env))

G_list_demo <- mandala_grm_prep(
  GRM = GRM,
  data = met_demo,
  geno_col = "geno",
  GRM_label = "GRM",
  verbose = FALSE
)

data.frame(
  n_obs = nrow(met_demo),
  n_geno = nlevels(met_demo$geno),
  n_env = nlevels(met_demo$env),
  n_rep = nlevels(met_demo$rep),
  n_block = nlevels(met_demo$block)
)
#>   n_obs n_geno n_env n_rep n_block
#> 1   200     20     5     2       3

For compact output, this helper keeps only the variance-component table.

show_vc <- function(fit, n = 8) {
  head(fit$varcomp, n)
}

Reference Model

Start with a familiar MET model. This is the baseline for the comparisons below. The baseline random-effect part is:

~ geno + geno:env + env:rep:block

The advanced examples keep this design logic where possible. Some structures, such as FA() and us(), replace the ordinary geno:env term because they are alternative covariance models for the same genotype-by-environment effect.

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

summary(fit_ref)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~geno + geno:env + env:rep:block, 
#>     data = met_demo, verbose = FALSE)
#> 
#> Variance Components:
#>      component  estimate  std.error  z.ratio bound %ch
#>           geno 0.4223813 0.16709007 2.527866     P  NA
#>       geno:env 0.1699904 0.08413409 2.020470     P  NA
#>  env:rep:block 0.1018692 0.05978506 1.703925     P  NA
#>       R.sigma2 0.5304641 0.08243951 6.434586     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect    estimate std.error     z.ratio
#>  (Intercept)  7.25704646 0.2481438 29.24532412
#>     envY1_L2 -0.01552629 0.2850537 -0.05446794
#>     envY1_L3 -0.72707803 0.2839072 -2.56097076
#>     envY1_L4 -0.93492962 0.2836142 -3.29648362
#>     envY1_L5 -0.26793694 0.2822549 -0.94927294
#> 
#> Converged: TRUE  |  Iterations: 9
#> 
#> 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.17596379 0.3028640 -3.8828111
#>    geno   G10 -0.12230449 0.3040986 -0.4021870
#>    geno   G11  1.11534227 0.3040685  3.6680623
#>    geno   G12 -0.05260894 0.3030515 -0.1735974
#>    geno   G13  0.72632632 0.3020293  2.4048208
#> 
#> logLik: -272.866   AIC: 553.732   BIC: 566.824   logLik_Trunc: -93.673
mandala_fixed_tests(fit_ref)
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: incremental_gls 
#> Denominator df: satterthwaite 
#> 
#>             Df    denDF       F.inc       Wald           Pr status
#> (Intercept)  1 26.94936 1617.748590 1617.74859 1.337913e-25     ok
#> env          4 24.46478    4.457298   17.82919 7.607607e-03     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_ref, "geno", verbose = FALSE))
#>   geno predicted_value std_error
#> 1   G1        5.597333 0.2326623
#> 2  G10        6.735803 0.2328194
#> 3  G11        8.073070 0.2328310
#> 4  G12        6.811109 0.2325921
#> 5  G13        7.652742 0.2323349
#> 6  G14        6.694200 0.2326769

G-Side Diagonal: Heterogeneous Genetic Variance

The model below allows the genetic variance to differ by environment. This is a G-side model because the heterogeneity belongs to the random genotype effects. The genotype main-effect term is retained.

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

show_vc(fit_g_diag, n = 12)
#>             component     estimate  std.error      z.ratio bound %ch
#> 1                geno 4.309310e-01 0.16811978 2.5632377308     P  NA
#> 2 geno_diag_env_Y1_L1 2.010610e-01 0.18849733 1.0666516470     P  NA
#> 3 geno_diag_env_Y1_L2 1.338838e-01 0.16796812 0.7970786249     P  NA
#> 4 geno_diag_env_Y1_L3 2.131463e-01 0.19232421 1.1082656368     P  NA
#> 5 geno_diag_env_Y1_L4 3.402035e-05 0.12935240 0.0002630052     P  NA
#> 6 geno_diag_env_Y1_L5 3.107569e-01 0.22303197 1.3933291577     P  NA
#> 7       env:rep:block 1.039429e-01 0.05999971 1.7323898126     P  NA
#> 8            R.sigma2 5.286062e-01 0.08212263 6.4367909209     P  NA

Use this when the question is:

Do genotypes express different amounts of genetic variation in different environments?

R-Side Diagonal: Heterogeneous Residual Variance

Now keep the ordinary genotype and genotype-by-environment random effects, but allow the residual variance to differ by environment.

fit_r_diag <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + geno:env + env:rep:block,
  R_formula = ~ diag(env),
  data   = met_demo,
  verbose = FALSE
)

show_vc(fit_r_diag, n = 14)
#>         component  estimate  std.error  z.ratio bound %ch
#> 1            geno 0.4231353 0.16676978 2.537242     P  NA
#> 2        geno:env 0.1624235 0.08296627 1.957705     P  NA
#> 3   env:rep:block 0.1084528 0.06311124 1.718439     P  NA
#> 4 R.var.env.Y1_L1 0.6437743 0.19120299 3.366968     P  NA
#> 5 R.var.env.Y1_L2 0.4766374 0.14946400 3.188978     P  NA
#> 6 R.var.env.Y1_L3 0.4835264 0.15699128 3.079957     P  NA
#> 7 R.var.env.Y1_L4 0.4673979 0.14316156 3.264828     P  NA
#> 8 R.var.env.Y1_L5 0.5913868 0.19153861 3.087559     P  NA

Use this when the question is:

Are some environments noisier or less precise than others?

The syntax is similar to the previous model, but the meaning is different:

random    = ~ geno:diag(env)  # G-side genetic heterogeneity
R_formula = ~ diag(env)       # R-side residual heterogeneity

R-Side Compound Symmetry

R-side compound symmetry is most useful for stage-two means or other data with one residual row per unit-by-environment cell. In the example below, genotype-environment means are formed from the demonstration data, then residuals for the same genotype are allowed to be correlated across environments.

met_means <- aggregate(yld ~ geno + env, data = met_demo, FUN = mean)
met_means$geno <- factor(met_means$geno)
met_means$env <- factor(met_means$env)

fit_r_cs <- mandala(
  fixed  = yld ~ env,
  random = ~ geno,
  R_formula = ~ comsym(env):id(geno),
  data   = met_means,
  verbose = FALSE
)

show_vc(fit_r_cs)
#>                     component  estimate std.error     z.ratio bound %ch
#> 1                        geno 0.3292411  43.58565 0.007553888     P  NA
#> 2 R.comsym.sigma2.env.by.geno 0.5473485  43.58554 0.012558031     P  NA
#> 3    R.comsym.rho.env.by.geno 0.1312822  71.29087 0.001841501     P  NA

The alias corrv() is also accepted on the R side:

fit_r_cs2 <- mandala(
  fixed  = yld ~ env,
  random = ~ geno,
  R_formula = ~ corrv(env):id(geno),
  data   = met_means,
  verbose = FALSE
)

Compound-Symmetry GxE

Compound symmetry is a simple correlated GxE structure. It assumes equal variance across environments and one common correlation among environment effects.

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

show_vc(fit_cs)
#>                     component   estimate  std.error   z.ratio bound %ch
#> 1                        geno  0.4132874 0.16366440  2.525212     P  NA
#> 2 G.comsym.sigma2.env.by.geno  0.1582592 0.08262142  1.915474     P  NA
#> 3    G.comsym.rho.env.by.geno -0.2499990 0.24633787 -1.014862     B  NA
#> 4               env:rep:block  0.1063228 0.06072490  1.750892     P  NA
#> 5                    R.sigma2  0.5293953 0.08227490  6.434469     P  NA

The alias corrv() is also available:

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

Use compound symmetry when the covariance pattern should be simple and interpretable: one variance and one correlation.

Factor-Analytic GxE

Factor-analytic models are more flexible than compound symmetry. They model the genotype-by-environment covariance using latent factors.

First fit FA without a separate genotype main-effect term.

fit_fa1_no_geno <- mandala(
  fixed  = yld ~ env,
  random = ~ FA(geno, env, k = 1) + env:rep:block,
  data   = met_demo,
  mme_trace_mode = "adaptive",
  verbose = FALSE
)

fit_fa1_no_geno$varcomp
#>             component     estimate  std.error  z.ratio bound %ch
#> 1  FA_lambda_Y1_L1_F1 8.649485e-01         NA       NA     P  NA
#> 2  FA_lambda_Y1_L2_F1 8.798488e-01         NA       NA     P  NA
#> 3  FA_lambda_Y1_L3_F1 6.234405e-01         NA       NA     P  NA
#> 4  FA_lambda_Y1_L4_F1 5.568106e-01         NA       NA     P  NA
#> 5  FA_lambda_Y1_L5_F1 3.641782e-01         NA       NA     P  NA
#> 6        FA_psi_Y1_L1 8.881221e-02         NA       NA     P  NA
#> 7        FA_psi_Y1_L2 2.721628e-05         NA       NA     P  NA
#> 8        FA_psi_Y1_L3 2.110495e-01         NA       NA     P  NA
#> 9        FA_psi_Y1_L4 3.664743e-02         NA       NA     P  NA
#> 10       FA_psi_Y1_L5 3.116913e-01         NA       NA     P  NA
#> 12      env:rep:block 1.035645e-01 0.05909276 1.752575     P  NA
#> 13           R.sigma2 5.278975e-01 0.06835186 7.723235     P  NA

Now include the genotype main-effect term explicitly.

fit_fa1 <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + FA(geno, env, k = 1) + env:rep:block,
  data   = met_demo,
  mme_trace_mode = "adaptive",
  verbose = FALSE
)

fit_fa1$varcomp
#>             component  estimate  std.error  z.ratio bound %ch
#> 1                geno 1.1604153 0.63450707 1.828845     P  NA
#> 2  FA_lambda_Y1_L1_F1 0.1179992         NA       NA     P  NA
#> 3  FA_lambda_Y1_L2_F1 0.1238991         NA       NA     P  NA
#> 4  FA_lambda_Y1_L3_F1 0.1297991         NA       NA     P  NA
#> 5  FA_lambda_Y1_L4_F1 0.1356990         NA       NA     P  NA
#> 6  FA_lambda_Y1_L5_F1 0.1415990         NA       NA     P  NA
#> 7        FA_psi_Y1_L1 0.1524021         NA       NA     P  NA
#> 8        FA_psi_Y1_L2 0.1524021         NA       NA     P  NA
#> 9        FA_psi_Y1_L3 0.1892380         NA       NA     P  NA
#> 10       FA_psi_Y1_L4 0.1524021         NA       NA     P  NA
#> 11       FA_psi_Y1_L5 0.3257428         NA       NA     P  NA
#> 13      env:rep:block 0.1905026 0.10425039 1.827356     P  NA
#> 14           R.sigma2 0.5415233 0.07634461 7.093143     P  NA

A two-factor FA model can be useful when environment similarity is not well-described by a single latent axis.

fit_fa2 <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + FA(geno, env, k = 2) + env:rep:block,
  data   = met_demo,
  mme_trace_mode = "adaptive",
  verbose = FALSE
)

fit_fa2$varcomp
#>             component  estimate  std.error  z.ratio bound %ch
#> 1                geno 0.9283322 0.52435228 1.770436     P  NA
#> 2  FA_lambda_Y1_L1_F1 0.1573322         NA       NA     P  NA
#> 3  FA_lambda_Y1_L2_F1 0.1612655         NA       NA     P  NA
#> 4  FA_lambda_Y1_L3_F1 0.1651988         NA       NA     P  NA
#> 5  FA_lambda_Y1_L4_F1 0.1691321         NA       NA     P  NA
#> 6  FA_lambda_Y1_L5_F1 0.1730654         NA       NA     P  NA
#> 7  FA_lambda_Y1_L2_F2 0.1769987         NA       NA     P  NA
#> 8  FA_lambda_Y1_L3_F2 0.1809321         NA       NA     P  NA
#> 9  FA_lambda_Y1_L4_F2 0.1848654         NA       NA     P  NA
#> 10 FA_lambda_Y1_L5_F2 0.1887987         NA       NA     P  NA
#> 11       FA_psi_Y1_L1 0.7426658         NA       NA     P  NA
#> 12       FA_psi_Y1_L2 0.2032027         NA       NA     P  NA
#> 13       FA_psi_Y1_L3 0.2032027         NA       NA     P  NA
#> 14       FA_psi_Y1_L4 0.2032027         NA       NA     P  NA
#> 15       FA_psi_Y1_L5 0.2032027         NA       NA     P  NA
#> 17      env:rep:block 0.2540034 0.13918755 1.824901     P  NA
#> 18           R.sigma2 0.5426559 0.07892638 6.875469     P  NA

FA tools can summarize the fitted environment covariance, environment correlation, and genotype scores.

fa <- mandala_fa_tools(fit_fa2, geno = "geno", env = "env", rotate = TRUE)
fa$Lambda
#> 
#> Loadings:
#>       F1    F2   
#> Y1_L1 0.157      
#> Y1_L2 0.160 0.178
#> Y1_L3 0.164 0.182
#> Y1_L4 0.168 0.186
#> Y1_L5 0.172 0.190
#> 
#>                   F1    F2
#> SS loadings    0.135 0.135
#> Proportion Var 0.027 0.027
#> Cumulative Var 0.027 0.054
fa$VCOV
#>            Y1_L1      Y1_L2      Y1_L3      Y1_L4      Y1_L5
#> Y1_L1 0.76741922 0.02537226 0.02599110 0.02660994 0.02722877
#> Y1_L2 0.02537226 0.26053787 0.05866562 0.05999612 0.06132662
#> Y1_L3 0.02599110 0.05866562 0.26322981 0.06138850 0.06274994
#> Y1_L4 0.02660994 0.05999612 0.06138850 0.26598363 0.06417326
#> Y1_L5 0.02722877 0.06132662 0.06274994 0.06417326 0.26879933
fa$Cor
#>            Y1_L1      Y1_L2      Y1_L3      Y1_L4      Y1_L5
#> Y1_L1 1.00000000 0.05674237 0.05782835 0.05889793 0.05995116
#> Y1_L2 0.05674237 1.00000000 0.22401687 0.22790838 0.23173920
#> Y1_L3 0.05782835 0.22401687 1.00000000 0.23200217 0.23590204
#> Y1_L4 0.05889793 0.22790838 0.23200217 1.00000000 0.24000074
#> Y1_L5 0.05995116 0.23173920 0.23590204 0.24000074 1.00000000
head(fa$geno_summary)
#>     geno genomic_main common_GE_mean total_GE_mean predicted_common_mean predicted_total_mean common_stability_rms specific_rms total_GE_rms        F1        F2
#> G16  G16            0     0.16799691    0.17111600              7.048659             7.051778           0.21586172    0.1894068    0.3337986  3.829695 -3.133306
#> G11  G11            0     0.15373549    0.15063826              7.034397             7.031300           0.08820925    0.1055959    0.2053496  2.105316 -1.305680
#> G25  G25            0     0.13963472    0.14371759              7.020296             7.024379           0.43504338    0.2267208    0.5116307  6.453839 -6.254532
#> G13  G13            0     0.13123754    0.13342662              7.011899             7.014088           0.22348260    0.1702730    0.3108904  3.693607 -3.231038
#> G15  G15            0     0.06449388    0.06642769              6.945156             6.947089           0.33747627    0.1812258    0.3891638  4.727800 -4.838450
#> G18  G18            0     0.02441086    0.02366545              6.905073             6.904327           0.25226783    0.1704133    0.3056384 -3.072262  3.594511
#>     rank_predicted_common rank_predicted_total
#> G16                     1                    1
#> G11                     2                    2
#> G25                     3                    3
#> G13                     4                    4
#> G15                     5                    5
#> G18                     6                    6

The fitted FA correlation matrix can also be used to cluster environments.

plot_mandala_fa_env_similarity(
  fa,
  use = "correlation",
  main = "FA environment clustering"
)

FA + GM

FA can also be combined with a genomic main-effect term. In this model, GM(geno, GRM) is the genotype main-effect covariance and FA(geno, env, k = 1) models genotype-by-environment covariance.

fit_fa_gm <- mandala(
  fixed  = yld ~ env,
  random = ~ FA(geno, env, k = 1) + GM(geno, GRM) + env:rep:block,
  data   = met_demo,
  matrix_list = G_list_demo,
  mme_trace_mode = "adaptive",
  verbose = FALSE
)

fit_fa_gm$varcomp
#>             component   estimate  std.error  z.ratio bound %ch
#> 1  FA_lambda_Y1_L1_F1 0.54716741         NA       NA     P  NA
#> 2  FA_lambda_Y1_L2_F1 0.75874115         NA       NA     P  NA
#> 3  FA_lambda_Y1_L3_F1 0.51611237         NA       NA     P  NA
#> 4  FA_lambda_Y1_L4_F1 0.21264660         NA       NA     P  NA
#> 5  FA_lambda_Y1_L5_F1 0.07964944         NA       NA     P  NA
#> 6        FA_psi_Y1_L1 0.18117252         NA       NA     P  NA
#> 7        FA_psi_Y1_L2 0.08572616         NA       NA     P  NA
#> 8        FA_psi_Y1_L3 0.15771705         NA       NA     P  NA
#> 9        FA_psi_Y1_L4 0.08572616         NA       NA     P  NA
#> 10       FA_psi_Y1_L5 0.25323264         NA       NA     P  NA
#> 12       GM(geno,GRM) 0.31423112 0.20751684 1.514244     P  NA
#> 13      env:rep:block 0.10715770 0.06596431 1.624480     P  NA
#> 14           R.sigma2 0.54250971 0.07599147 7.139087     P  NA

In practice, avoid adding both geno and GM(geno, GRM) as competing genotype main-effect terms unless that comparison is intentional. GM(geno, GRM) is the genomic version of the genotype main-effect covariance.

Environment-Specific Spatial Residuals

Spatial field trend is often most naturally represented on the residual side. The by(env):ar1(row):ar1(col) residual structure fits a separate row-column AR1 residual pattern within each environment, while keeping genotype and GxE on the G side.

fit_spatial_r <- mandala(
  fixed  = yld ~ env,
  random = ~ geno + geno:env,
  R_formula = ~ by(env):ar1(row):ar1(col),
  data   = met_demo,
  verbose = FALSE
)

show_vc(fit_spatial_r, n = 12)
#>              component   estimate std.error z.ratio bound %ch
#> 1                 geno  0.2990856        NA      NA     P  NA
#> 2             geno:env  0.8363919        NA      NA     P  NA
#> 3   R.sigma2.env.Y1_L1  4.6828767        NA      NA     P  NA
#> 4  R.rho.row.env.Y1_L1  0.9430000        NA      NA     P  NA
#> 5  R.rho.col.env.Y1_L1 -0.5430000        NA      NA     P  NA
#> 6   R.sigma2.env.Y1_L2  0.6617576        NA      NA     P  NA
#> 7  R.rho.row.env.Y1_L2 -0.4383962        NA      NA     P  NA
#> 8  R.rho.col.env.Y1_L2  0.4717172        NA      NA     P  NA
#> 9   R.sigma2.env.Y1_L3  1.2677645        NA      NA     P  NA
#> 10 R.rho.row.env.Y1_L3  0.2358411        NA      NA     P  NA
#> 11 R.rho.col.env.Y1_L3  0.3804420        NA      NA     P  NA
#> 12  R.sigma2.env.Y1_L4  0.7553618        NA      NA     P  NA

This type of residual structure is useful for environment-specific spatial patterns. It should be added only when the design or diagnostics suggest that independent residual variation is too simple.

Unstructured GxE

An unstructured covariance is the most flexible GxE model shown here. It estimates all environment variances and covariances for the genotype-by- environment effect.

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

show_vc(fit_us, n = 12)
#>                       component    estimate std.error    z.ratio bound %ch
#> 1                          geno  0.41656001 0.1719842  2.4220823     P  NA
#> 2  G.us.env.by.geno.Y1_L1.Y1_L1  0.10447875 0.1959643  0.5331520     P  NA
#> 3  G.us.env.by.geno.Y1_L2.Y1_L1 -0.15040609 0.1380002 -1.0898974     P  NA
#> 4  G.us.env.by.geno.Y1_L2.Y1_L2  0.10447875 0.1948188  0.5362869     P  NA
#> 5  G.us.env.by.geno.Y1_L3.Y1_L1  0.13718616 0.1552147  0.8838475     P  NA
#> 6  G.us.env.by.geno.Y1_L3.Y1_L2 -0.09930617 0.1513577 -0.6561027     P  NA
#> 7  G.us.env.by.geno.Y1_L3.Y1_L3  0.23820619 0.2418446  0.9849555     P  NA
#> 8  G.us.env.by.geno.Y1_L4.Y1_L1 -0.03460073 0.1347036 -0.2568657     P  NA
#> 9  G.us.env.by.geno.Y1_L4.Y1_L2  0.06445315 0.1359603  0.4740585     P  NA
#> 10 G.us.env.by.geno.Y1_L4.Y1_L3 -0.04295336 0.1504928 -0.2854181     P  NA
#> 11 G.us.env.by.geno.Y1_L4.Y1_L4  0.10447875 0.1954488  0.5345583     P  NA
#> 12 G.us.env.by.geno.Y1_L5.Y1_L1 -0.05951104 0.1475108 -0.4034350     P  NA

Use us(env):geno only when the number of environments is modest. The number of covariance parameters grows quickly as the number of environments increases.

Compare the Models

Information criteria are useful for rough comparison when the models are fitted to the same response and fixed-effect structure.

model_compare <- data.frame(
  model = c(
    "reference",
    "G-side diag(env)",
    "R-side diag(env)",
    "compound symmetry",
    "FA k=1 without geno",
    "FA k=1",
    "FA k=2",
    "FA + GM",
    "R-side AR1xAR1 by env",
    "US GxE"
  ),
  logLik = c(
    fit_ref$logLik,
    fit_g_diag$logLik,
    fit_r_diag$logLik,
    fit_cs$logLik,
    fit_fa1_no_geno$logLik,
    fit_fa1$logLik,
    fit_fa2$logLik,
    fit_fa_gm$logLik,
    fit_spatial_r$logLik,
    fit_us$logLik
  ),
  AIC = c(
    fit_ref$AIC,
    fit_g_diag$AIC,
    fit_r_diag$AIC,
    fit_cs$AIC,
    fit_fa1_no_geno$AIC,
    fit_fa1$AIC,
    fit_fa2$AIC,
    fit_fa_gm$AIC,
    fit_spatial_r$AIC,
    fit_us$AIC
  ),
  BIC = c(
    fit_ref$BIC,
    fit_g_diag$BIC,
    fit_r_diag$BIC,
    fit_cs$BIC,
    fit_fa1_no_geno$BIC,
    fit_fa1$BIC,
    fit_fa2$BIC,
    fit_fa_gm$BIC,
    fit_spatial_r$BIC,
    fit_us$BIC
  )
)

model_compare
#>                    model    logLik       AIC       BIC
#> 1              reference -272.8662  553.7325  566.8245
#> 2       G-side diag(env) -271.8293  559.6586  585.8426
#> 3       R-side diag(env) -272.3860  560.7720  586.9560
#> 4      compound symmetry -271.5495  553.0990  569.4640
#> 5    FA k=1 without geno -268.8503  563.7006  606.2496
#> 6                 FA k=1 -275.9114  579.8228  625.6448
#> 7                 FA k=2 -278.6559  593.3118  652.2258
#> 8                FA + GM -269.0085  566.0171  611.8391
#> 9  R-side AR1xAR1 by env -825.2592 1674.5184 1713.7944
#> 10                US GxE -268.8499  543.6998  553.5188

A practical workflow is:

  1. Fit a simple reference model.
  2. Add heterogeneity only when it answers a real design or biological question.
  3. Use compound symmetry when a simple common-correlation GxE model is enough.
  4. Use FA when environment covariance is scientifically important.
  5. Use us() only for modest numbers of environments.
  6. Check predictions and convergence, not only AIC or BIC.

Selected Literature

Advanced variance structures are used when the ordinary independent random effect is too restrictive. Diagonal, compound-symmetry, unstructured, and factor-analytic covariance models provide different levels of flexibility for genotype-by-environment interaction. Genomic and spatial structures extend the same mixed-model framework by replacing simple identity covariance assumptions with relationship or field-layout covariance information.