# Linear Models with Correlated Errors

Now that we have a sense of how to model spatial variation, the next step is to incorporate that into a linear model. The starting point is the linear mixed model. In RCBD design, often the treatments are treated as fixed and the block effect as random.

$$Y_ij = \mu + \alpha_i + \beta_j + \epsilon_{ij}$$

$Y_ij$ is the independent variable

$\mu$ is the overall mean

$\alpha_i$ is the effect due to the $i^{th}$ treatment

$\beta_j$ is the effect due to the $j^{th}$ block

$\epsilon_{ij}$ are the error terms distributed as $N ~\sim (0,\sigma)$

Here is an expanded version of the last term:

$$ \epsilon_{ij} ~\sim N \Bigg( 0,
\left[ { \begin{array}{ccc} \sigma & \cdots & 0 \\

\vdots & \ddots & \vdots \\

0 & \cdots & \sigma \end{array} } \right] \Bigg) $$

This is a mathematically representation of **iid**, independent and identically distributed, an assumption of linear models. When there is spatial autocorrelation, observations closer to one another are correlated, so the off-diagonals in the variance-covariance matrix are not zero.

Spatial models seek to mathematically model this covariance so it is properly accounted for during hypothesis testing and prediction.

### Code for this section

*The following scripts build upon work done in previous section(s).*

## R

```
library(emmeans); library()
# (nlme and gstat should already be loaded)
library(spaMM) # for running `corMatern()`
# standard linear model
nin_lme <- lme(yield ~ gen, random = ~1|rep,
data = Nin,
na.action = na.exclude)
# extract the esimated marginal means for variety
preds_lme <- as.data.frame(emmeans(nin_lme, "gen"))
# use information from the variogram fitting for intialising the parameters
nugget <- Nin_variofit_gau$psill[1]
range <- Nin_variofit_gau$range[2]
sill <- sum(Nin_variofit_gau$psill)
nugget.effect <- nugget/sill
# initalise the covariance structure (from the nlme package)
cor.gaus <- corSpatial(value = c(range, nugget.effect),
form = ~ row.length + col.width,
nugget = T, fixed = F,
type = "gaussian",
metric = "euclidean")
# update the rcbd model
nin_gaus <- update(nin_lme, corr = cor.gaus)
# extract predictions for 'gen'
preds_gaus <- as.data.frame(emmeans(nin_gaus, "gen")
# a similar procedure can be follow for other models
# but we are going to take a shortcut and not specify the parameters
# exponential
cor.exp <- corSpatial(form = ~ row.length + col.width,
nugget = T, fixed = F)
nin_exp <- update(nin_lme, corr = cor.exp)
preds_exp <- as.data.frame(emmeans(nin_exp, "gen"))
# Matern structure
cor.mat <- corMatern(form = ~ row.length + col.width,
nugget = T, fixed = F)
nin_matern <- update(nin_lme, corr = cor.mat)
preds_mat <- as.data.frame(emmeans(nin_matern, "gen")
```

## SAS

```
proc mixed data=alliance ;
class entry rep;
model yield = entry ;
random rep;
lsmeans entry/cl;
ods output LSMeans=NIN_RCBD_means;
title1 'NIN data: RCBD';
run;
proc mixed data=alliance maxiter=150;
class entry;
model yield = entry /ddfm=kr;
repeated/subject=intercept type=sp(gau) (Row Col) local;
parms (11) (22) (19);
lsmeans entry/cl;
ods output LSMeans=NIN_Spatial_means;
title1 'NIN data: Gaussian Spatial Adjustment';
run;
```