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
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[]))) : 'names' attribute  must be the same length as the vector 
Any advice on how I can properly encode this? (I found a related, but unanswered question here).