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.
## [,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.000A 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)\). Fifteen plots are then deleted to mimic missing data.
set.seed(3)
nRange <- 16; nRow <- 24; nGeno <- 96
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), 15), ]
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 G29 50.06152 1 1 F1
## 2 2 1 G57 46.12666 2 1 F1
## 3 3 1 G43 55.73973 3 1 F1
## 4 4 1 G02 47.48929 4 1 F1
## 5 5 1 G65 55.61240 5 1 F1
## 6 6 1 G05 47.60453 6 1 F1All models share the same fixed effects (an intercept) and a random genotype effect; they differ only in how the field is described.
fits <- list(
# M0: no spatial model
iid = mmes(yield ~ 1, random=~geno, data=field, verbose=FALSE),
# M1: independent random range and row effects (large-scale strips)
rowcol = mmes(yield ~ 1, random=~geno + rangef + rowf, data=field, 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, verbose=FALSE),
# M3: correlated residuals (R side)
ar1Res = mmes(yield ~ 1, random=~geno,
rcov=~vsm(ar1m(rangef), ar1m(rowf), ism(units)),
data=field, 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, 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 | -141.227 | 284.454 | 6.979 | 0.890 |
| rowcol | -119.592 | 241.185 | 7.469 | 0.908 |
| spline | -108.224 | 218.449 | 6.899 | 0.915 |
| ar1Res | -81.485 | 164.971 | 7.355 | 0.943 |
| ar1Fld | -79.873 | 161.746 | 7.232 | 0.941 |
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.
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 7.232 |
| ar1m(rangef) | variance | 7.134 |
| ar1m(rangef) | rho | 0.666 |
| ar1m(rowf) | rho | 0.468 |
| sigma2 | sigma2 | 1.007 |
Compare with the truth: genetic variance 9 (the 96 values actually
drawn have variance 6.7, which is what the model can recover), spatial
variance 8, \(\rho_\text{range}=0.7\),
\(\rho_\text{row}=0.4\), nugget 1. 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
design matrices, so mmes() reproduces the SpATS fit. The
classic Yates oats trial is used here, with SpATS output shown for
reference.
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(14,21), degree = c(3,3), penaltyord = c(2,2)
)
mix <- mmes(Y~V, henderson = TRUE,
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, verbose=FALSE,
data=M$data)
summary(mix)$varcomp## term parameter estimate StdError Zratio
## 1 vsm(ism(R)) sigma2 107.0392029 88.92513 1.2037002441
## 2 vsm(ism(C)) sigma2 177.1312530 176.53843 1.0033580445
## 3 vsm(ism(M$fC)) sigma2 0.4439506 1448.07385 0.0003065801
## 4 vsm(ism(M$fR)) sigma2 0.4967557 62.69920 0.0079228401
## 5 vsm(ism(M$fC.R)) sigma2 1.4262622 38.30944 0.0372300497
## 6 vsm(ism(M$C.fR)) sigma2 0.4035155 20.04298 0.0201325114
## 7 vsm(ism(M$fC.fR)) sigma2 0.8736894 60.13291 0.0145293046
## 8 vsm(ism(units)) sigma2 490.7116628 105.11925 4.6681426172The 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.
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:
| 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 239Note 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.
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, verbose=FALSE)
knitr::kable(covparams_mmes(mDsum)[, c("factor", "section", "parameter", "estimate")],
digits=3)| factor | section | parameter | estimate |
|---|---|---|---|
| sigma2 | NA | sigma2 | 7.432 |
| trial | T1 | variance | 7.974 |
| trial | T2 | variance | 16.734 |
| trial | T3 | variance | 3.715 |
| ar1m(range) | T1 | rho | 0.807 |
| ar1m(range) | T2 | rho | 0.213 |
| ar1m(range) | T3 | rho | 0.439 |
| ar1m(row) | T1 | rho | 0.208 |
| ar1m(row) | T2 | rho | 0.655 |
| ar1m(row) | T3 | rho | 0.505 |
The section column tells which trial each parameter
belongs to, and the estimates now track the truth in every 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. The two
can therefore 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
## 1.057576e+02 4.000000e+00 5.840487e-22As 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, verbose=FALSE)
knitr::kable(covparams_mmes(mLevels)[, c("factor", "section", "parameter", "estimate")],
digits=3)| factor | section | parameter | estimate |
|---|---|---|---|
| sigma2 | NA | sigma2 | 7.446 |
| trial | T1 | variance | 7.791 |
| trial | T2 | variance | 16.975 |
| trial | T3 | variance | 3.498 |
| ar1m(range) | T1 | rho | 0.806 |
| ar1m(range) | T2 | rho | 0.226 |
| ar1m(row) | T1 | rho | 0.242 |
| ar1m(row) | T2 | rho | 0.654 |
## used available
## 817 817The 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
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))).
residuals(fit) against range and row) show
whether spatial patterns remain.ar1m().henderson = FALSE (direct inversion) also
works, 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.