I'm working with OTU count data (essentially counts of clustered sequence reads), and trying to determine differential abundance of counts between two groups. As of now, we're using a method defined in Romero, et al's paper titled "The composition and stability of the vaginal microbiota of normal pregnant women is different from that of non-pregnant women", whereby we fit the model using Poisson, NB, and ZINB (applying mixed effects when we have repeated measures), find the model with the best fit (using AIC/BIC), and use those coefficients and p-values to generate q-values and determine significantly-differently abundant OTUs between the two groups that we're comparing.
My questions are as follows:
(1) If my model looks like
OTU ~ treatment or
OTU ~ treatment + (1|studyid) is AIC or BIC better? The same individuals will be in all comparisons, so I don't think that I'm concerned about differences in sample size. Romero uses AIC, but there also seems to be some argument for BIC, and I'm not sure in which instances one is more correct than the other.
(2) Along the same lines, I assume that I can use AIC/BIC for the mixed-effects versions of the Poisson, NB, and ZINB? I'm using the versions implemented in R using the
(3) Some of my models won't converge, and to this point, I've taken the models that don't converge and consider them to be non-biologically-relevant anyways. I therefore toss out the resulting coefficients/p-values, and only compare the converging models. Does this make sense? Or do I need to consider that the models that didn't converge may still be useful?
I know that these questions can't be validated on a bench, but I'm not sure how to create a script that will be as true as possible to the bacteria that we find in our samples.