library(mandala)
library(SpATS)
library(dplyr)
library(ggplot2)Comparison: Mandala vs SpATS
P-spline spatial analysis for field trials
Overview
This tutorial compares the mandala package with the SpATS package for spatial analysis of field trials. Both packages implement P-spline (PSANOVA) spatial modeling, making them directly comparable for spatial trend correction.
NoteWhat you’ll learn
- How Mandala’s
pspline2D()compares to SpATS’sPSANOVA() - Side-by-side spatial model fitting
- Variogram diagnostics in both packages
- Fitted value comparison and validation
- When to use each package
Background
SpATS Package
- Purpose: Spatial analysis of field trials using P-splines
- Key Features:
- PSANOVA spatial smoothing
- Built-in spatial diagnostics
- Genotype as fixed or random
- Effective dimension calculations
- Primary Use: Spatial correction in breeding trials
Mandala Package
- Purpose: Comprehensive field trial analysis
- Key Features:
pspline2D()for SpATS-equivalent P-spline modeling- AR1×AR1 spatial correlation structures
- Interactive variogram diagnostics
- Integration with GBLUP and FA models
- Primary Use: Field trial analysis with multiple spatial options
Loading Packages
Example 1: Wheat Trial Data
We’ll use the built-in wheatdata from SpATS to demonstrate equivalence.
Data Setup
data("wheatdata")
# Prepare data for both packages
wheatdata$row.num <- as.numeric(wheatdata$row)
wheatdata$col.num <- as.numeric(wheatdata$col)
wheatdata$yield <- as.numeric(wheatdata$yield)
wheatdata$row <- as.factor(wheatdata$row)
wheatdata$col <- as.factor(wheatdata$col)
wheatdata$rep <- as.factor(wheatdata$rep)
str(wheatdata)'data.frame': 330 obs. of 9 variables:
$ yield : num 483 526 557 564 498 510 344 600 466 370 ...
$ geno : Factor w/ 107 levels "1","2","3","4",..: 4 10 17 16 21 32 33 34 72 74 ...
$ rep : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ row : Factor w/ 22 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ col : Factor w/ 15 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ rowcode: Factor w/ 2 levels "1","2": 2 2 1 2 2 1 2 2 1 2 ...
$ colcode: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$ row.num: num 1 2 3 4 5 6 7 8 9 10 ...
$ col.num: num 1 1 1 1 1 1 1 1 1 1 ...head(wheatdata) yield geno rep row col rowcode colcode row.num col.num
1 483 4 1 1 1 2 1 1 1
2 526 10 1 2 1 2 1 2 1
3 557 17 1 3 1 1 1 3 1
4 564 16 1 4 1 2 1 4 1
5 498 21 1 5 1 2 1 5 1
6 510 32 1 6 1 1 1 6 1Visualize Field Layout
ggplot(wheatdata, aes(x = col.num, y = row.num, fill = yield)) +
geom_tile() +
scale_fill_gradient2(low = "#2E6B8A", mid = "#FAF8F5", high = "#9B5B5B",
midpoint = mean(wheatdata$yield)) +
coord_equal() +
theme_minimal() +
labs(title = "Wheat Trial - Raw Yield",
x = "Column", y = "Row", fill = "Yield")
SpATS Analysis
spats_fit <- SpATS(
response = "yield",
genotype = "geno",
genotype.as.random = FALSE,
spatial = ~ PSANOVA(row.num, col.num,
nseg = c(10, 20),
degree = c(3, 3),
pord = c(2, 2)),
data = wheatdata
)Effective dimensions
-------------------------
It. Deviance f(row.num) f(col.num)f(row.num):col.numrow.num:f(col.num)f(row.num):f(col.num)
1 1089936.925340 2.301 5.377 2.212 4.669 13.093
2 2259.016100 2.092 6.669 1.312 3.887 10.805
3 2244.207735 1.929 8.244 0.731 3.234 8.987
4 2217.738426 1.824 10.179 0.275 2.702 7.796
5 2177.035873 1.788 11.936 0.065 2.341 7.516
6 2155.880798 1.768 12.643 0.013 2.181 8.128
7 2154.137260 1.704 12.762 0.003 2.102 8.794
8 2154.011058 1.632 12.777 0.001 2.055 9.299
9 2153.953448 1.566 12.779 0.000 2.026 9.675
10 2153.918651 1.508 12.779 0.000 2.009 9.957
11 2153.897223 1.459 12.780 0.000 2.000 10.168
12 2153.884068 1.418 12.780 0.000 1.997 10.329
13 2153.876082 1.386 12.780 0.000 1.997 10.450
14 2153.871289 1.361 12.780 0.000 2.000 10.543
15 2153.868429 1.342 12.780 0.000 2.003 10.615
16 2153.866719 1.328 12.780 0.000 2.008 10.670
17 2153.865682 1.318 12.781 0.000 2.013 10.714
18 2153.865039 1.310 12.781 0.000 2.017 10.748
Timings:
SpATS 0.109 seconds
All process 0.135 secondssummary(spats_fit)
Spatial analysis of trials with splines
Response: yield
Genotypes (as fixed): geno
Spatial: ~PSANOVA(row.num, col.num, nseg = c(10, 20), degree = c(3, 3), pord = c(2, 2))
Number of observations: 330
Number of missing data: 0
Effective dimension: 136.86
Deviance: 2153.865
Dimensions:
Effective Model Nominal Ratio Type
geno 106.0 106 106 1.00 F
Intercept 1.0 1 1 1.00 F
row.num 1.0 1 1 1.00 S
col.num 1.0 1 1 1.00 S
col.numrow.num 1.0 1 1 1.00 S
f(row.num) 1.3 11 11 0.12 S
f(col.num) 12.8 21 21 0.61 S
f(row.num):col.num 0.0 11 11 0.00 S
row.num:f(col.num) 2.0 21 21 0.10 S
f(row.num):f(col.num) 10.7 231 231 0.05 S
Total 136.9 405 405 0.34
Residual 193.1
Nobs 330
Type codes: F 'Fixed' R 'Random' S 'Smooth/Semiparametric'SpATS Diagnostics
plot(spats_fit)
SpATS Variogram
plot(variogram(spats_fit))
Mandala Analysis
P-spline Model (SpATS Equivalent)
# Mandala pspline2D - equivalent to SpATS PSANOVA
spline_fit <- mandala(
fixed = yield ~ geno,
random = ~ row.num + col.num + row.num:col.num +
pspline2D(row.num, col.num,
nseg = c(10, 20),
degree = c(3, 3),
pord = c(2, 2)),
data = wheatdata
)Initial data rows: 330
Final data rows after NA handling: 330 Explicit param_df Check:
name type vc_idx init lower upper origin
1 row.num var 1 9510.162 1e-06 Inf G
2 col.num var 2 9510.162 1e-06 Inf G
3 row.num:col.num var 3 9510.162 1e-06 Inf G
4 sf(row.num) var 4 9510.162 1e-06 Inf G
5 sf(col.num) var 5 9510.162 1e-06 Inf G
6 sf(row.num):sf(col.num) var 6 9510.162 1e-06 Inf G
7 R.sigma2 R_var 0 11887.703 1e-06 Inf R
--- Starting EM Warmup ---
EM Warmup: starting theta:
row.num 9510 col.num 9510 row.num:col.num 9510 sf(row.num) 9510 sf(col.num) 9510 sf(row.num):sf(col.num) 9510 R.sigma2 11890
EM iter 1: logLik=-1439.2375; theta: row.num 6558 col.num 7632 row.num:col.num 7345 sf(row.num) 7553 sf(col.num) 7838 sf(row.num):sf(col.num) 8992 R.sigma2 5211
EM iter 2: logLik=-1395.8907; theta: row.num 4729 col.num 6285 row.num:col.num 6368 sf(row.num) 5989 sf(col.num) 6476 sf(row.num):sf(col.num) 8558 R.sigma2 5500
EM iter 3: logLik=-1389.8149; theta: row.num 3428 col.num 5273 row.num:col.num 5337 sf(row.num) 4730 sf(col.num) 5358 sf(row.num):sf(col.num) 8102 R.sigma2 5377
EM iter 4: logLik=-1381.0810; theta: row.num 2544 col.num 4521 row.num:col.num 4543 sf(row.num) 3732 sf(col.num) 4448 sf(row.num):sf(col.num) 7666 R.sigma2 5410
EM iter 5: logLik=-1375.0492; theta: row.num 1930 col.num 3958 row.num:col.num 3910 sf(row.num) 2945 sf(col.num) 3708 sf(row.num):sf(col.num) 7249 R.sigma2 5431
EM Warmup: ending theta:
row.num 1930 col.num 3958 row.num:col.num 3910 sf(row.num) 2945 sf(col.num) 3708 sf(row.num):sf(col.num) 7249 R.sigma2 5431
--- EM Warmup Complete. Starting AI-REML ---Starting AI-REML logLik = -1370.3290
Iter LogLik Sigma2 DF wall Restrained
[STABILIZING] Adding nugget 62477 to V diagonal...
1 -1354.4069 4435.083 223 09:48:08 ( 0 restrained)
[STABILIZING] Adding nugget 54329 to V diagonal...
2 -1337.6924 3370.663 223 09:48:09 ( 0 restrained)
[STABILIZING] Adding nugget 45498 to V diagonal...
3 -1304.6970 0.000 223 09:48:09 ( 1 restrained)
[STABILIZING] Adding nugget 38674 to V diagonal...
4 -1301.3150 0.000 223 09:48:09 ( 1 restrained)
[STABILIZING] Adding nugget 34157 to V diagonal...
5 -1299.7197 0.000 223 09:48:10 ( 1 restrained)
[STABILIZING] Adding nugget 29957 to V diagonal...
6 -1299.5786 0.000 223 09:48:10 ( 1 restrained)
[STABILIZING] Adding nugget 25465 to V diagonal...
7 -1298.4056 0.000 223 09:48:10 ( 1 restrained)
[STABILIZING] Adding nugget 21646 to V diagonal...
8 -1298.4056 0.000 223 09:48:11 ( 1 restrained)
Converged at iter 8 by tolerance(s): delta_theta=0.00e+00, delta_loglik=0.00e+00
Main AI-REML Loop finished. Combining results.
solveMME_cpp: adding nugget 1 to MME diagonal (scaled from median diag)summary(spline_fit)Model statement:
mandala(fixed = yield ~ geno, random = ~row.num + col.num + row.num:col.num +
pspline2D(row.num, col.num, nseg = c(10, 20), degree = c(3,
3), pord = c(2, 2)), data = wheatdata)
Variance Components:
component estimate std.error z.ratio bound %ch
row.num 608.638806 1098.2321 5.541987e-01 P NA
col.num 4975.620738 3144.3384 1.582406e+00 P NA
row.num:col.num 1834.368421 450487.7576 4.071961e-03 P NA
sf(row.num) 855.627687 640.1351 1.336636e+00 P NA
sf(col.num) 1450.689594 1103.8222 1.314242e+00 P NA
sf(row.num):sf(col.num) 3283.372796 1733.3344 1.894252e+00 P NA
R.sigma2 0.000001 450487.7576 2.219816e-12 B NA
Fixed Effects (BLUEs) [first 5]:
effect estimate std.error z.ratio
(Intercept) 528.536279 0.3323532 1590.285018
geno2 49.322810 0.5446576 90.557464
geno3 15.190728 0.5450680 27.869417
geno4 -3.484060 0.5532713 -6.297201
geno5 3.021166 0.5392727 5.602298
Converged: TRUE | Iterations: 8
Random Effects (BLUPs) [first 5]:
random level estimate std.error z.ratio
row.num 1 4.6462522 0.4284627 10.844007
row.num 2 17.6825541 0.5847960 30.237133
row.num 3 12.9493980 0.5452031 23.751510
row.num 4 0.9004335 0.5298191 1.699511
row.num 5 33.2697537 0.5453053 61.011243
logLik: -1298.406 AIC: 2610.811 BIC: 2634.661 logLik_Trunc: -1093.482Mandala Spatial Diagnostics
# Spatial fitted value plot
plot_mandala_spatial(spline_fit, wheatdata,
response = "yield",
geno = "geno",
row = "row",
col = "col")
Notes:
- The predictions are obtained by averaging across the hypertable
calculated from model terms constructed solely from factors in
the averaging and classify sets.
- The simple averaging set: (none)
- The ignored set: row.num:col.num, sf(row.num), sf(col.num), sf(row.num):sf(col.num)
Mandala Variogram Options
# Static 3D mesh
mandala_variogram(spline_fit, data = wheatdata,
lag_width = 1, plot_type = "3d")
$variogram_grid
[,1] [,2] [,3] [,4] [,5]
[1,] 0.000000e+00 1.415379e-09 1.430687e-09 1.333888e-09 1.240534e-09
[2,] 1.235693e-09 1.292795e-09 1.294928e-09 1.310815e-09 1.203838e-09
[3,] 1.160733e-09 1.268368e-09 1.202739e-09 1.282324e-09 1.245733e-09
[4,] 1.223067e-09 1.233555e-09 1.267561e-09 1.355977e-09 1.108196e-09
[5,] 1.211086e-09 1.185239e-09 1.152668e-09 1.213955e-09 1.251256e-09
[6,] 1.226811e-09 1.226442e-09 1.251363e-09 1.233839e-09 1.303506e-09
[7,] 1.171979e-09 1.298923e-09 1.197121e-09 1.269247e-09 1.254007e-09
[8,] 1.352956e-09 1.270902e-09 1.247089e-09 1.216865e-09 1.365988e-09
[9,] 1.191194e-09 1.222003e-09 1.329110e-09 1.293242e-09 1.182023e-09
[10,] 1.517854e-09 1.299598e-09 1.222671e-09 1.392423e-09 1.272549e-09
[11,] 1.444526e-09 1.170462e-09 1.211212e-09 1.216454e-09 1.379795e-09
[12,] 1.436770e-09 1.257300e-09 1.424160e-09 1.363518e-09 1.274405e-09
[13,] 1.555287e-09 1.205762e-09 1.320926e-09 1.280443e-09 1.455255e-09
[14,] 1.505886e-09 1.245308e-09 1.284572e-09 1.238574e-09 1.472073e-09
[15,] 1.710579e-09 1.227555e-09 1.486569e-09 1.371403e-09 1.264984e-09
[16,] 1.452844e-09 1.316361e-09 1.326178e-09 1.487645e-09 1.468671e-09
[17,] 1.628715e-09 1.401217e-09 1.435485e-09 1.524233e-09 1.434129e-09
[18,] 1.644693e-09 1.517678e-09 1.366269e-09 1.659635e-09 1.510176e-09
[19,] 2.043148e-09 1.523457e-09 1.369255e-09 1.453316e-09 1.744001e-09
[20,] 2.070180e-09 1.415160e-09 1.117956e-09 1.462838e-09 1.531079e-09
[21,] 2.235135e-09 1.448555e-09 1.430268e-09 1.414585e-09 1.912342e-09
[22,] 3.010280e-09 1.671620e-09 1.529402e-09 1.990422e-09 2.042772e-09
[,6] [,7] [,8] [,9] [,10]
[1,] 1.372455e-09 1.641692e-09 1.440314e-09 1.187832e-09 1.504353e-09
[2,] 1.389735e-09 1.173762e-09 1.394752e-09 1.247825e-09 1.301921e-09
[3,] 1.333247e-09 1.323334e-09 1.321732e-09 1.289119e-09 1.300471e-09
[4,] 1.430812e-09 1.395821e-09 1.397911e-09 1.241286e-09 1.361014e-09
[5,] 1.298608e-09 1.244489e-09 1.107706e-09 1.280326e-09 1.304627e-09
[6,] 1.376957e-09 1.254288e-09 1.332091e-09 1.199960e-09 1.181521e-09
[7,] 1.234561e-09 1.229512e-09 1.418167e-09 1.199028e-09 1.149202e-09
[8,] 1.283370e-09 1.245107e-09 1.294881e-09 1.185669e-09 1.240078e-09
[9,] 1.315455e-09 1.334104e-09 1.355438e-09 1.266213e-09 1.165576e-09
[10,] 1.339337e-09 1.397790e-09 1.378745e-09 1.264569e-09 1.170504e-09
[11,] 1.358070e-09 1.392566e-09 1.212808e-09 1.399050e-09 1.297269e-09
[12,] 1.299221e-09 1.427130e-09 1.210902e-09 1.141366e-09 1.315217e-09
[13,] 1.458664e-09 1.338722e-09 1.321035e-09 1.408979e-09 1.389589e-09
[14,] 1.208938e-09 1.475392e-09 1.216206e-09 1.493588e-09 1.320839e-09
[15,] 1.548171e-09 1.305491e-09 1.324200e-09 1.476568e-09 1.210071e-09
[16,] 1.537234e-09 1.567284e-09 1.371132e-09 1.507973e-09 1.326011e-09
[17,] 1.550678e-09 1.407743e-09 1.606780e-09 1.421875e-09 1.107474e-09
[18,] 1.442106e-09 1.844614e-09 1.341661e-09 1.637447e-09 1.538289e-09
[19,] 1.926919e-09 1.466432e-09 1.553880e-09 1.772400e-09 2.092342e-09
[20,] 1.656689e-09 1.334946e-09 1.105715e-09 1.970283e-09 1.380275e-09
[21,] 2.106461e-09 1.180584e-09 1.773978e-09 1.800105e-09 1.425896e-09
[22,] 2.419452e-09 1.151728e-09 1.967037e-09 2.944455e-09 1.154318e-09
[,11] [,12] [,13] [,14] [,15]
[1,] 1.417020e-09 1.593466e-09 1.197541e-09 1.756001e-09 2.460439e-09
[2,] 1.293161e-09 1.223192e-09 1.127496e-09 1.358754e-09 2.417085e-09
[3,] 1.221575e-09 1.425743e-09 1.203372e-09 1.686548e-09 2.255278e-09
[4,] 1.209075e-09 1.159023e-09 1.105120e-09 1.275272e-09 1.983199e-09
[5,] 1.301165e-09 1.413059e-09 1.339680e-09 1.883536e-09 2.215789e-09
[6,] 1.182431e-09 1.146402e-09 1.365616e-09 1.033900e-09 2.182882e-09
[7,] 1.496602e-09 1.454652e-09 1.208977e-09 1.495779e-09 2.083055e-09
[8,] 1.050996e-09 1.378295e-09 1.237963e-09 1.486452e-09 2.261698e-09
[9,] 1.112376e-09 1.423636e-09 1.533661e-09 1.624115e-09 1.819842e-09
[10,] 1.132530e-09 1.474030e-09 1.305355e-09 1.497481e-09 1.512854e-09
[11,] 1.111780e-09 1.224051e-09 1.480681e-09 1.389746e-09 1.546400e-09
[12,] 1.161898e-09 1.171641e-09 1.325254e-09 1.747512e-09 1.976757e-09
[13,] 1.244084e-09 1.245158e-09 1.584193e-09 1.372380e-09 1.460893e-09
[14,] 8.820791e-10 1.323314e-09 1.431705e-09 1.538489e-09 9.763371e-10
[15,] 1.333289e-09 1.126313e-09 1.434829e-09 1.715955e-09 1.536108e-09
[16,] 9.797865e-10 1.682397e-09 1.324397e-09 1.462896e-09 1.353299e-09
[17,] 1.187142e-09 1.556004e-09 1.496400e-09 1.047419e-09 9.779765e-10
[18,] 1.049745e-09 1.957633e-09 1.302816e-09 1.781815e-09 1.021572e-09
[19,] 1.336043e-09 1.472882e-09 2.070271e-09 1.503229e-09 8.233582e-10
[20,] 1.423774e-09 1.876202e-09 2.206044e-09 3.237776e-09 9.427809e-10
[21,] 1.415285e-09 1.788031e-09 2.420264e-09 1.275346e-09 NA
[22,] 1.772348e-09 2.250581e-09 2.800100e-09 NA NA
$counts
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 1 308 286 264 242 220 198 176 154 132 110 88 66
[2,] 315 588 546 504 462 420 378 336 294 252 210 168 126
[3,] 300 560 520 480 440 400 360 320 280 240 200 160 120
[4,] 285 532 494 456 418 380 342 304 266 228 190 152 114
[5,] 270 504 468 432 396 360 324 288 252 216 180 144 108
[6,] 255 476 442 408 374 340 306 272 238 204 170 136 102
[7,] 240 448 416 384 352 320 288 256 224 192 160 128 96
[8,] 225 420 390 360 330 300 270 240 210 180 150 120 90
[9,] 210 392 364 336 308 280 252 224 196 168 140 112 84
[10,] 195 364 338 312 286 260 234 208 182 156 130 104 78
[11,] 180 336 312 288 264 240 216 192 168 144 120 96 72
[12,] 165 308 286 264 242 220 198 176 154 132 110 88 66
[13,] 150 280 260 240 220 200 180 160 140 120 100 80 60
[14,] 135 252 234 216 198 180 162 144 126 108 90 72 54
[15,] 120 224 208 192 176 160 144 128 112 96 80 64 48
[16,] 105 196 182 168 154 140 126 112 98 84 70 56 42
[17,] 90 168 156 144 132 120 108 96 84 72 60 48 36
[18,] 75 140 130 120 110 100 90 80 70 60 50 40 30
[19,] 60 112 104 96 88 80 72 64 56 48 40 32 24
[20,] 45 84 78 72 66 60 54 48 42 36 30 24 18
[21,] 30 56 52 48 44 40 36 32 28 24 20 16 12
[22,] 15 28 26 24 22 20 18 16 14 12 10 8 6
[,14] [,15]
[1,] 44 22
[2,] 84 42
[3,] 80 40
[4,] 76 38
[5,] 72 36
[6,] 68 34
[7,] 64 32
[8,] 60 30
[9,] 56 28
[10,] 52 26
[11,] 48 24
[12,] 44 22
[13,] 40 20
[14,] 36 18
[15,] 32 16
[16,] 28 14
[17,] 24 12
[18,] 20 10
[19,] 16 8
[20,] 12 6
[21,] 8 4
[22,] 4 2
$row_lags
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
$col_lags
[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14# 2D lag vs semivariance
mandala_variogram(spline_fit, data = wheatdata,
lag_width = 2, plot_type = "2d")
lag_distance semivariance n_pairs
1 1 1.303695e-09 1797
2 3 1.262350e-09 4697
3 5 1.263675e-09 7015
4 7 1.286642e-09 7599
5 9 1.294353e-09 8583
6 11 1.302808e-09 6867
7 13 1.304620e-09 6751
8 15 1.423123e-09 4535
9 17 1.451389e-09 3169
10 19 1.495849e-09 2093
11 21 1.639754e-09 977
12 23 1.846818e-09 186
13 25 2.007106e-09 16# Interactive Plotly variogram
library(plotly)
mandala_variogram(spline_fit, data = wheatdata,
lag_width = 1,
plot_type = "dynamic",
empirical_na = "smooth",
overlay = "ar1xar1")Alternative: AR1 Spatial Model
Mandala also supports AR1 correlation structures, which SpATS does not:
# AR1×AR1 spatial correlation
ar1_fit <- mandala(
fixed = yield ~ geno,
random = ~ ar1(row) + ar1(col) + ar1(row):ar1(col),
data = wheatdata
)Initial data rows: 330
Final data rows after NA handling: 330
Zlst[[3]] dimensions: 330x330
Explicit param_df Check:
name type vc_idx init lower upper origin
1 ar1(row) var 1 9510.162 1.0e-06 Inf G
2 rho.ar1(row) rho 1 0.500 -9.5e-01 0.95 G
3 ar1(col) var 2 9510.162 1.0e-06 Inf G
4 rho.ar1(col) rho 2 0.500 -9.5e-01 0.95 G
5 ar1(row):ar1(col) var 3 9510.162 1.0e-06 Inf G
6 rho.row.col.1 rho1 3 0.500 -9.5e-01 0.95 G
7 rho.row.col.2 rho2 3 0.500 -9.5e-01 0.95 G
8 R.sigma2 R_var 0 11887.703 1.0e-06 Inf R
--- Starting EM Warmup ---
EM Warmup: starting theta:
ar1(row) 9510 rho.ar1(row) 0.5 ar1(col) 9510 rho.ar1(col) 0.5 ar1(row):ar1(col) 9510 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 11890
EM iter 1: logLik=-1433.8367; theta: ar1(row) 5607 rho.ar1(row) 0.5 ar1(col) 8766 rho.ar1(col) 0.5 ar1(row):ar1(col) 7368 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5382
EM iter 2: logLik=-1388.7742; theta: ar1(row) 3735 rho.ar1(row) 0.5 ar1(col) 8281 rho.ar1(col) 0.5 ar1(row):ar1(col) 6511 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5759
EM iter 3: logLik=-1383.5202; theta: ar1(row) 2535 rho.ar1(row) 0.5 ar1(col) 8077 rho.ar1(col) 0.5 ar1(row):ar1(col) 5574 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5597
EM iter 4: logLik=-1374.4362; theta: ar1(row) 1824 rho.ar1(row) 0.5 ar1(col) 7994 rho.ar1(col) 0.5 ar1(row):ar1(col) 4844 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5657
EM iter 5: logLik=-1368.3727; theta: ar1(row) 1378 rho.ar1(row) 0.5 ar1(col) 7989 rho.ar1(col) 0.5 ar1(row):ar1(col) 4248 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5696
EM Warmup: ending theta:
ar1(row) 1378 rho.ar1(row) 0.5 ar1(col) 7989 rho.ar1(col) 0.5 ar1(row):ar1(col) 4248 rho.row.col.1 0.5 rho.row.col.2 0.5 R.sigma2 5696
--- EM Warmup Complete. Starting AI-REML ---Starting AI-REML logLik = -1363.3120
Iter LogLik Sigma2 DF wall Restrained
1 -1342.2029 4329.129 223 09:48:23 ( 1 restrained)
2 -1341.1299 3290.138 223 09:48:24 ( 0 restrained)
3 -1319.1257 1540.043 223 09:48:24 ( 0 restrained)
4 -1306.8953 759.618 223 09:48:24 ( 1 restrained)
5 -1305.3234 768.621 223 09:48:24 ( 0 restrained)
6 -1304.5211 1185.441 223 09:48:25 ( 1 restrained)
7 -1304.5211 1185.441 223 09:48:25 ( 1 restrained)
Converged at iter 7 by tolerance(s): delta_theta=0.00e+00, delta_loglik=0.00e+00
Main AI-REML Loop finished. Combining results.
solveMME_cpp: adding nugget 2.19553e-09 to MME diagonal (scaled from median diag)summary(ar1_fit)Model statement:
mandala(fixed = yield ~ geno, random = ~ar1(row) + ar1(col) +
ar1(row):ar1(col), data = wheatdata)
Variance Components:
component estimate std.error z.ratio bound %ch
ar1(row) 809.1309566 388.5417862 2.082481 P NA
rho.ar1(row) 0.3656464 0.2634864 1.387724 P NA
ar1(col) 9733.0895992 1617.0734685 6.018953 P NA
rho.ar1(col) 0.9500000 0.0752517 12.624300 B NA
ar1(row):ar1(col) 2089.2092851 540.4125308 3.865953 P NA
rho.row.col.1 0.2113991 0.1805995 1.170541 P NA
rho.row.col.2 0.6646716 0.1307389 5.083961 P NA
R.sigma2 1185.4409293 462.0376181 2.565681 P NA
Fixed Effects (BLUEs) [first 5]:
effect estimate std.error z.ratio
(Intercept) 560.93380 90.44528 6.201913
geno2 -37.26906 39.78923 -0.936662
geno3 -79.61364 42.83845 -1.858462
geno4 -124.64943 42.60802 -2.925492
geno5 -97.46382 42.18750 -2.310253
Converged: TRUE | Iterations: 7
Random Effects (BLUPs) [first 5]:
random level estimate std.error z.ratio
ar1(row) 1 -6.052136 16.25523 -0.3723193
ar1(row) 2 -11.321963 15.93199 -0.7106434
ar1(row) 3 -25.261593 16.32902 -1.5470363
ar1(row) 4 -34.106965 16.43449 -2.0753284
ar1(row) 5 -6.774268 16.55608 -0.4091710
logLik: -1304.521 AIC: 2625.042 BIC: 2652.300 logLik_Trunc: -1099.598# Compare variograms
mandala_variogram(ar1_fit, data = wheatdata,
lag_width = 1,
plot_type = "dynamic",
overlay = "ar1xar1")Comparing Fitted Values
When both models are fit, you can validate their equivalence:
# Extract fitted values
fitted_spats <- as.numeric(spats_fit$fitted)
fitted_mandala <- as.numeric(spline_fit$fitted)
# Correlation
cor_fit <- cor(fitted_spats, fitted_mandala)
cat(sprintf("Correlation: %.4f\n", cor_fit))Correlation: 0.9695# RMSE of difference
rmse_diff <- sqrt(mean((fitted_spats - fitted_mandala)^2))
cat(sprintf("RMSE difference: %.4f\n", rmse_diff))RMSE difference: 37.7682# Scatter plot
comp_df <- data.frame(
SpATS = fitted_spats,
Mandala = fitted_mandala
)
ggplot(comp_df, aes(x = SpATS, y = Mandala)) +
geom_point(alpha = 0.5, color = "#5D4E6D") +
geom_abline(slope = 1, intercept = 0, color = "#9B5B5B", linetype = "dashed") +
theme_minimal() +
labs(title = "SpATS vs Mandala Fitted Values",
x = "SpATS Fitted", y = "Mandala Fitted")
Key Differences Summary
| Aspect | Mandala | SpATS |
|---|---|---|
| P-spline syntax | pspline2D() |
PSANOVA() |
| AR1 spatial | ar1(row) + ar1(col) |
Not supported |
| Variograms | 2D, 3D, interactive | 2D only |
| Genotype effects | Fixed or random | Fixed or random |
| Genomic models | GM() for GBLUP |
Not supported |
| FA G×E models | FA() |
Not supported |
| MET integration | by(env): syntax |
Single trial focus |
| Diagnostics | mandala_diagnostic_plots() |
plot() |
When to Use Each
Use Mandala when:
- Need both P-spline AND AR1 options
- Integrating spatial with GBLUP
- Multi-environment trials with spatial
- Interactive variogram exploration
- Factor-analytic G×E + spatial
Use SpATS when:
- Single-trial P-spline analysis
- Established SpATS workflow
- Effective dimensions important
- Need SpATS-specific outputs
Additional Resources
- SpATS: https://cran.r-project.org/package=SpATS
- Mandala: https://listoagriculture.com/products
- Rodríguez-Álvarez et al. (2018) - SpATS methodology
Session Information
sessionInfo()R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] plotly_4.12.0 ggplot2_4.0.2 dplyr_1.2.0 SpATS_1.0-19 mandala_1.0.1
loaded via a namespace (and not attached):
[1] Matrix_1.7-4 gtable_0.3.6 jsonlite_2.0.0
[4] compiler_4.5.2 maps_3.4.3 tidyselect_1.2.1
[7] Rcpp_1.1.1 splines2_0.5.4 tidyr_1.3.2
[10] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0
[13] lattice_0.22-7 R6_2.6.1 labeling_0.4.3
[16] generics_0.1.4 knitr_1.51 htmlwidgets_1.6.4
[19] dotCall64_1.2 tibble_3.3.1 pillar_1.11.1
[22] RColorBrewer_1.1-3 rlang_1.1.7 xfun_0.56
[25] S7_0.2.1 lazyeval_0.2.2 otel_0.2.0
[28] viridisLite_0.4.2 cli_3.6.5 withr_3.0.2
[31] magrittr_2.0.4 crosstalk_1.2.2 digest_0.6.39
[34] grid_4.5.2 spam_2.11-3 lifecycle_1.0.5
[37] fields_17.1 vctrs_0.7.1 evaluate_1.0.5
[40] glue_1.8.0 data.table_1.18.2.1 farver_2.1.2
[43] purrr_1.2.1 httr_1.4.7 rmarkdown_2.30
[46] tools_4.5.2 pkgconfig_2.0.3 htmltools_0.5.9