This Presentation:

https://github.com/IdahoAgStats/lattice-spatial-analysis-talk

A longer tutorial:

https://idahoagstats.github.io/guide-to-field-trial-spatial-analysis

• Perceived lack of need
• Unsure of benefits
• No training in the topic/intimidated by the statistical methodology
• Limited time to devote to statistical analysis
• Unclear what would happen to blocking if spatial stats are used
• very few resources for easy implementation

This is not how this works. Blocking is compatible with spatial analysis and recommended for most (all?) field trials.

• Experimental treatments
• fully crossed effects
• Blocking scheme along the expected direction of field variation

A typical linear model:

\(Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}\)

Response = trial mean + treatment effect + block effect + leftover error

We Assume:

1. The error terms, or residuals, are independent of another with a shared distribution:

\[\epsilon_i \sim N(0,\sigma_e)\]

1. Each block captures variation unique to that block and there is no other variation related to spatial position of the experimental plots.

• 36 soft white winter wheat cultivars
• 4 blocks
• 2 missing data points
• the linear model:

\(Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}\)

library(nlme)
lm1 <- lme(yield ~ cultivar, random = ~ 1|block, data = mydata, na.action = na.exclude)

plot(lm1)

\(H_0\): There is no spatial autocorrelation
\(H_a:\) There is spatial autocorrelation!

This uses a simple weighting matrix that weights all neighbors that share a plot border (the chess-based “rook” formation) equally.

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mydata$residuals 
## weights: weights 
## omitted: 88, 97 
## number of simulations + 1: 1000 
## 
## statistic = 0.15869, observed rank = 999, p-value = 0.001
## alternative hypothesis: greater

Areal data = finite region divided into discrete sub-regions (plots) with aggregated outcomes

Options:

1. model row and column trends
   • good for known gradients (hill slope, salinity patterns)
2. assume plots close together are more similar than plots far apart. The errors terms can be modelled based on proximity, but there is no trial-wide trend
   • autoregressive models (AR)
   • models utilizing “gaussian random fields” for continuously varying data (e.g. point data)
   • Smoothing splines
   • nearest neighbor

\[Y_{ij} = \mu + A_i + \epsilon_{ij}\] \[\epsilon_i \sim N(0,\sigma)\]

If N = 4:

\[e_i ~\sim N \Bigg( 0, \left[ {\begin{array}{ccc} \sigma^2 & 0 & 0 & 0\\ 0 & \sigma^2 & 0 & 0\\ 0 & 0 & \sigma^2 & 0\\ 0 & 0 & 0 & \sigma^2 \end{array} } \right] \Bigg) \]

The variance-covariance matrix indicates a shared variance and all off-diagonals are zero, that is, the errors are uncorrelated.

Same linear model: \[Y_{ij} = \mu + A_i + \epsilon_{ij}\]

Different variance structure:

\[e_i ~\sim N \Bigg( 0, = \sigma^2 \left[ {\begin{array}{cc} 1 & \rho & \rho^2 & \rho^3 \\ \rho & 1 & \rho & \rho^2 \\ \rho^2 & \rho & 1 & \rho \\ \rho^3 & \rho^2 & \rho & 1 \\ \end{array} } \right] \Bigg) \]

• \(\rho\) is a correlation parameter ranging from -1 to 1 where 0 is no correlation and values approaching 1 indicate spatial correlation.
• The “one” in AR1 means that only the next most adjacent point is considered. There can be AR2, AR3, …, ARn models.

• From a statistical standpoint, it’s one of the more intuitive models
• The implementation in R is a little shaky
  • several packages, all hard to use and incompatible with other R packages
• It is implemented in SAS
• Some proprietary software implements this (AsREML), others do not (Agrobase)

A measure of spatial correlation based on all pairwise correlations in a data set, binned by distance apart:

\(\gamma^2(h) = \frac{1}{2} Var[Z(s+h)-Z(s)]\)
\(Z(s)\) = observed data at point \(s\).
\(Z(s)\) = observed data at another point \(h\) distance from point \(s\).

For a data set with \(N\) observation, there are this many pairwise points:

