Spatial Analysis

Spatial field variation is common in breeding trials. Mandala can fit simple row-column mixed models and then provide spatial diagnostic plots from the fitted model. The examples use a small augmented single-location trial.

Data Preparation

library(mandala)

aug <- read.csv(
  system.file("extdata", "augmented_single_n200.csv", package = "mandala"),
  stringsAsFactors = FALSE,
  check.names = FALSE
)

for (v in c("geno", "row", "col")) {
  aug[[v]] <- factor(aug[[v]])
}

aug$yld <- as.numeric(aug$yld)

This dataset has many unreplicated genotypes and a smaller set of replicated check genotypes. Row and column describe the field position of each plot.

data.frame(
  plots = nrow(aug),
  genotypes = length(unique(aug$geno)),
  replicated_genotypes = sum(table(aug$geno) > 1),
  unreplicated_genotypes = sum(table(aug$geno) == 1),
  rows = length(unique(aug$row)),
  columns = length(unique(aug$col))
)
#>   plots genotypes replicated_genotypes unreplicated_genotypes rows columns
#> 1   200       160                   20                    140   14      15

Model 1: Genotype Baseline

The first model estimates genotype effects without using field-position information. This gives a simple baseline for judging spatial adjustment.

fit_base <- mandala(
  yld ~ 1,
  ~ geno,
  data = aug,
  verbose = FALSE
)

summary(fit_base)
#> Model statement:
#> mandala(fixed = yld ~ 1, random = ~geno, data = aug, verbose = FALSE)
#> 
#> Variance Components:
#>  component  estimate  std.error  z.ratio bound %ch
#>       geno 0.6686315 0.14727894 4.539899     P  NA
#>   R.sigma2 0.4571173 0.09963492 4.587923     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect estimate std.error  z.ratio
#>  (Intercept) 7.149171 0.0819997 87.18534
#> 
#> Converged: TRUE  |  Iterations: 5
#> 
#> Model Notes:
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>  random level   estimate std.error     z.ratio
#>    geno    G1 -0.0259458 0.3585428 -0.07236459
#>    geno   G10 -0.6446108 0.3585428 -1.79786301
#>    geno  G100  0.4141553 0.5233297  0.79138499
#>    geno  G101 -0.1874388 0.5233297 -0.35816584
#>    geno  G102  0.6048431 0.5233297  1.15575913
#> 
#> logLik: -286.519   AIC: 577.037   BIC: 583.624   logLik_Trunc: -103.650

Model 2: Row-Column Spatial Adjustment

The second model adds row and column effects. This is a simple row-column spatial adjustment: it does not assume a full AR1 residual model, but it is a useful first step for field trends that align with rows or columns.

fit_spatial <- mandala(
  yld ~ 1,
  ~ geno + row + col,
  data = aug,
  verbose = FALSE
)

summary(fit_spatial)
#> Model statement:
#> mandala(fixed = yld ~ 1, random = ~geno + row + col, data = aug, 
#>     verbose = FALSE)
#> 
#> Variance Components:
#>  component   estimate  std.error  z.ratio bound %ch
#>       geno 0.42983930 0.06126367 7.016218     P  NA
#>        row 0.29599940 0.12308129 2.404910     P  NA
#>        col 0.30309231 0.12172672 2.489941     P  NA
#>   R.sigma2 0.06052303 0.01945609 3.110749     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect estimate std.error  z.ratio
#>  (Intercept) 7.186965 0.2108871 34.07968
#> 
#> Converged: TRUE  |  Iterations: 11
#> 
#> Model Notes:
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>  random level    estimate std.error    z.ratio
#>    geno    G1 -0.34350610 0.1730044 -1.9855338
#>    geno   G10 -0.69072325 0.1787749 -3.8636481
#>    geno  G100  0.42780962 0.2863221  1.4941552
#>    geno  G101  0.06766663 0.2869923  0.2357785
#>    geno  G102  0.37495894 0.2763080  1.3570326
#> 
#> logLik: -222.107   AIC: 452.214   BIC: 465.387   logLik_Trunc: -39.238

Alternative Spatial Models

Mandala can also fit residual AR1 structures through R_formula. A one-margin AR1 model is written by crossing an identity margin with an AR1 margin. A separable two-dimensional residual model is written as ar1(row):ar1(col).

fit_ar1_col <- mandala(
  yld ~ 1,
  ~ geno,
  data = aug,
  R_formula = ~ id(row):ar1(col),
  verbose = FALSE
)

fit_ar1_row <- mandala(
  yld ~ 1,
  ~ geno,
  data = aug,
  R_formula = ~ ar1(row):id(col),
  verbose = FALSE
)

fit_ar1_2d <- mandala(
  yld ~ 1,
  ~ geno,
  data = aug,
  R_formula = ~ ar1(row):ar1(col),
  verbose = FALSE
)

