0

I have a community matrix (fish gut microbiome) from a laboratory experiment consisting of 2 treatments, each with 3 independent aquaria, which were sampled 3 times. At each sampling event 6 fish were collected for microbiome sequencing. Since the fish that were sampled from the same tank are not independent from each other, my plan was to use a mixed effects model to test for the effects of treatment and time on microbiome composition with aquarium as random factor, e.g. Y ~ treatment * time + (1|aquarium).

However, while the original sampling design was balanced, it was not possible to obtain data from all sampled fish. Therefore, the number of fish per sampling event and aquarium ranges between 4 - 6. While this does not seem to be a problem in R for a univariate analysis (e.g. using nlme::lme), I could not find a solution for a multivariate response variable.

In the R package vegan, restricted permutation tests for e.g. RDA and PERMANOVA are implemented to account for a dependence of observations by defining which samples will be permuted using the permute::how command. However this only works for balanced data sets, i.e. the same number of observations for each level of the random factor (here: aquaria).

The only option I am aware of to solve this problem, which I would really like avoid, is to delete observation to achieve a balanced design again. Do you know of any other method that also works with unbalanaced data? I welcome any ideas of how to solve this issue.

Thanks!

Christiane Hassenrück
  • 1
  • 1

2 Answers2

1

I only know of one package capable of fitting MV mixed effect regressions; asreml-r, which is a paid package and relatively expensive - I am pretty sure that vegan::adonis2 or the permanova functions are not capable of this; they only do nesting, which is not the same. You could include aquarium as a fixed term to see how much variation actually exists between aquarium though.

The differences in sample size aren't too extraordinary. However, if you are really worried about them you can use multiple imputation (here is a tutorial to do it in R) to generate the missing observations.

What kind of community matrix are you dealing with? Presence-Absence?

Sven Hohenstein
  • 6,285
  • 25
  • 30
  • 39
André.B
  • 1,290
  • 6
  • 20
  • Thanks for the advice! I ran the (permanova) model with aquarium as fixed term (~treatment * time + aquarium). Indeed, this factor is quite important. In the end, I calculated the F statistic for the fixed effect of treatment manually: SumOfSqs(treatment)/df(treatment) divided by SumOfSqs(aquarium)/df(aquarium). If I am not much mistaken, this should provide the correct F statistic as expected in a mixed model. – Christiane Hassenrück Dec 11 '18 at 10:26
  • To answer your other question: I did not try to impute the missing values, as I am not sure how reliable this will be for a matrix with 102 observations (samples) of ~1000 variables (OTUs), where I am missing the values for all variables for 6 observations. I work with sequence counts (relative sequence abundances). – Christiane Hassenrück Dec 11 '18 at 10:30
  • That isn't quite how a mixed model works - they are not calculated in the same way as fixed effects. Here (https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used) is a discussion on this topic that might illuminate what you are after, but they are not easily calculated by hand. Fair enough regarding imputation, but I wouldn't worry about the small difference. Given you are dealing with relative abundances you might actually want a CCA rather than a redundancy analysis. – André.B Dec 12 '18 at 18:13
0

After fighting with the problem some more, in the end I did use an imputation approach: I adjusted the number of samples for all combinations of aquarium (this is the random factor) and sampling timepoint to 6, by randomly picking one (or two) of the available observations depending on how many were missing. Then I ran adonis2 with restricted permutations. I repeated this process 100 times, and then used the average of the R2, F, p output as final result.

Christiane Hassenrück
  • 1
  • 1