\(\frac{N(N-1)}{2}\)

This uses semivariance to mathematically relate spatial correlations with distance

range = distance up to which is there is spatial correlation sill = uncorrelated variance of the variable of interest nugget = measurement error, or short-distance spatial variance and other unaccounted for variance

• There are many difference mathematical models for explaining semivariance:
  • exponential, Gaussian, Matérn, spherical, …
    
• It is usually used for kriging, or prediction of a new point through spatial interpolation
  
• It can also be used in a linear model where local observations are used to predict a data point in addition to treatment effects
  
• Bonus: R and SAS are really good at this!

Copy data into new object so we can assign it a new class (and remove missing data):

library(gstat); library(sp); library(dplyr)
mydata_sp <- mydata %>% filter(!is.na(yield))

Establish coordinates for data set to make it an sp object (“spatial points”):

coordinates(mydata_sp) <- ~ row + range

Set the maximum distance for looking at pairwise correlations:

max_dist <- 0.5*max(dist(coordinates(mydata_sp)))

Calculate a sample variogram:

semivar <- variogram(yield ~ block + cultivar, data = mydata_sp,
                        cutoff = max_dist, width = max_dist/12)
nugget_start <- min(semivar$gamma)

The empirical variogram:

plot(semivar)

Set up models for fitting variograms:

vgm1 <- vgm(model = "Exp", nugget = nugget_start) # exponential
vgm2 <- vgm(model = "Sph", nugget = nugget_start) # spherical
vgm3 <- vgm(model = "Gau", nugget = nugget_start) # Gaussian

Fit the variogram models to the data:

variofit1 <- fit.variogram(semivar, vgm1)
variofit2 <- fit.variogram(semivar, vgm2)
variofit3 <- fit.variogram(semivar, vgm3)

Look at the error terms to see which model is the best at minimizing error.

## [1] "exponential: 26857.3"

## [1] "spherical: 26058.3"

## [1] "Gaussian: 41861.0"

The spherical model is the best at minimizing error.

plot(semivar, variofit2, main = "Spherical model")

Extract the nugget and sill information from the spherical variogram:

nugget <- variofit2$psill[1] 
range <- variofit2$range[2] 
sill <- sum(variofit2$psill) 
nugget.effect <- nugget/sill  # the nugget/sill ratio

Build a correlation structure in nlme:

cor.sph <- corSpatial(value = c(range, nugget.effect), 
                  form = ~ row + range, 
                  nugget = T, fixed = F,
                  type = "spherical", 
                  metric = "euclidean")

Update the Model:

lm_sph <- update(lm1, corr = cor.sph)

logLik(lm1)

## 'log Lik.' -489.0572 (df=38)

logLik(lm_sph)

## 'log Lik.' -445.4782 (df=40)

anova(lm1)[2,]

##          numDF denDF F-value p-value
## cultivar    35   103  1.6411   0.029

anova(lm_sph)[2,]

##          numDF denDF  F-value p-value
## cultivar    35   103 2.054749  0.0028

library(emmeans)
lme_preds <- as.data.frame(emmeans(lm1, "cultivar")) %>% mutate(model = "mixed model")
sph_preds <- as.data.frame(emmeans(lm_sph, "cultivar")) %>% 
  mutate(model = "mixed model + spatial")
preds <- rbind(lme_preds, sph_preds)

Highest yielding wheat: ‘Stephens’ (released in 1977)

• When models omit blocking, the predictions may be unchanged or they may worsen. This varies by the agronomic field, but in general, blocking a field trial and including block in the statistical model improves your experimental power and controls experimental error.
• There is no single spatial model that fits all
• However, using any spatial model is usually better than none at all
• When you use spatial covariates, your estimates are better and more precise. This really does help you!

• Track row and range information in your trial data set.
• Look at the tutorial! (we will also add SAS code)
• Try out a few models and see how it impacts your results.

Statistical consulting to support the College of Agriculture and Life Sciences.

Bill Price, Director, bprice@uidaho.edu, AgSci307

Julia Piaskowski, jpiaskowski@uidaho.edu, AgSci 305