In oharar/LatentINLA: Package to Write GLLVMs in INLA

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

This is an analysis of plant data from central Norway. The data consist of counts of 90 plant species from quadrats throughout Levanger, a town on the Trondheim Fjord. Each site also has 32 environmental variables.

First, set up and get the data:

library(LatentINLA)
data("Levanger")
# remove rare species & relatively empty site
useSpecies <- colMeans(Levanger$species>0)>0.5

Y <- Levanger$species[, useSpecies]

We will use 1 latent variable. First we fit the model (of course), using a normal distribution, which does not really fit the structure of the data.

#Lev <- FitGLLVM(Y=Y, Family="gaussian", nLVs = 1, 
Lev <- FitGLLVM(Y=Y, Family="poisson", nLVs = 1, 
                RowEffPriorsd = 1, ColEffPriorsd = 1,
                RowEff = "random", ColEff = "random", ColScorePriorsd=1)

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(Lev$rowterm) <- rownames(Lev$comm)
par(mar=c(4.1,4,1,1))
plotterm(Lev$rowterm, xlab="Site Mean Effect", ylab="")
mtext("Site", 2, line =5)

And then the species:

par(mar=c(4.1,6,1,1))
rownames(Lev$colterm) <- colnames(Levanger$species)[useSpecies] # might not be right
plotterm(Lev$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 <- Levanger$site[, c("management", "soil_depth")]
LevCov <- FitGLLVM(Y=Y, 
                      X=X, 
#                   W=forminifera$traits[useSpecies, ],
#                   Family="gaussian", nLVs = 1, 
                   Family="poisson", nLVs = 1, 
                RowEffPriorsd = 1, ColEffPriorsd = 1,
                RowEff = "random", ColEff = "random", ColScorePriorsd=1)

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

par(mar=c(2.1,6,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))
  rownames(ToPlot) <- gsub(".", "\n", rownames(ToPlot), fixed = TRUE)
  plotterm(ToPlot, xlab="", ylab="")
  mtext(env, 3, line=1, cex=1.3)
  abline(v=0)
}, fixed=LevCov$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(LevCov$colterm) <- gsub("column", "", rownames(LevCov$colterm))
plotterm(LevCov$colterm, xlab="Species Mean Effect", ylab="")
mtext("Species", 2, line =5)