Environmental correlation in gllvm

[Leitemfa]

Hi,

I see the function for extracting the residual correlation from the model. Is it also possible to get the shared environmental response, similarly as it is done in the boral package?

[fhui28]

Please see the getResidualCor.gllvm and getResidualCov.gllvm functions for the former. For latter, as far as I know you will need to manually construct the estimated linear predictors due to the X variables, before calculating applying cov and then cov2cor

[Leitemfa]

Hi @fhui28 Thanks for you help, but I'm still struggling with that. How do I manually construct this linear predictor? Can't I take it from the model$params$Xcoef?

Would you mind to walk me through this construction?

[fhui28]

Here is a worked example based on the first example in the gllvm package.

library(gllvm)
library(tidyverse)

data(microbialdata)
X <- microbialdata$Xenv
y <- microbialdata$Y[, order(colMeans(microbialdata$Y > 0), decreasing = TRUE)[21:40]]
fit <- gllvm(y, X, formula = ~ pH + Phosp, family = poisson())
fit$logL

# Residual correlation
getResidualCor.gllvm(fit)

# Environmental correlation
modelmat <- model.matrix(fit$formula, data = X) %>% 
   as.matrix %>% 
   {.[,-1]} # Remove intercept
linpred <- tcrossprod(modelmat, fit$params$Xcoef)
envircor <- cov2cor(cov(linpred))

There might well be an easier way to do this [not up to date with the updates!], but I did not spot it.

Good luck!

[Leitemfa]

Thank you very much @fhui28 , It works amazingly! I just have few final questions now: is it possible to check the effect of one environmental variable at time? For the example you show, does it make sense to get the subset from the model.matrix or should I run a new model with the variables separated?

modelmat <- model.matrix(fit$formula, data = X) %>% 
   as.matrix %>% 
   {.[,-1]} # Remove intercept
linpred.ph <- tcrossprod(modelmat[,1], fit$params$Xcoef[,1])
envircor.ph <- cov2cor(cov(linpred.ph))
linpred.phosp <- tcrossprod(modelmat[,1], fit$params$Xcoef[,1])
envircor.phosp <- cov2cor(cov(linpred.phosp))

Another question is how do I check the confidence interval of these correlation values?

[fhui28]

is it possible to check the effect of one environmental variable at time? For the example you show, does it make sense to get the subset from the model.matrix or should I run a new model with the variables separated?

In principle you can, but it's kind of meaningless to do as the correlation between any two species will have then directly and simply proportional to the difference between the two estimated coefficients of that covariate for those two species.

Another question is how do I check the confidence interval of these correlation values?

Sounds like something for simulation @jebyrnes?

[BertvanderVeen]

I will leave a note here now that with the addition of random predictor slopes with num.RR it would be possible to extract species correlations due to the environment more explicitly (and with randomB="P" they would differ for predictors). Though functionality for this is not currently available.

It would not be too difficult to write something along the lines of getResidualEnvCov.gllvm and getResidualEnvCor.gllvm, but I do not have time to do that in the next month or so.