AICc vs. Model Diagnostics

[KemenesT]

Hi,

I'm doing an analysis on a community of bacteria at the Phylum level. I find that if I include 2 latent variables my AICc is lower than if I include 0. But looking at the diagnostics plots it looks to me like the model with 0 latent variables is better. The results are otherwise very similar. This made me wonder, if there are 0 latent variables, is that equivalent to doing separate regressions for all the phyla? My understanding was that the latent variables are the link that allows the coefficient estimation for each Phylum to influence each other.

Do you have any advice on the best approach here? Why does including the latent variables make such a difference to the outcomes? Any advice is very much appreciated :)

The groups have the following distributions (log scale):

These are the models:

Formula and model AICc's

formula_obj <- formula(~ fDepth + PO4s + NH4s + NO3s + NO2s)

M1.beta df = 36 AICc = -1512.3421
M2.beta df = 45 AICc = -1716.0608

Model 1

M1.beta <- gllvm(y = community_beta/100,
X = Environment,
studyDesign = StudyDesign,
row.eff = ~ (1 | fSiteCode),
formula = formula_obj,
num.lv = 0,
family = "beta",
method = "LA",
control.start = list(n.init = 5,
jitter.var = 0.1))

Model 2

M2.beta <- gllvm(y = community_beta/100,
X = Environment,
studyDesign = StudyDesign,
row.eff = ~ (1 | fSiteCode),
formula = formula_obj,
num.lv = 2,
family = "beta",
method = "LA",
control.start = list(n.init = 5,
jitter.var = 0.1))

Thank you!

[AlainZuur]

I'm sure Bert can give you better advice...but to answer part of your question...you have a random effect that is shared between the species. So...when num.lv = 0, they still share a random effect...hence it is not the same as 5 independent GLMMs. Without the random effect, it would be. And I think that Bert will tell you that 2 LVs may be an overkill for your data....try 1. You can see in your residual plot that two species have very small residuals...that is a bit suspicious..and I guess you get a perfect fit for those species?

Anyway...my 2 cents.

Alain

[KemenesT]

Thanks Alain!

[BertvanderVeen]

I largely agree with Alain his answer.

In principle you are right: the latent variables here are what "connect" the species, and Alain is right that the random site effect also does that.

Two latent variables with a few species is probably overkill. That is also what looks like is happening from your residuals vs. species plot where all the residuals for the second and fourth species lay on a single point. How many observations per species do you have, and how many covariates are you including in the model? From the number of parameters it looks like around 5-6 covariates.

If you really want to create an ordination plot (which sometimes you just do, and that is OK with me), looking at information criteria might not be the right way to go. AIC(c) helps you to find a model that might predict well, but it will not help you (much) in finding a model for inference. You will have to find a way to improve the residual plots though, by adjusting your model.

[KemenesT]

Thanks so much for your help both Alain and Bert!

I see, no I'm not too concerned with the ordination. I am trying rather to examine the different phyla in relation to the covariates, but I am trying to avoid doing multiple separate regressions because all these bacteria inhabit the same mats. The random effects linking them together addresses my main concern there then.
I have 91 samples and all the species occur in all of them, so 91 observations per species. In some samples they do occur in very low abundances, just not 0. Then I have 5 covariates and a site random effect with 13 sites. Per site I have about 7 samples, thought they're not perfectly balanced.

So it seems to me that M1 with no latent variables is the way to go to answer my question, and most robust? (I did try 1 LV but it still affects the fit poorly)

[BertvanderVeen]

Yes, it does seem like you do not need the LVs here. As long as you have refitted the model a few times, and with LVs it always fits worse, than omitting the LVs should be fine :)

[gerverska]

Thanks @KemenesT for starting this discussion, and thanks @BertvanderVeen and @AlainZuur for your thoughts and answers. I'm curious what @KemenesT and others might do in this situation to obtain residual covariances and residual correlations. I'm assuming that it follows something along the line of this solution, but it seems this will also require the application of the appropriate, distribution-specific adjustments used in getResidCov() / getResidCor().

The solution I've linked to seems very lme4-specific, so what would the equivalent approach be for gllvm? I'm currently going through the getResidCov() code to see if I can adapt something for random effects-only situations.