rbind(
  ar1_col = fit_ar1_col$varcomp,
  ar1_row = fit_ar1_row$varcomp,
  ar1_2d = fit_ar1_2d$varcomp
)
#>           component  estimate  std.error   z.ratio bound %ch
#> ar1_col.1      geno 0.3678125 0.07474901  4.920634     P  NA
#> ar1_col.2  R.sigma2 0.7903385 0.07467178 10.584166     P  NA
#> ar1_col.3 R.rho.col 0.8550000 0.34497086  2.478470     P  NA
#> ar1_row.1      geno 0.4952079 0.08177405  6.055807     P  NA
#> ar1_row.2  R.sigma2 0.5176285 0.05140675 10.069272     P  NA
#> ar1_row.3 R.rho.row 0.8575177 0.25926688  3.307471     P  NA
#> ar1_2d.1       geno 0.4007287 0.06639072  6.035915     P  NA
#> ar1_2d.2   R.sigma2 0.3149987 0.05321163  5.919735     P  NA
#> ar1_2d.3  R.rho.row 0.4457830 0.06824734  6.531873     P  NA
#> ar1_2d.4  R.rho.col 0.7079273 0.13492079  5.246984     P  NA

A two-dimensional spline can be fitted as a random spatial smooth. The compact call pspline2D(row, col) creates internal numeric coordinates and adds the spatial plane used by Mandala’s spline model.

fit_spline <- mandala(
  yld ~ 1,
  ~ geno + pspline2D(row, col),
  data = aug,
  verbose = FALSE
)

fit_spline$varcomp
#>               component    estimate  std.error    z.ratio bound %ch
#> 1                  geno 0.445503646 0.06530119 6.82228944     P  NA
#> 2            f(row.num) 0.357005638 0.41565919 0.85889027     P  NA
#> 3            f(col.num) 0.224737496 0.27146797 0.82786007     P  NA
#> 4    f(row.num):col.num 0.006806629 0.17435814 0.03903820     P  NA
#> 5    row.num:f(col.num) 0.006806629 0.07019476 0.09696776     P  NA
#> 6 f(row.num):f(col.num) 0.010268015 0.02003222 0.51257493     P  NA
#> 7                   row 0.006806629 0.01855800 0.36677592     P  NA
#> 8                   col 0.064312123 0.04266959 1.50721224     P  NA
#> 9              R.sigma2 0.056571369 0.01992362 2.83941213     P  NA

Model Evaluation

For models fitted to the same response and data, logLik, AIC, and BIC provide a quick model comparison. Lower AIC or BIC is preferred, while a higher logLik indicates better likelihood fit.

fit_list <- list(
  baseline = fit_base,
  row_col = fit_spatial,
  ar1_col = fit_ar1_col,
  ar1_row = fit_ar1_row,
  ar1_row_col = fit_ar1_2d,
  pspline2D = fit_spline
)

model_compare <- data.frame(
  model = c("baseline", "row_col", "ar1_col", "ar1_row", "ar1_row_col", "pspline2D"),
  logLik = c(
    fit_base$logLik,
    fit_spatial$logLik,
    fit_ar1_col$logLik,
    fit_ar1_row$logLik,
    fit_ar1_2d$logLik,
    fit_spline$logLik
  ),
  AIC = c(
    fit_base$AIC,
    fit_spatial$AIC,
    fit_ar1_col$AIC,
    fit_ar1_row$AIC,
    fit_ar1_2d$AIC,
    fit_spline$AIC
  ),
  BIC = c(
    fit_base$BIC,
    fit_spatial$BIC,
    fit_ar1_col$BIC,
    fit_ar1_row$BIC,
    fit_ar1_2d$BIC,
    fit_spline$BIC
  )
)

model_compare
#>         model    logLik      AIC      BIC
#> 1    baseline -286.5187 577.0375 583.6241
#> 2     row_col -222.1069 452.2138 465.3870
#> 3     ar1_col -240.7168 487.4337 497.3136
#> 4     ar1_row -244.3509 494.7018 504.5818
#> 5 ar1_row_col -230.0119 468.0237 481.1970
#> 6   pspline2D -212.7031 443.4061 472.9092
best_model <- model_compare$model[which.min(model_compare$AIC)]
best_fit <- fit_list[[best_model]]
best_model
#> [1] "pspline2D"

In this augmented example, the row-column and spline models should improve the fit over the baseline model because the dataset contains a clear field trend. The best model is still a statistical choice: compare fit, parsimony, residual diagnostics, and whether the spatial terms make sense for the field layout.

Fixed-Effect Test

This example has only an intercept fixed effect, so the fixed-effect test is mostly a workflow demonstration. In trials with treatments, environments, or covariates, the same function tests those fixed terms.

