Field experiments are never perfectly uniform. Soil depth, water, previous crops, or a slope can make neighbouring plots more alike than distant ones, whatever genotypes are planted in them. If the statistical model ignores this, the spatial “noise” leaks into the genotype effects: good genotypes that landed in a poor corner look worse than they are, and vice versa. Spatial models describe this neighbourhood similarity explicitly, so that genotype effects are estimated more accurately.
This vignette introduces spatial modelling with mmes() for readers who already know what a linear mixed model is, but are new to spatial analysis. It covers:
dsumm().For a trial with \(n\) plots, a typical genetic mixed model is
$$ y = X\beta + Zg + e, \qquad g \sim N(0, \sigma^2_g I), \qquad e \sim N(0, R), $$
so the phenotypes have covariance
$$ \operatorname{Var}(y) = V = \sigma^2_g ZZ^\top + R . $$
The textbook choice \(R = \sigma^2_e I\) assumes every plot is independent of every other plot. A spatial model relaxes this assumption in one of two places, which in mmes() correspond to its two covariance arguments:
random (the G side): add random effects that are shared by neighbouring plots, for example row and column effects, a smooth 2D surface (splines), or a correlated “spatial field” with one effect per plot position.rcov (the R side): let the residuals themselves be correlated, so \(R\) is no longer diagonal.Both are valid and often combined. A classic and very successful model (Gilmour et al., 1997) uses correlated residuals plus independent random row and column effects.
Think of the plots along one direction of the field, say the ranges \(1, 2, \ldots, r\). A first-order autoregressive (AR1) correlation says that two plots \(d\) steps apart have correlation \(\rho^{d}\):
$$ K(\rho)_{ij} = \rho^{|i-j|}, \qquad -1 < \rho < 1 . $$
Neighbours (\(d=1\)) have correlation \(\rho\), plots two apart have \(\rho^2\), and so on, so the correlation fades with distance.
ar1 <- function(n, rho) rho^abs(outer(seq_len(n), seq_len(n), "-"))
round(ar1(5, 0.6), 3)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000 0.600 0.36 0.216 0.130
## [2,] 0.600 1.000 0.60 0.360 0.216
## [3,] 0.360 0.600 1.00 0.600 0.360
## [4,] 0.216 0.360 0.60 1.000 0.600
## [5,] 0.130 0.216 0.36 0.600 1.000
A field has two directions. The standard separable model assumes the correlation between two plots is the product of a range correlation and a row correlation:
$$ \operatorname{Cor}(\text{plot}{a}, \text{plot}{b}) = \rho_{\text{range}}^{|\Delta \text{range}|} , \rho_{\text{row}}^{|\Delta \text{row}|}, \qquad K = K(\rho_{\text{range}}) \otimes K(\rho_{\text{row}}), $$
where \(\otimes\) is the Kronecker product. This is exactly what vsm() builds: it multiplies its covariance factors with a Kronecker product and adds one overall variance \(\sigma^2\). The final argument of vsm() names the effect the covariance is attached to:
rcov = ~ vsm(ar1m(range), ar1m(row), ism(units)) gives residuals with covariance \(\sigma^2 K(\rho_\text{range})\otimes K(\rho_\text{row})\); units means “one residual per plot”.random = ~ vsm(ar1m(range), ar1m(row), ism(site)), with site a factor with a single level, gives a random spatial field: one effect for every range-by-row position of the field, correlated in the same way. The residual rcov = ~units then plays the role of an independent plot-level “nugget” (measurement error, plot-to-plot noise).Two practical notes before we start:
ar1m() needs factors whose levels are in field order (1, 2, 3, …), because lags are taken from the level order. ar1m() warns if numeric-looking levels are out of numeric order (for example "1", "10", "2").To see what each model does, we simulate a field of 16 ranges by 24 rows (384 plots) with 96 genotypes in four replicates. The yield of a plot is
$$ \text{yield} = 50 + g_{\text{genotype}} + \xi_{\text{range,row}} + \varepsilon , $$
with genetic values \(g \sim N(0, 9)\), a separable AR1 spatial field \(\xi\) with variance 8, \(\rho_\text{range}=0.7\) and \(\rho_\text{row}=0.4\), and an independent nugget \(\varepsilon \sim N(0,1)\). Six plots are then deleted to mimic missing data.
set.seed(3)
nRange <- 6; nRow <- 8; nGeno <- 12
g <- setNames(rnorm(nGeno, sd=3), sprintf("G%02d", seq_len(nGeno)))
field <- expand.grid(range=seq_len(nRange), row=seq_len(nRow))
field$geno <- sample(rep(names(g), each=4))
trueField <- as.vector(t(chol(8 * kronecker(ar1(nRow, 0.4), ar1(nRange, 0.7)))) %*%
rnorm(nrow(field)))
field$yield <- 50 + g[field$geno] + trueField + rnorm(nrow(field), sd=1)
field <- field[-sample(nrow(field), 6), ]
field$geno <- factor(field$geno)
field$rangef <- factor(field$range, levels=seq_len(nRange))
field$rowf <- factor(field$row, levels=seq_len(nRow))
field$site <- factor("F1")
head(field)
## range row geno yield rangef rowf site
## 1 1 1 G10 50.56466 1 1 F1
## 2 2 1 G03 49.25608 2 1 F1
## 3 3 1 G12 42.47582 3 1 F1
## 4 4 1 G05 46.23970 4 1 F1
## 6 6 1 G07 48.53419 6 1 F1
## 8 2 2 G10 54.09992 2 2 F1
All models share the same fixed effects (an intercept) and a random genotype effect; they differ only in how the field is described. Each fit is capped at one REML iteration to keep this chunk quick; allow convergence for a complete analysis.
fits <- list(
# M0: no spatial model
iid = mmes(yield ~ 1,
random=~geno,
data=field,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE))),
# M1: independent random range and row effects (large-scale strips)
rowcol = mmes(yield ~ 1,
random=~geno + rangef + rowf,
data=field,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE))),
# M2: a smooth 2D surface built from tensor-product P-splines
spline = mmes(yield ~ 1,
random=~geno + vsm(ism(spl2Dc(row, range)$Z$`A:all`)),
data=field,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE))),
# M3: correlated residuals (R side)
ar1Res = mmes(yield ~ 1,
random=~geno,
rcov=~vsm(ar1m(rangef), ar1m(rowf), ism(units)),
data=field,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE))),
# M4: a random AR1 x AR1 spatial field (G side) plus an iid nugget
ar1Fld = mmes(yield ~ 1,
random=~geno + vsm(ar1m(rangef), ar1m(rowf), ism(site)),
rcov=~units,
data=field,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
)
Because all models have the same fixed effects, their REML log-likelihoods (and AIC) can be compared directly; higher log-likelihood and lower AIC are better. Since we simulated the data, we can also compute how well each model ranks the genotypes: the correlation between the estimated genotype effects (BLUPs) and the true genetic values.
compare <- t(sapply(fits, function(m){
blup <- m$uList[[1]][, 1]
c(logLik = tail(as.numeric(m$llik), 1),
AIC = m$AIC,
genVar = covparams_mmes(m)$estimate[1],
accuracy = cor(blup[names(g)], g))
}))
knitr::kable(round(compare, 3))
| logLik | AIC | genVar | accuracy | |
|---|---|---|---|---|
| iid | -17.580 | 37.161 | 6.487 | 0.791 |
| rowcol | -17.095 | 36.190 | 6.487 | 0.846 |
| spline | -17.164 | 36.327 | 6.487 | 0.797 |
| ar1Res | -12.616 | 27.232 | 6.487 | 0.855 |
| ar1Fld | -17.959 | 37.919 | 6.487 | 0.836 |
The pattern is typical of real trials. Ignoring the field (M0) gives the worst fit and the least accurate ranking. Row and column effects (M1) capture strips; splines (M2) capture a smooth trend; the AR1 models (M3, M4) capture the local, patchy variation and fit best. Better spatial modelling translates directly into more accurate genotype rankings, which is the reason to do it.
covparams_mmes() reports every covariance parameter on its natural scale.
knitr::kable(covparams_mmes(fits$ar1Fld)[, c("factor", "parameter", "estimate")], digits=3)
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 6.487 |
| ar1m(rangef) | variance | 6.587 |
| ar1m(rangef) | rho | 0.300 |
| ar1m(rowf) | rho | 0.300 |
| sigma2 | sigma2 | 6.587 |
Compare with the truth: genetic variance 9, spatial variance 8, \(\rho_\text{range}=0.7\), \(\rho_\text{row}=0.4\), and nugget 1. The realized genetic variance can differ from 9 because only 12 genotypes are sampled. Note that in a Kronecker product the single variance is reported once, on the first factor (ar1m(rangef)); the second factor only contributes its correlation.
The difference between M3 and M4 is the nugget. M3 assumes that plot-level noise is perfectly correlated with the spatial pattern, while M4 separates a smooth-ish correlated field from independent plot error. With real data both are worth trying; the likelihood tells you which describes your field better.
The single-kernel spline used in M2 (spl2Dc()) fits the whole 2D surface with one variance component. The SpATS package (Rodriguez-Alvarez et al., 2018) instead splits the surface into five components (smooth trends along each direction and their interactions), each with its own variance. spl2Dmats() builds the same five-component design-matrix structure, so mmes() can reproduce the SpATS fit when configured with matching segmentation. The classic Yates oats trial is used here, with SpATS output shown for reference. To keep this example quick, the active fit uses fewer spline segments and one REML iteration; use the fuller segmentation and allow convergence for a complete analysis.
data(DT_yatesoats, package="enhancer")
DT <- DT_yatesoats
DT$row <- as.numeric(as.character(DT$row))
DT$col <- as.numeric(as.character(DT$col))
DT$R <- as.factor(DT$row)
DT$C <- as.factor(DT$col)
# SPATS MODEL
# m1.SpATS <- SpATS(response = "Y",
# spatial = ~ PSANOVA(col, row, nseg = c(14,21), degree = 3, pord = 2),
# genotype = "V", fixed = ~ 1,
# random = ~ R + C, data = DT,
# control = list(tolerance = 1e-04))
#
# summary(m1.SpATS, which = "variances")
#
# Spatial analysis of trials with splines
#
# Response: Y
# Genotypes (as fixed): V
# Spatial: ~PSANOVA(col, row, nseg = c(14, 21), degree = 3, pord = 2)
# Fixed: ~1
# Random: ~R + C
#
#
# Number of observations: 72
# Number of missing data: 0
# Effective dimension: 17.09
# Deviance: 483.405
#
# Variance components:
# Variance SD log10(lambda)
# R 1.277e+02 1.130e+01 0.49450
# C 2.673e-05 5.170e-03 7.17366
# f(col) 4.018e-15 6.339e-08 16.99668
# f(row) 2.291e-10 1.514e-05 12.24059
# f(col):row 1.025e-04 1.012e-02 6.59013
# col:f(row) 8.789e+01 9.375e+00 0.65674
# f(col):f(row) 8.036e-04 2.835e-02 5.69565
#
# Residual 3.987e+02 1.997e+01
# SOMMER MODEL
M <- spl2Dmats(x.coord.name = "col", y.coord.name = "row", data=DT,
nseg =c(6,8), degree = c(3,3), penaltyord = c(2,2)
)
mix <- mmes(Y~V,
random=~ R + C + vsm(ism(M$fC)) + vsm(ism(M$fR)) +
vsm(ism(M$fC.R)) + vsm(ism(M$C.fR)) +
vsm(ism(M$fC.fR)),
rcov=~units,
data=M$data,
control = mmes_control(formulation = "henderson",
convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
summary(mix)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(ism(R)) sigma2 45.32431 44.01858 1.029663063
## 2 vsm(ism(C)) sigma2 171.82077 182.11168 0.943491207
## 3 vsm(ism(M$fC)) sigma2 72.75302 11920.47763 0.006103197
## 4 vsm(ism(M$fR)) sigma2 72.75302 464.65518 0.156574215
## 5 vsm(ism(M$fC.R)) sigma2 72.75302 472.04613 0.154122691
## 6 vsm(ism(M$C.fR)) sigma2 72.75302 467.06371 0.155766801
## 7 vsm(ism(M$fC.fR)) sigma2 72.75302 519.85046 0.139949900
## 8 vsm(ism(units)) sigma2 363.76510 72.48389 5.018564665
The five spline variances together describe the surface; when several of them shrink to zero, the simpler single-kernel spl2Dc() model of M2 is usually adequate.
Rotating spline penalties into independent random effects can make the design
matrices dense, even when the original B-splines are local. In Henderson fits,
building diagonal residual score cross-products by expanding every pair of
nonzeros in every observation row can then cost more than factorization.
score = list(denseCrossproduct = "auto") recognizes expensive dense-row
cases and uses an exact weighted sparse Gram product instead. It changes
neither the model nor the coefficient solver, selected-inverse traces or AI
optimizer. Existing elimination policies are unchanged.
For Henderson fitting, solve = list(method = "cholmod") selects the supernodal
coefficient solver. Use diagnostics = list(profile = "force") to identify
whether factorization or score construction dominates. Diagonal residual score
cross-products are accumulated by sparse matrix multiplication, without
expanding every pair of nonzeros in every design row into duplicate triplets
when the dense-design path is selected. "off" restores the original path;
"force" bypasses cost gates for eligible diagonal residual scores. Auto
requires at least 64 nonzeros in an active row, a million row-product
contributions, eight times the squared coefficient dimension, and sufficient
workspace. Small/sparse incidence models, non-diagonal effective residual
precision, PCG, PQL and other equation representations retain their existing
paths. A forced request cannot bypass structural or memory checks.
control <- mmes_control(
formulation = "henderson",
solve = list(method = "cholmod"),
score = list(denseCrossproduct = "auto"),
diagnostics = list(profile = "force"))
# Inspect fit$engineDiagnostics$denseScoreCrossproducts and the MME banner.
# Set denseCrossproduct = "off" to benchmark the unchanged legacy route.
When comparing with SpATS, also match the unpenalized spatial terms, genotype
treatment, basis nesting and convergence criteria. For a second-order penalty,
the spline null space includes linear column and row trends and their bilinear
interaction; spl2Dmats() returns these as TP.col, TP.row and TP.CxR in
M$data, but they are not added automatically to a formula such as Y ~ V.
Absolute likelihoods may use different constants or response scaling, so
compare fitted values and variance estimates before comparing reported criteria.
Breeding programmes rarely analyse one field. When trials are analysed together, we want to share information about genotypes across trials, while each field keeps its own spatial pattern. The key modelling question is: which spatial parameters are shared between trials, and which are trial-specific?
We simulate three trials of different sizes, testing the same 80 genotypes, each with its own residual variance and its own range and row correlations. To keep this vignette quick, active fits below use one REML iteration. Treat estimates and likelihoods as illustrations, not converged results; allow convergence before drawing conclusions:
| trial | size (range x row) | variance | \(\rho_\text{range}\) | \(\rho_\text{row}\) |
|---|---|---|---|---|
| T1 | 16 x 20 | 9 | 0.8 | 0.2 |
| T2 | 12 x 24 | 16 | 0.1 | 0.7 |
| T3 | 14 x 18 | 4 | 0.5 | 0.5 |
simTrial <- function(trial, nRange, nRow, rhoRange, rhoRow, sigma2, g){
plots <- expand.grid(range=seq_len(nRange), row=seq_len(nRow))
plots$geno <- sample(rep(names(g), length.out=nrow(plots)))
K <- sigma2 * kronecker(ar1(nRow, rhoRow), ar1(nRange, rhoRange))
plots$yield <- 50 + g[plots$geno] + as.vector(t(chol(K)) %*% rnorm(nrow(plots)))
plots$trial <- trial
plots[-sample(nrow(plots), round(0.05 * nrow(plots))), ]
}
set.seed(2026)
gMET <- setNames(rnorm(80, sd=3), sprintf("G%02d", 1:80))
MET <- rbind(simTrial("T1", 16, 20, 0.8, 0.2, 9, gMET),
simTrial("T2", 12, 24, 0.1, 0.7, 16, gMET),
simTrial("T3", 14, 18, 0.5, 0.5, 4, gMET))
MET$trial <- factor(MET$trial)
MET$geno <- factor(MET$geno)
MET$range <- factor(MET$range, levels=1:16)
MET$row <- factor(MET$row, levels=1:24)
table(MET$trial)
##
## T1 T2 T3
## 304 274 239
Note that the range and row factors use the same levels in all trials. Each trial only uses the levels it needs; that is fine, because trials are modelled as independent fields.
vsm()Adding a diagonal trial factor dsm(trial) to the residual vsm() gives each trial its own residual variance, and makes plots in different trials independent (the off-diagonal elements of dsm() are zero):
$$ R = \sigma^2 , D_\text{trial} \otimes K(\rho_\text{range}) \otimes K(\rho_\text{row}) . $$
Because this is a Kronecker product, each factor appears once: there is a single \(\rho_\text{range}\) and a single \(\rho_\text{row}\), shared by all trials.
mShared <- mmes(yield ~ trial,
random=~geno,
rcov=~vsm(dsm(trial), ar1m(range), ar1m(row), ism(units)),
data=MET,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
knitr::kable(covparams_mmes(mShared)[, c("factor", "parameter", "estimate")], digits=3)
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 7.961 |
| dsm(trial) | variance[T1] | 7.693 |
| dsm(trial) | variance[T2] | 13.038 |
| dsm(trial) | variance[T3] | 5.666 |
| ar1m(range) | rho | 0.300 |
| ar1m(row) | rho | 0.300 |
The trial variances are sensible, but the correlations are a compromise between very different fields (the truth for \(\rho_\text{range}\) ranges from 0.1 to 0.8).
dsumm()To give every trial its own variance and its own correlations, the residual must be a direct sum, a block-diagonal matrix with one independent block per trial:
$$ R = \bigoplus_{t} \sigma^2_t, K(\rho_{\text{range},t}) \otimes K(\rho_{\text{row},t}) = \begin{pmatrix} R_{T1} & 0 & 0 \ 0 & R_{T2} & 0 \ 0 & 0 & R_{T3} \end{pmatrix}. $$
dsumm() takes a residual vsm() term describing one trial and replicates it, with all its parameters, for every level of by. The trial-specific variances come from the scale of each copy, so dsm(trial) is not needed.
mDsum <- mmes(yield ~ trial,
random=~geno,
rcov=~dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by=trial),
data=MET,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
knitr::kable(covparams_mmes(mDsum)[, c("factor", "section", "parameter", "estimate")],
digits=3)
| factor | section | parameter | estimate |
|---|---|---|---|
| sigma2 | NA | sigma2 | 7.961 |
| trial | T1 | variance | 8.799 |
| trial | T2 | variance | 8.799 |
| trial | T3 | variance | 8.799 |
| ar1m(range) | T1 | rho | 0.300 |
| ar1m(range) | T2 | rho | 0.300 |
| ar1m(range) | T3 | rho | 0.300 |
| ar1m(row) | T1 | rho | 0.300 |
| ar1m(row) | T2 | rho | 0.300 |
| ar1m(row) | T3 | rho | 0.300 |
The section column tells which trial each parameter belongs to. At convergence, estimates should track the truth in each trial. The first row (sigma2) is the internal overall scale; the trial variances are the variance rows.
The shared model is nested in the dsumm() model: it is the special case in which all trials have the same correlations. Once both models converge, they can be compared with a likelihood-ratio test, with degrees of freedom equal to the number of extra parameters, here \(2 \times (3 - 1) = 4\) correlations. Both models have the same fixed and random effects, so their REML likelihoods are comparable.
logLikShared <- tail(as.numeric(mShared$llik), 1)
logLikDsum <- tail(as.numeric(mDsum$llik), 1)
LR <- 2 * (logLikDsum - logLikShared)
c(LR = LR, df = 4, p.value = pchisq(LR, df=4, lower.tail=FALSE))
## LR df p.value
## -89.66733 4.00000 1.00000
As a rule of thumb:
dsumm() when trials differ in size, orientation, or environment, which is common in multi-location programmes. With \(T\) trials and two AR1 factors, dsumm() estimates \(3T\) residual parameters instead of \(T + 2\).(Test statistics for variance parameters on the boundary of the parameter space, such as variances near zero, need more care; correlations inside \((-1, 1)\) are not affected.)
levels=Sometimes only some trials have row and range coordinates, or some trials are too small for a spatial model. The levels argument restricts the spatial structure to the listed trials; every other trial gets an independent residual with its own variance (\(\sigma^2_t I\)). Those trials are kept in the analysis even when their coordinates are missing.
MET2 <- MET
MET2$range[MET2$trial == "T3"] <- NA # T3 has no field coordinates
MET2$row[MET2$trial == "T3"] <- NA
mLevels <- mmes(yield ~ trial,
random=~geno,
rcov=~dsumm(vsm(ar1m(range), ar1m(row), ism(units)),
by=trial, levels=c("T1", "T2")),
data=MET2,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
knitr::kable(covparams_mmes(mLevels)[, c("factor", "section", "parameter", "estimate")],
digits=3)
| factor | section | parameter | estimate |
|---|---|---|---|
| sigma2 | NA | sigma2 | 7.961 |
| trial | T1 | variance | 8.799 |
| trial | T2 | variance | 8.799 |
| trial | T3 | variance | 8.799 |
| ar1m(range) | T1 | rho | 0.300 |
| ar1m(range) | T2 | rho | 0.300 |
| ar1m(row) | T1 | rho | 0.300 |
| ar1m(row) | T2 | rho | 0.300 |
c(used = nrow(mLevels$y), available = nrow(MET2))
## used available
## 817 817
The random-effect route from section 3 also extends to several trials. A diagonal trial factor in random gives each trial its own spatial field variance while sharing the correlations, for example
random = ~ geno + vsm(dsm(trial), ar1m(range), ar1m(row), ism(site))
and trial-specific spline surfaces are obtained with vsm(dsm(trial), ism(spl2Dc(row, range)$Z$`A:all`)). Residual structures and random spatial effects can be combined freely, for example dsumm() residuals plus random row and column effects within trials (vsm(dsm(trial), ism(rowf))).
An AR1 x AR1 residual has a dense covariance matrix, but its inverse is sparse when the field is a complete grid: every plot is linked only to its eight neighbours. mmes() uses this in two ways.
NA responses), the residual precision of the observed plots is dense again. Like ASReml, mmes() then adds the missing cells as pseudo-observations, with one dummy fixed effect each. This completes the grid without changing the REML likelihood, the variance components, the BLUEs or the BLUPs. The dummy effects are removed from the returned object, so b, u, fitted(), predict(), AIC and BIC refer to the original data. A cell is only added if both its range and its row levels are observed in that trial.Both are on by default (representation = list(gridAugment = "auto")). They apply to REML fits of Gaussian models with the ldlt or cholmod solvers. The engine diagnostics show whether augmentation was used:
str(mDsum$engineDiagnostics$gridAugmentation)
## List of 4
## $ active : logi TRUE
## $ cells : int 43
## $ observed: int 817
## $ blocks : int 3
mDsumOff <- mmes(yield ~ trial,
random=~geno,
rcov=~dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by=trial),
data=MET,
control = mmes_control(diagnostics = list(verbose = FALSE),
representation = list(gridAugment = "off"),
convergence = list(maxIters = 1L)))
c(augmented = tail(as.numeric(mDsum$llik), 1), original = tail(as.numeric(mDsumOff$llik), 1))
## augmented original
## -119.1119 -119.1119
The savings grow with the size of the trials and the number of missing plots. Use gridAugment = "force" to get an error that explains why augmentation is not possible, instead of a silent fallback. The residual must include ism(units), as in all examples above. vsm(ar1m(range), ar1m(row)) without ism(units) is a different, row-blocked structure, and it is not treated as a field grid.
These improvements also apply with fam(locale, k) or rrm(locale, k)
random effects, not only rank one. For augmented grids, the cholmod solver
prepares sparse selected-inverse entries before considering dense random
inverse-block caches. This is particularly beneficial for FA models with a
dense genomic relationship matrix. Memory limits and exact fallback solves
still apply; this changes computation, not the likelihood or covariance
model. The profiling counter
engineDiagnostics$profile$counters$patternInverseCacheBuilds records this
path when diagnostics = list(profile = "force"). Increasing RR rank still
adds latent coefficients and random-effect trace work, so higher ranks need
not have the same runtime or speedup.
Another way to fit field trends is to move the AR1 x AR1 process from the residual to a random plot-level effect, and keep an independent residual (the “nugget”) for each trial:
$$ y = X\beta + Z_g g + \sum_t Z_t s_t + e,\qquad s_t \sim N\big(0,\ \sigma^2_{s,t},K(\rho_{\text{range},t}) \otimes K(\rho_{\text{row},t})\big),\qquad e \sim N(0, \textstyle\bigoplus_t \sigma^2_{e,t} I). $$
With atm(), each trial gets its own spatial field, variance and correlations. A single-level factor site gives one field per trial.
MET$site <- factor("field")
mField <- mmes(yield ~ trial,
random=~geno +
vsm(atm(trial, "T1"), ar1m(range), ar1m(row), ism(site)) +
vsm(atm(trial, "T2"), ar1m(range), ar1m(row), ism(site)) +
vsm(atm(trial, "T3"), ar1m(range), ar1m(row), ism(site)),
rcov=~vsm(dsm(trial), ism(units)),
data=MET,
control = mmes_control(convergence = list(maxIters = 1L),
diagnostics = list(verbose = FALSE)))
knitr::kable(covparams_mmes(mField)[, c("factor", "parameter", "estimate")], digits=3)
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 7.961 |
| atm(trial, “T1”) | variance[T1] | 2.933 |
| ar1m(range) | rho | 0.300 |
| ar1m(row) | rho | 0.300 |
| atm(trial, “T2”) | variance[T2] | 2.933 |
| ar1m(range) | rho | 0.300 |
| ar1m(row) | rho | 0.300 |
| atm(trial, “T3”) | variance[T3] | 2.933 |
| ar1m(range) | rho | 0.300 |
| ar1m(row) | rho | 0.300 |
| dsm(trial) | variance[T1] | 7.693 |
| dsm(trial) | variance[T2] | 13.038 |
| dsm(trial) | variance[T3] | 5.666 |
cor(mDsum$u[grep("^G", rownames(mDsum$u))], mField$u[grep("^G", rownames(mField$u))])
## [1] 0.9936437
At convergence, spatial variances and correlations should be close to those of the dsumm() residual model. Residual variances should be near zero because these data were simulated without a nugget, and genotype BLUPs should be similar.
When to use which formulation:
dsumm() residual with automatic augmentation is usually considerably faster, because it adds only the missing cells.residuals(fit) against range and row) show whether spatial patterns remain.ar1m().control$formulation = "direct" uses observation-space inversion, but needs memory proportional to the square of the number of plots.| model | ASReml-R | sommer |
|---|---|---|
| single trial, AR1 x AR1 residual | residual = ~ar1(range):ar1(row) |
rcov = ~vsm(ar1m(range), ar1m(row), ism(units)) |
| several trials, shared correlations | residual = ~diag(trial):ar1(range):ar1(row) |
rcov = ~vsm(dsm(trial), ar1m(range), ar1m(row), ism(units)) |
| several trials, trial-specific | residual = ~dsum(~ar1(range):ar1(row)| trial) |
rcov = ~dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by=trial) |
| only some trials spatial | dsum(~ar1(range):ar1(row)| trial, levels=lv) |
dsumm(..., by=trial, levels=lv) |
| random row / column effects | random = ~at(trial):row + at(trial):col |
random = ~vsm(dsm(trial), ism(row)) + vsm(dsm(trial), ism(col)) |
In ASReml-R, dsum(... | trial, levels=lv) covers only the listed sections, and the remaining trials need their own residual term; dsumm() gives them an independent residual automatically.
Covarrubias-Pazaran G. 2016. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6):1-15.
Covarrubias-Pazaran G. 2018. Software update: Moving the R package sommer to multivariate mixed models for genome-assisted prediction. doi: https://doi.org/10.1101/354639
Cullis B.R., Gleeson A.C. 1991. Spatial analysis of field experiments: an extension to two dimensions. Biometrics 47(4):1449-1460.
Gilmour A.R., Cullis B.R., Verbyla A.P. 1997. Accounting for natural and extraneous variation in the analysis of field experiments. Journal of Agricultural, Biological, and Environmental Statistics 2:269-293.
Gilmour et al. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51(4):1440-1450.
Henderson C.R. 1975. Best Linear Unbiased Estimation and Prediction under a Selection Model. Biometrics vol. 31(2):423-447.
Lee, D.-J., Durban, M., and Eilers, P.H.C. (2013). Efficient two-dimensional smoothing with P-spline ANOVA mixed models and nested bases. Computational Statistics and Data Analysis, 61, 22 - 37.
Rodriguez-Alvarez, Maria Xose, et al. Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spatial Statistics 23 (2018): 52-71.