R2GLMM sensu Nakagawa et al. 2017

[gerverska]

Hi everyone,

Related to #189, I have been having difficulties getting stable output from models including LVs, but models that include random row effects seem more easily reproduced. While I'm still not entire sure how to go about obtaining output similar to getResidCov() + getResidCor() that only work with LVs, I was also thinking it'd be helpful to see if I could at least produce an R2 for random row effects + fixed effect models.

Digging into the supplementary code of Nakagawa and Schielzeth 2013 and Nakagawa et al. 2017, I think I have adapted a solution for gllvm, but I'm not sure if there are other considerations for the multivariate context or whether I'm accessing the appropriate gllvm parameters/outputs.

# Adapted from Nakagawa et al. 2017, Appendix S6 ####
# Original code has been retained to help cross-reference my changes (below the line of original code)

# Load data ####
data('spider')

# Fit null and full models ####
model.null <- gllvm(spider$abund,
                    row.eff = 'random',
                    family = 'negative.binomial',
                    num.lv = 0,
                    gradient.check = T,
                    control = list(reltol = 1e-16),
                    control.start = list(n.init = 5))

model.full <- gllvm(spider$abund,
                    spider$x[1:2],
                    row.eff = 'random',
                    family = 'negative.binomial',
                    num.lv = 0,
                    gradient.check = T,
                    control = list(reltol = 1e-16),
                    control.start = list(n.init = 5))

# Calculation of the variance in fitted values ####
# VarF <- var(as.vector(model.matrix(parmodGLMERf) %*% fixef(parmodGLMERf)))
VarF <- var(as.vector(model.full$params$Xcoef %*% t(model.full$X.design)))

# lambda <- as.numeric(exp(fixef(parmodGLMERr) +
#                              0.5 * (as.numeric(VarCorr(parmodGLMERr)$Population) +
#                                         as.numeric(VarCorr(parmodGLMERr)$Container))))
lambda <- as.numeric(exp(model.null$params$sigma^2))

# Full model ####
# thetaF <- getME(parmodGLMERf, "glmer.nb.theta")
thetaF <- model.full$params$inv.phi

VarOdF <- 1/lambda + 1/thetaF # the delta method
VarOlF <- log(1 + (1/lambda) + (1/thetaF)) # log-normal approximation
VarOtF <- trigamma((1/lambda + 1/thetaF)^(-1)) # trigamma function
# comparing the three ####
c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)

# R2[GLMM(m)] - marginal R2[GLMM] ####
# R2glmmM <- VarF/(VarF + sum(as.numeric(VarCorr(parmodGLMERf))) + VarOtF)
R2glmmM <- VarF/(VarF + sum(as.numeric(model.full$params$sigma^2)) + sum(VarOlF))

Does this seem right? It seems like it can generate pretty high R2 for the spider data. Also, it seems like the log-normal approximation creates slightly higher R2.

[BertvanderVeen]

If you previously installed a development version of the repository - try reinstalling. There was a but in the repeated fitting procedure that messed things up a bit. If that is not the cause of the instability you are witnessing, it would be helpful to get some more details.

I am not quite sure what you expect me to comment on for your code: at least half of the code is for GLMMs fitted with lme4/glmmTMB, and thus not applicable to gllvm. Or am I missing something here?

[gerverska]

Hi Bert. The dev version I have is about a week or two old, but I can try re-downloading. The main issue I've encountered with LVs in some cases is that I would fit a model with 1-2 LVs with n.init = 5 up to n.init = 9, sometimes with jitter.var = 0.1, and I would see that AICc/BIC would jump around by 100 units or so and the ordiplot() would change run to run.

I would always set gradient.check = T, and most of the fits I would make run to run seemed to converge without warnings. I know you can set seed once you have refit the model a few times and achieved similar solutions with each refit, but I've found that row.eff = 'random' seemed to be better at achieving this repeatability after multiple refits than num.lv > 0.

When I was fitting 1-2 LVs in my null and full models and comparing their residual covariances (as in the first two gllvm vignettes) to obtain an R2, I would find that the R2 would jump around quite a bit and sometimes go negative, even if both null and full models converged without large gradients. Given the stability I think I'm seeing with row.eff = 'random, this motivated the idea to figure out a way to get an R2 for row effects only models.

Apologies for not explaining the code arrangement very well: the hashed out/commented lines of code do contain the original lme4 code to show how I was adapting the concept for gllvm output, shown in the un-hashed line below it. Sorry if this wasn't a helpful arrangement, but I was hoping it'd let you know that I wasn't coding this out of thin air! The lme4 version and the gllvm version obviously differ in the univariate/multivariate sense, but they also differ in how certain values are extracted from model objects, so I wasn't sure if I was proceeding correctly.

Main differences (lme4 vs gllvm):

• VarCorr(parmodGLMERr)$Population vs model.null$params$sigma^2: the lme4 version outputs the variance associated with the random effect when accessed like this, while the standard deviation for gllvm random effects (thus the squaring) is kept in params$sigma, correct?
• getME(parmodGLMERf, "glmer.nb.theta") vs model.full$params$inv.phi: this one was a little more difficult for me to track down the equivalent term, but I think these two should be the same. Because each taxon response has it's own inv.phi, the calculations of VarOdF, VarOlF, and VarOtF are output as a vectors, and one of these last three (I used the lognormal approx VarOlF) is summed in the final calculation of R2GLMM. I'm assuming this is the appropriate way to handle the multiple inv.phis?

Here's a link to the Nakagawa et al. 2017 appendix I'm referring to, for reference.

[gerverska]

And here's the code, removing the lme4 references:

# Adapted from Nakagawa et al. 2017, Appendix S6 ####

# Load data ####
data('spider')

# Fit null and full models ####
model.null <- gllvm(spider$abund,
                    row.eff = 'random',
                    family = 'negative.binomial',
                    num.lv = 0,
                    gradient.check = T,
                    control = list(reltol = 1e-16),
                    control.start = list(n.init = 5))

model.full <- gllvm(spider$abund,
                    spider$x[1:2],
                    row.eff = 'random',
                    family = 'negative.binomial',
                    num.lv = 0,
                    gradient.check = T,
                    control = list(reltol = 1e-16),
                    control.start = list(n.init = 5))

# Calculation of the variance in fitted values ####
VarF <- var(as.vector(model.full$params$Xcoef %*% t(model.full$X.design)))

# Getting the observation-level variance Null model ####
lambda <- as.numeric(exp(model.null$params$sigma^2))

# Full model ####
thetaF <- model.full$params$inv.phi

VarOdF <- 1/lambda + 1/thetaF # the delta method
VarOlF <- log(1 + (1/lambda) + (1/thetaF)) # log-normal approximation
VarOtF <- trigamma((1/lambda + 1/thetaF)^(-1)) # trigamma function

# Comparing the three ####
c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)

# R2[GLMM(m)] - marginal R2[GLMM] ####
R2glmmM <- VarF/(VarF + sum(as.numeric(model.full$params$sigma^2)) + sum(VarOlF))