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.