In oharar/LatentINLA: Package to Write GLLVMs in INLA

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Forminifera Data GLLVM

This is an analysis the much-used spider data. Everyone else uses is, so we will too.

First, set up and get the data:

library(LatentINLA)
data("forminifera")
# remove rare species & relatively empty site
useSite <- 1:29
useSpecies <- 1:22
Y <- forminifera$comm[useSite, useSpecies]

We will use 2 latent variables. First we fit the model (of course). We will ignore the warnings this produces.

formin <- FitGLLVM(Y=Y, Family="poisson", nLVs = 2, 
                   RowEff = "random", ColEff = "random", ColScorePriorsd=1,
                   control.inla = list(int.strategy = "eb"))

We can admire our wonderful ordination plot, and wonder what sp22 is.

biplot(formin)

We have also estimated our site and species main effects. First the sites. Using a random effect helps the estimation.

plotterm <- function(t, ...) {
  plot(t$mean, 1:nrow(t), xlim=c(min(t$'0.025quant'), max(t$'0.975quant')),
       yaxt="n", ...)
  segments(t$'0.025quant', 1:nrow(t), t$'0.975quant', 1:nrow(t))
  axis(2, rownames(t), at=1:nrow(t), las=1)
}
rownames(formin$rowterm) <- rownames(formin$comm)
par(mar=c(4.1,4,1,1))
plotterm(formin$rowterm, xlab="Site Mean Effect", ylab="")
mtext("Site", 2, line =5)

And then the species:

par(mar=c(4.1,6,1,1))
rownames(formin$colterm) <- colnames(forminifera$comm)[1:22] # might not be right
plotterm(formin$colterm, xlab="Species Mean Effect", ylab="")
mtext("Species", 2, line =5)

From this we can conclude that some species are rarer than others.

Adding covariates

Now we can add some environmental covariates.

X <- apply(forminifera$envir[useSite, 1:3], 2, scale)
forminCov <- FitGLLVM(Y=Y, 
                      X=X, 
#                   W=forminifera$traits[useSpecies, ],
                   Family="poisson", nLVs = 2, 
                   RowEff = "random", ColEff = "fixed", ColEffPriorsd=1)

Whic also give us a boplot:

biplot(forminCov)

But we can also see the effects of the covariates on each species:

par(mar=c(2.1,5,2.4,1), mfrow=c(1,3), oma=c(2,2,0,0))
sapply(colnames(X), function(env, fixed) {
  Use <- rownames(fixed)!="Intercept" & grepl(env, rownames(fixed))
  ToPlot <- fixed[Use,]
  rownames(ToPlot) <- gsub(".*:column", "", rownames(ToPlot))
  plotterm(ToPlot, xlab="", ylab="")
  mtext(env, 3, line=1, cex=1.3)
  abline(v=0)
}, fixed=forminCov$fixed)
mtext("Species", 2, line = 0.1, outer=TRUE)
mtext("Coefficient", 1, line = 1, outer=TRUE)

And then the species main effects, which suggest that the mean abundances are the same.

par(mar=c(4.1,6,1,1))
rownames(forminCov$colterm) <- gsub("column", "", rownames(forminCov$colterm))
plotterm(forminCov$colterm, xlab="Species Mean Effect", ylab="")
mtext("Species", 2, line =5)

Overall, sp22 is weird.

A Constrained Model (if it works)

Now we can add some environmental covariates.

X <- apply(forminifera$envir[useSite, 1:3], 2, scale)
forminCov <- FitConstrainedGLLVM(Y=Y, 
                      X=X, 
#                   W=forminifera$traits[useSpecies, ],
                   Family="poisson", nLVs = 2)

#                   RowEff = "random", ColEff = "fixed", ColEffPriorsd=1)