0
$\begingroup$

I'm working on a modified example from a glmmTMB vignette (found here) 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):

#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":

#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:

#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:

#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

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

#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

This is a far cry from the from predictions made by the example model in the vignette, which are much closer to the original volcano dataset:

enter image description here

$\endgroup$
0

2 Answers 2

1
$\begingroup$

I figured it out:

  1. ranef.glmmTMB was re-ordering the factor levels, and includes the intercept, which I didn't realize
  2. I mixed up rows/columns in predict_cols
#How do models do at estimating spatial random effect?
checkRE1 <- data.frame(s,re=ranef(mod1)[[1]][[1]][,1])
#For some reason, ranef.glmmTMB includes the intercept
checkRE2 <- data.frame(s[with(s,order(y,x)),],re=unlist(ranef(mod2,condVar = FALSE)[[1]][[1]])+fixef(mod2)$cond[[1]])

par(mfrow=c(2,1)) #Spatial RE estimates are closer to actual
with(checkRE1,plot(z,re,xlab='Actual',ylab='Estimated',main='Non-spatial RE')); abline(0,1)
with(checkRE2,plot(z,re,xlab='Actual',ylab='Estimated',main='Spatial RE')); abline(0,1)
par(mfrow=c(1,1))

Non-spatial vs Spatial Random Effects

#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(1:nrow(volcano),i))) 
  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)

Compare results to input

$\endgroup$
0
$\begingroup$

I'm not an expert with glmmTMB and your error structure is not initially intuitive to me... but here is a recommendation made to me when I was experimenting with error structure in mixed models: instead of only comparing those two models you've provided, try creating even more models with with all possible and justifiable error structures, and compare their performance using this helpful package. Again I am not an expert and this may not be what you need... but it just might help you progress.

$\endgroup$
1
  • $\begingroup$ I'm not really interested in picking "optimal error structures", since I generated the data and I know what the error structure is. I'm interested in why the spatial random effect performed worse than the non-spatial random effect. $\endgroup$ Jan 5, 2021 at 21:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.