Spatial random variable with 9,738 coordinates causes R to crash

[devilevie]

Hi all,

I'm trying to fit a HMSC model for the presence-absence of 23 mammal species across 9,738 points (random points across species ranges). I managed to get the model working with point ID as a random effect, but when I try to change it to spatial coordinates, it crashes. Here's some of the relevant code:

studyDesign = data.frame(Coords = as.factor(1:nrow(Y)))

#SPATIAL EFFECTS
xycoords = XData[, 2:3]
xycoords = as.matrix(xycoords)
rownames(xycoords) = 1:9738
rL.spatial = HmscRandomLevel(sData = xycoords)

m = Hmsc(Y=Y, XData = XData, XFormula = XFormula,
phyloTree = taxonomicTree,
distr="probit",
studyDesign = studyDesign, ranLevels = list(Coords = rL.spatial))

models = list(m)
names(models) = c("presence-absence model")

Runs fine until here, but then when I run the following to test the model it causes R to crash:

for(i in 1:length(models)){
print(i)
sampleMcmc(models[[i]],samples=2)
}

Any idea what I'm doing wrong? Thanks for your help.

[gtikhonov]

You have absolutely too many spatial locations in your model for vanilla Gaussian Process (GP) -priored latent factors to work well. Try formulations of HmscRandomLevel(...) that rely on approximate GP structures. For instance,

rL.spatial = HmscRandomLevel(sData = xycoords, sMethod ="NNGP")

See the help page of HmscRandomLevel(...) for more details or this paper (https://doi.org/10.1002/ecy.2929) for statistical presentation of implemented methods.

[devilevie]

Thanks, that's really helpful. I'll try the new code and have a look at the paper, hopefully I can get it working.