Fit Models for Hierarchical Agricultural Experiments
agriDesignR contributors
Source:vignettes/mixed-models.Rmd
mixed-models.Rmd
library(agriDesignR)
data("wheat_splitplot", package = "agriDesignR")Audit topology before fitting
The included wheat data represent four blocks, two whole-plot
temperature levels, and five genotype sub-plots. Tell
check_design() how units were randomized so the audit can
distinguish an RCBD factorial from a split plot.
design <- check_design(
data = wheat_splitplot,
response = "grain_yield",
main_plot = "temperature",
sub_plot = "genotype",
block = "block"
)
design$design_type
#> [1] "Split-Plot Design (Parcelas Divididas)"
design$balance_status
#> [1] "balanced"suggest_model() turns topology into inspectable formulas
and an engine recommendation. Treat recommendation as a transparent
starting point: the estimand, randomization, and scientific question
still determine final model.
suggestion <- suggest_model(design)
suggestion$recommended_engine
#> [1] "lme4"
suggestion$fixed_formula_str
#> [1] "grain_yield ~ temperature * genotype"
suggestion$random_formula_lme4
#> [1] "(1 | block) + (1 | block:temperature)"For an RCBD with only a few blocks, the default
block_as = "auto" keeps the historical fixed-block
analysis. If blocks are sampled from a wider population, request the
mixed specification explicitly; this keeps the block decision
independent from the computational engine:
fit_fixed <- fit_experiment(
data = wheat_splitplot,
response = "grain_yield",
main_plot = "temperature",
sub_plot = "genotype",
block = "block",
block_as = "fixed",
engine = "lm"
)
fit_random <- fit_experiment(
data = wheat_splitplot,
response = "grain_yield",
main_plot = "temperature",
sub_plot = "genotype",
block = "block",
block_as = "random",
engine = "lme4"
)
#> Warning: Random block effect requested with 4 levels; variance estimates may be
#> unstable or singular.
#> Registered S3 method overwritten by 'lme4':
#> method from
#> na.action.merMod car
fit_random$fixed_formula
#> grain_yield ~ temperature * genotype
#> <environment: 0x555823775ea8>
fit_random$random_formula
#> [1] "(1 | block) + (1 | block:temperature)"With four or fewer random-effect levels, inspect the singularity diagnostic and variance estimate before treating the random-block model as reliable.
Fit a mixed model
fit_experiment() fits the suggested model and stores the
data, call, ANOVA table, and assumption audit in an
agri_fitted_model object. Default
engine = "auto" uses design information; choose an engine
explicitly when a preregistered analysis requires one.
fit <- fit_experiment(
data = wheat_splitplot,
response = "grain_yield",
main_plot = "temperature",
sub_plot = "genotype",
block = "block"
)
class(fit)
#> [1] "agri_fitted_model" "list"
fit$engine
#> [1] "lme4"
fit$anova_table
#> Type III Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> temperature 359.95 359.95 1 3 75.346 0.003217 **
#> genotype 1584.58 396.15 4 24 82.923 1.086e-13 ***
#> temperature:genotype 238.75 59.69 4 24 12.494 1.238e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1S3 methods provide concise summaries while keeping fitted model available for diagnostics and downstream inference:
summary(fit)
#>
#> ======================================================================
#> agriDesignR: Fitted Experimental Model Engine
#> ======================================================================
#> Design Layout : Split-Plot Design (Parcelas Divididas)
#> Model Class : LMM (Linear Mixed-Effects Model)
#> Engine / Family : lme4 (gaussian)
#> Fixed Formula : grain_yield ~ temperature * genotype
#> Random Terms : (1 | block) + (1 | block:temperature)
#>
#> Health Check : [!] MILD DEVIATIONS DETECTED (Robust under ANOVA)
#> - Normality : [PASSED] Residuals follow a normal distribution (Shapiro-Wilk W = 0.98, p = 0.6964).
#> - Homoscedasticity: [ACCEPTABLE] Moderate variance spread across groups (Variance ratio = 21.36x, Levene p = 0.18). ANOVA F-tests remain robust when group sizes are balanced.
#>
#> Type III ANOVA / Deviance Significance Table
#> ----------------------------------------------------------------------
#> Type III Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> temperature 359.95 359.95 1 3 75.346 0.003217 **
#> genotype 1584.58 396.15 4 24 82.923 1.086e-13 ***
#> temperature:genotype 238.75 59.69 4 24 12.494 1.238e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ======================================================================
#> Tip: Use 'summary(fit)' for parameter estimates, 'plot(fit)' for 4-panel residual graphs,
#> or 'remedy(fit)' if assumption violations require transformations or varIdent.
#>
#>
#> Underlying Model Parameter Estimates
#> ----------------------------------------------------------------------
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: full_formula
#> Data: raw_data
#>
#> REML criterion at convergence: 162.1
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -1.3829 -0.5391 -0.1395 0.5613 1.7512
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> block:temperature (Intercept) 6.750 2.598
#> block (Intercept) 9.148 3.025
#> Residual 4.777 2.186
#> Number of obs: 40, groups: block:temperature, 8; block, 4
#>
#> Fixed effects:
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 57.763 2.273 6.906 25.407 4.44e-08 ***
#> temperatureHeat_36C -18.133 2.401 6.514 -7.553 0.000188 ***
#> genotypeG2 4.885 1.546 24.000 3.161 0.004223 **
#> genotypeG3 -3.765 1.546 24.000 -2.436 0.022643 *
#> genotypeG4 7.947 1.546 24.000 5.142 2.90e-05 ***
#> genotypeG5 1.117 1.546 24.000 0.723 0.476632
#> temperatureHeat_36C:genotypeG2 -0.185 2.186 24.000 -0.085 0.933249
#> temperatureHeat_36C:genotypeG3 -5.735 2.186 24.000 -2.624 0.014878 *
#> temperatureHeat_36C:genotypeG4 9.408 2.186 24.000 4.304 0.000244 ***
#> temperatureHeat_36C:genotypeG5 1.985 2.186 24.000 0.908 0.372814
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Correlation of Fixed Effects:
#> (Intr) tmH_36C gntyG2 gntyG3 gntyG4 gntyG5 tH_36C:G2 tH_36C:G3
#> tmprtrH_36C -0.528
#> genotypeG2 -0.340 0.322
#> genotypeG3 -0.340 0.322 0.500
#> genotypeG4 -0.340 0.322 0.500 0.500
#> genotypeG5 -0.340 0.322 0.500 0.500 0.500
#> tmpH_36C:G2 0.240 -0.455 -0.707 -0.354 -0.354 -0.354
#> tmpH_36C:G3 0.240 -0.455 -0.354 -0.707 -0.354 -0.354 0.500
#> tmpH_36C:G4 0.240 -0.455 -0.354 -0.354 -0.707 -0.354 0.500 0.500
#> tmpH_36C:G5 0.240 -0.455 -0.354 -0.354 -0.354 -0.707 0.500 0.500
#> tH_36C:G4
#> tmprtrH_36C
#> genotypeG2
#> genotypeG3
#> genotypeG4
#> genotypeG5
#> tmpH_36C:G2
#> tmpH_36C:G3
#> tmpH_36C:G4
#> tmpH_36C:G5 0.500
anova(fit)
#> Type III Analysis of Variance Table with Kenward-Roger's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> temperature 359.95 359.95 1 3 75.346 0.003217 **
#> genotype 1584.58 396.15 4 24 82.923 1.086e-13 ***
#> temperature:genotype 238.75 59.69 4 24 12.494 1.238e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Decompose an interaction with simple effects
When an interaction is scientifically relevant, compare genotypes
within each temperature level instead of interpreting averaged genotype
means alone. agri_posthoc() returns estimated means and
adjusted pairwise contrasts.
post_simple <- agri_posthoc(
fit,
trt = "genotype",
by = "temperature",
method = "tukey"
)
#> Note: adjust = "tukey" was changed to "sidak"
#> because "tukey" is only appropriate for one set of pairwise comparisons
post_simple$means_table
#> temperature = Control_22C:
#> genotype emmean SE df lower.CL upper.CL Group
#> G1 57.7625 2.273493 6.91 52.37172 63.15328 ab
#> G2 62.6475 2.273493 6.91 57.25672 68.03828 cd
#> G3 53.9975 2.273493 6.91 48.60672 59.38828 a
#> G4 65.7100 2.273493 6.91 60.31922 71.10078 c
#> G5 58.8800 2.273493 6.91 53.48922 64.27078 bd
#>
#> temperature = Heat_36C:
#> genotype emmean SE df lower.CL upper.CL Group
#> G1 39.6300 2.273493 6.91 34.23922 45.02078 a
#> G2 44.3300 2.273493 6.91 38.93922 49.72078 b
#> G3 30.1300 2.273493 6.91 24.73922 35.52078 c
#> G4 56.9850 2.273493 6.91 51.59422 62.37578 d
#> G5 42.7325 2.273493 6.91 37.34172 48.12328 ab
#>
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
post_simple$contrasts_table
#> temperature = Control_22C:
#> contrast estimate SE df t.ratio p.value
#> G1 - G2 -4.8850 1.545519 24 -3.161 0.0312
#> G1 - G3 3.7650 1.545519 24 2.436 0.1400
#> G1 - G4 -7.9475 1.545519 24 -5.142 0.0003
#> G1 - G5 -1.1175 1.545519 24 -0.723 0.9491
#> G2 - G3 8.6500 1.545519 24 5.597 <0.0001
#> G2 - G4 -3.0625 1.545519 24 -1.982 0.3047
#> G2 - G5 3.7675 1.545519 24 2.438 0.1395
#> G3 - G4 -11.7125 1.545519 24 -7.578 <0.0001
#> G3 - G5 -4.8825 1.545519 24 -3.159 0.0314
#> G4 - G5 6.8300 1.545519 24 4.419 0.0016
#>
#> temperature = Heat_36C:
#> contrast estimate SE df t.ratio p.value
#> G1 - G2 -4.7000 1.545519 24 -3.041 0.0407
#> G1 - G3 9.5000 1.545519 24 6.147 <0.0001
#> G1 - G4 -17.3550 1.545519 24 -11.229 <0.0001
#> G1 - G5 -3.1025 1.545519 24 -2.007 0.2927
#> G2 - G3 14.2000 1.545519 24 9.188 <0.0001
#> G2 - G4 -12.6550 1.545519 24 -8.188 <0.0001
#> G2 - G5 1.5975 1.545519 24 1.034 0.8373
#> G3 - G4 -26.8550 1.545519 24 -17.376 <0.0001
#> G3 - G5 -12.6025 1.545519 24 -8.154 <0.0001
#> G4 - G5 14.2525 1.545519 24 9.222 <0.0001
#>
#> Degrees-of-freedom method: kenward-roger
#> P value adjustment: tukey method for comparing a family of 5 estimatesThe letter family is defined by trt and by,
the p-value adjustment, and the alpha level.
plot_publication() uses the mapped grouping factor within
each x-axis (and facet) level by default; pass trt and
by to override that visual comparison when the scientific
question requires a different direction.
For a classical analysis of a genuinely balanced design, pass
engine = "lm" explicitly and report why fixed-effects error
structure matches the randomization. For missing plots or hierarchical
randomization, retain the mixed-model structure and report the number of
independent units.