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)\). 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 F1
All 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)
)
## Solver selected: ldlt
## Solver selected: ldlt
## Solver selected: ldlt
## Solver selected: ldlt
## Solver selected: ldlt
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.981 | 0.890 |
| rowcol | -119.592 | 241.185 | 7.468 | 0.908 |
| spline | -108.224 | 218.449 | 6.899 | 0.915 |
| ar1Res | -81.485 | 164.971 | 7.356 | 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.
knitr::kable(covparams_mmes(fits$ar1Fld)[, c("factor", "parameter", "estimate")], digits=3)
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 7.232 |
| ar1m(rangef) | variance | 7.133 |
| ar1m(rangef) | rho | 0.666 |
| ar1m(rowf) | rho | 0.469 |
| sigma2 | sigma2 | 1.009 |
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)
## Solver selected: cholmod
summary(mix)$varcomp
## term parameter estimate StdError Zratio
## 1 vsm(ism(R)) sigma2 100.29557789 87.23660 1.1496960283
## 2 vsm(ism(C)) sigma2 180.28017638 181.56042 0.9929486638
## 3 vsm(ism(M$fC)) sigma2 0.53564980 1478.86764 0.0003622027
## 4 vsm(ism(M$fR)) sigma2 0.02082692 55.25000 0.0003769578
## 5 vsm(ism(M$fC.R)) sigma2 0.01522377 35.37195 0.0004303911
## 6 vsm(ism(M$C.fR)) sigma2 0.01843282 18.99536 0.0009703859
## 7 vsm(ism(M$fC.fR)) sigma2 0.01444448 59.60175 0.0002423499
## 8 vsm(ism(units)) sigma2 502.76569161 108.59251 4.6298377445
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.
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 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, verbose=FALSE)
## Solver selected: ldlt
knitr::kable(covparams_mmes(mShared)[, c("factor", "parameter", "estimate")], digits=3)
| factor | parameter | estimate |
|---|---|---|
| sigma2 | sigma2 | 7.416 |
| dsm(trial) | variance[T1] | 6.820 |
| dsm(trial) | variance[T2] | 20.797 |
| dsm(trial) | variance[T3] | 4.435 |
| ar1m(range) | rho | 0.589 |
| ar1m(row) | rho | 0.480 |
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, verbose=FALSE)
## Solver selected: ldlt
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.980 |
| trial | T2 | variance | 16.732 |
| trial | T3 | variance | 3.715 |
| ar1m(range) | T1 | rho | 0.807 |
| ar1m(range) | T2 | rho | 0.212 |
| 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.840525e-22
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, verbose=FALSE)
## Solver selected: ldlt
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.798 |
| trial | T2 | variance | 16.974 |
| trial | T3 | variance | 3.497 |
| ar1m(range) | T1 | rho | 0.806 |
| ar1m(range) | T2 | rho | 0.226 |
| ar1m(row) | T1 | rho | 0.242 |
| ar1m(row) | T2 | rho | 0.654 |
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))).
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.