I'm working on a modified example from a `glmmTMB` vignette (found [here](https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html)) using spatial covariance structures. I'm trying to show the difference between "standard" random intercepts and "spatial" random intercepts models, similar to models that ecologists or geologists might be interested in (e.g. independent and dependent variable are measured at many sites across a landscape).

My main questions are:

1. Am I specifying the `glmmTMB` model correctly?
2. If so, what is causing the poor performance of the spatial model compared to the non-spatial model? I would have thought that the spatial information would *improve* its performance

First, I generate a model matrix and coefficients for the fixed effects (single dependent variable **x**):

```r
#Load glmmTMB
library(glmmTMB)

set.seed(1234)

Nsites <- 100 #Number of sites
N <- 300 #Number of samples (3 per site)

x <- runif(N,-1,1) #Single dependent variable "x"
fixMM <- model.matrix(~x) #Model matrix for fixed effects
fixCoefs <- c(0,1) #Fixed effect coefficients
sites <- rep(factor(1:Nsites),length.out=N) #Assign to sites 

```

Then I use the `volcano` dataset (scaled) to get some spatially-correlated random intercepts for each "site":

```r
#Spatial random effect (from "volcano" elevation dataset)
volcano <- (volcano-mean(volcano))/sd(volcano) #Rescale elevation values for model
s <- data.frame(z = as.vector(volcano),
                x = as.vector(row(volcano)),
                y = as.vector(col(volcano)))
s <- s[sample(nrow(s),Nsites),] #Sample sites
s$siteName <- factor(1:Nsites) #Label sites

```

Then I generate the data and assemble into a dataframe, suitable for fitting in `glmmTMB`:

```r
#Generate data from "fixed" and site effect
fixedEff <- fixMM %*% fixCoefs #Fixed effect
siteEff <- model.matrix(~sites-1) %*% s$z #Site effect
yhat <- fixedEff + siteEff #Predicted value
y <- yhat + rnorm(N,0,0.5) #Data (Residual SD = 0.5)

#Add spatial coordinates as E and N
siteE <- model.matrix(~sites-1) %*% s$x 
siteN <- model.matrix(~sites-1) %\*% s$y 

#Set up spatial data for use in glmmTMB model
siteCoords <- numFactor(siteE,siteN) #Coordinates
siteGroup <- factor(rep(1,N)) #Group

#Assemble dataframe
dat <- data.frame(y,x,sites,siteE,siteN,siteCoords,siteGroup) 

```

Then I fit the two types of models, and see how they do at predicting the fixed and random effects:

```r
#Model using only site ID
mod1 <- glmmTMB(y~x+(1|sites),data=dat,family=gaussian())

#Model using site location
mod2 <- glmmTMB(y~x+exp(siteCoords + 0|siteGroup),data=dat,family=gaussian())

#How do models do at estimating fixed effects (int = 0, slope = 1)
data.frame(NonSpatial=confint(mod1)[1,],Spatial=confint(mod2)[1,]) #Intercept
data.frame(NonSpatial=confint(mod1)[2,],Spatial=confint(mod2)[2,]) #Slope

#How do models do at estimating spatial random effect?
checkRE <- data.frame(s,ranefM1=ranef(mod1)[[1]][[1]][,1],
                     ranefM2=unlist(ranef(mod2)[[1]][[1]]))

par(mfrow=c(2,1)) #Spatial model does a poor job
with(checkRE,plot(z,ranefM1,xlab='Actual',ylab='Estimated',main='Non-spatial RE')); abline(0,1)
with(checkRE,plot(z,ranefM2,xlab='Actual',ylab='Estimated',main='Spatial RE')); abline(0,1)
par(mfrow=c(1,1))

```

As shown below, the spatial model (`mod2`) does a terrible job at estimating the random intercepts for each site:

[![Random effects from mod1 and mod2](https://i.stack.imgur.com/pE48w.png)](https://i.stack.imgur.com/pE48w.png)

After extracting the results for the whole field, it is again clear that the spatial model is doing a bad job:

```r
#Can the spatial model reconstruct "volcano" image?
predict_col <- function(i) { #Function to construct 1 column at a time (more efficient)
  newdata <- data.frame( siteCoords = numFactor(expand.grid(x=i,y=1:nrow(volcano)))) #i=cols,j=rows
  newdata$siteGroup <- factor(rep(1,nrow(newdata)))
 newdata$x <- 0 #Marginalize across x
  predict(mod2, newdata=newdata, type="response", allow.new.levels=TRUE)
}

#Takes a few minutes to predict all responses 
pred <- sapply(1:ncol(volcano),predict_col)

#Plot predictions + errors
par(mfrow=c(2,2))

#Original data
image(volcano,main='Original')

#Data from sites
volcano.data <- array(NA, dim(volcano))
volcano.data[cbind(s$x, s$y)] <- s$z
image(volcano.data, main="Sampled data", useRaster=TRUE)

#Plot data from predictions
image(pred, main="Predicted data", useRaster=TRUE)

#How wrong are the predictions?
error.data <- abs(volcano-pred)
image(error.data, main="abs(Error)", useRaster=TRUE)

```

[![Prediction errors from original](https://i.stack.imgur.com/JOX9L.png)](https://i.stack.imgur.com/JOX9L.png)

This is a far cry from the from predictions made by the example model in the [vignette](https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html), which are much closer to the original `volcano` dataset:

[![enter image description here](https://i.stack.imgur.com/3hX2D.png)](https://i.stack.imgur.com/3hX2D.png)