mandala_fixed_tests(best_fit, type = "selected", denDF = "residual")
#> 
#> Fixed-effect term tests
#> -----------------------
#> Method: selected 
#> Denominator df: residual 
#> 
#>                 Df denDF       F.inc        Wald           Pr status
#> (Intercept)      1   196 820.8275192 820.8275192 5.396274e-72     ok
#> row.num          1   196  30.8677667  30.8677667 8.928691e-08     ok
#> col.num          1   196  12.7333661  12.7333661 4.516097e-04     ok
#> row.num:col.num  1   196   0.2266391   0.2266391 6.345574e-01     ok
#> 
#> Note:
#> - Residual df were used for all fixed-effect terms.

Standard Diagnostics

mandala_diagnostic_plots(best_fit, response = "yld")

Spatial Diagnostics

mandala_spatial_diagnostics() shows observed values, fitted values, residuals, and a smoothed spatial surface. This is often the fastest way to see whether spatial pattern remains after fitting the model.

mandala_spatial_diagnostics(
  best_fit,
  data = aug,
  response = "yld",
  row = "row",
  col = "col"
)

Empirical Variogram

The empirical variogram summarizes how residual similarity changes with field distance. Strong structure in the variogram can indicate remaining spatial pattern.

vg_2d <- mandala_variogram(
  best_fit,
  data = aug,
  response = "yld",
  row = "row",
  col = "col",
  plot_type = "2d"
)

The same empirical variogram can also be displayed as a 3-D surface.

vg_3d <- mandala_variogram(
  best_fit,
  data = aug,
  response = "yld",
  row = "row",
  col = "col",
  plot_type = "3d"
)

Genotype Prediction

Finally, predict adjusted genotype values from the spatial model.

pred_aug <- mandala_predict(best_fit, "geno", present = TRUE, verbose = FALSE)
head(pred_aug)
#>   geno predicted_value std_error
#> 1   G1        6.884606 0.1725286
#> 2  G10        6.492034 0.1826274
#> 3 G100        7.522207 0.2754595
#> 4 G101        7.260529 0.2822562
#> 5 G102        7.467809 0.2672455
#> 6 G103        7.029064 0.2638560

Quick SpATS Comparison

If SpATS is installed, a comparable genotype-random spatial spline can be fitted and compared with Mandala’s spline model through fitted values.

if (requireNamespace("SpATS", quietly = TRUE)) {
  aug$row_num <- as.numeric(as.character(aug$row))
  aug$col_num <- as.numeric(as.character(aug$col))

  fit_spats <- SpATS::SpATS(
    response = "yld",
    genotype = "geno",
    genotype.as.random = TRUE,
    spatial = ~ SpATS::SAP(row_num, col_num, nseg = c(8, 8), degree = 3, pord = 2),
    data = aug,
    control = SpATS::controlSpATS(tolerance = 1e-03, monitoring = 0)
  )

  fitted_spats <- as.numeric(fit_spats$fitted)
  fitted_mandala <- as.numeric(fit_spline$fitted)
  comp_df <- data.frame(
    yld = aug$yld,
    fitted_spats = fitted_spats,
    fitted_mandala = fitted_mandala
  )
  r_fit <- cor(comp_df$fitted_spats, comp_df$fitted_mandala)

  plot(
    comp_df$fitted_spats, comp_df$fitted_mandala,
    xlab = "SpATS fitted values",
    ylab = "Mandala fitted values",
    main = "SpATS vs Mandala fitted values"
  )
  abline(0, 1, col = "red", lty = 2)
  legend(
    "topleft",
    legend = sprintf("r = %.4f", r_fit),
    bty = "n"
  )

  data.frame(fitted_correlation = r_fit)
} else {
  data.frame(note = "Install SpATS to run this comparison.")
}

#>   fitted_correlation
#> 1          0.9981549

Selected Literature

Spatial adjustment is a central component of modern field-trial analysis because fertility gradients, management zones, and local field trends can inflate residual variation or bias genotype comparisons. Row-column random effects, separable AR1 residuals, empirical variograms, and two-dimensional spline smooths are common tools for diagnosing and modeling these trends.

Gilmour, A. R., Cullis, B. R. and 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.
Smith, A., Cullis, B. and Thompson, R. (2001). Analyzing variety by environment data using multiplicative mixed models and adjustments for spatial field trend. Biometrics, 57, 1138-1147.
Cullis, B. R., Smith, A. B. and Coombes, N. E. (2006). On the design of early generation variety trials with correlated data. Journal of Agricultural, Biological, and Environmental Statistics, 11, 381-393.
Eilers, P. H. C. and Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11, 89-121.
Rodríguez-Álvarez, M. X., Boer, M. P., van Eeuwijk, F. A. and Eilers, P. H. C. (2018). Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spatial Statistics, 23, 52-71.