# How to specify a model with spatial correlation AND crossed random intercepts in R (nlme)

I want to model repeated measurements (let's say yearly) from multiple sites. Let's say after modeling these data with random intercepts for site and year we find there remains spatial correlation in residuals between sites. Here's an example dataset...sorry no actual structure in the data, which I hope is fine because I'm just trying to get the nlme syntax correct.

#set number of sites and years
n.sites=10
n.years=25

#populate dataframe
n.obs<- n.sites*n.years
df<- data.frame(
site=rep(letters[1:n.sites], each=n.years),
year=paste("year", rep(1:n.years, n.sites)),
y=rnorm(n.obs),
coord.x=rep(runif(n.sites), each=n.years),
coord.y=rep(runif(n.sites), each=n.years)
)


To my knowledge, lme() in the nlme package is my best option (outside of moving to WinBUGS or similar).

There is a good worked example of how to model spatially correlated data in nlme here, however because my sites all have the same lat-lon in every year I will get the error below if I first try a simpler model without crossed random intercepts.

lme(y~1, random= ~1|site,
correlation = corGaus(form = ~ coord.x + coord.y|site),
data=df)
Error in getCovariate.corSpatial(object, data = data) :
cannot have zero distances in "corSpatial"


Which I believe is because after grouping observations by site, there are still repeated lat lon (i.e., "zero distances") for each year in a given site. So I clearly need to include the crossed random intercepts as the grouping factor in the correlation= argument.

Specifying crossed random effects in lme() is tricky (at least compared to lme4), but here is how I do it:

df\$dummy<- 1
lme(y~1, random=list(dummy=pdBlocked(list(pdIdent(~year-1),
pdIdent(~site-1)))),
data=df)


However, when I try to copy & paste this to the correlation= argument I get an error which believe is telling me my syntax is not correct.

lme(y~1, random=list(dummy=pdBlocked(list(pdIdent(~year-1),
pdIdent(~site-1)))),
correlation = corGaus(form = ~ coord.x + coord.y|list(dummy=pdBlocked(list(pdIdent(~year-1),                                                                              pdIdent(~site-1))))),
data=df)
Error in names(val) <- unlist(lapply(val, function(el) deparse(el[[2]]))) :
'names' attribute [2] must be the same length as the vector [1]


Any advice on how I can properly encode this? (I found a related, but unanswered question here).

• A simple hack would be to make the different years an arbitrarily far distance apart (such adding some multiple of the year to the original XY coordinates). I believe you could use this and pass a maximum distance to the value option that is smaller than the induced distance for the fake coordinates to achieve what you want. – Andy W Mar 6 '18 at 23:37