I have a response Y that is a percentage ranging between 0-1. My data is nested by taxonomy or evolutionary relationship say phylum/genus/family/species and I have one continuous covariate temp and one categorial covariate fac with levels fac1 & fac2.

I am interested in estimating:

1. is there a difference in Y between fac1 and fac2 (intercept) and how much variance is explained by that
2. does each level of fac responds differently in regard to temp (linearly so slope)
3. is there a difference in Y for each level of my taxonomy and how much variance is explained by those (see varcomp)
4. does each level of my taxonomy responds differently in regard to temp (linearly so slope)

A brute force idea would be to split my data into the lowest taxonomy here species, do a linear beta regression for each species i as betareg(Y(i)~temp) . Then extract slope and intercepts for each speies and group them to a higher taxonomic level per fac and compare the distribution of slopes (intercepts) say, via Kullback-Leibler divergence to a distribution that I get when bootstrapping my Y values. Or compare the distribution of slopes (or interepts) just between taxonomic levels or my factor fac respectively.Or just compare mean slopes and intercepts between taxonomy levels or my factor levels. Not sure is this is a good idea. And also not sure of how to answer the question of how many variance is explained by my taxonomy level, like in nested random mixed effect models.

Another option may be just those mixed models, but how can I include all the aspects I want to test in one model

say I could use the "gamlss" package to do:

library(gamlss)

model<-gamlss(Y~temp*fac re(random=~1|phylum/genus/family/species),family=BE)

But here I see no way to incorporate a random slope or can I do:

model<-gamlss(Y~re(random=~temp*fac|phylum/genus/family/species),family=BE)

but the internal call to lme has some trouble with that and guess this is not the right notation anyways. Is there any way to achive what I want to test, not necessarily with gamlss but any other package that inlcuded nested structures and beta regressions? Thanks!

CodePudding user response：

In glmmTMB, if you have no exact 0 or 1 values in your response, something like this should work:

library(glmmTMB)
glmmTMB(Y ~ temp*fac   (1   temp | phylum/genus/family/species),
           data = ...,
           family = beta_family)

• if you have zero values, you will need to do something . For example, you can add a zero-inflation term in glmmTMB; brms can handle zero-one-inflated Beta responses; you can "squeeze" the 0/1 values in a little bit (see the appendix of Smithson and Verkuilen's paper on Beta regression). If you have only a few 0/1 values it won't matter very much what you do. If you have a lot, you'll need to spend some serious time thinking about what they mean, which will influence how you handle them. Do they represent censoring (i.e. values that aren't exactly 0/1 but are too close to the borders to measure the difference)? Are they a qualitatively different response? etc. ...)
• As I said in my comment, computing variance components for GLMMs is pretty tricky - there's not necessarily an easy decomposition, e.g. see here. However, you can compute the variances of intercept and slope at each taxonomic level and compare them (and you can use the standard deviations to compare with the magnitudes of the fixed effects ...)
• The model given here might be pretty demanding, depending on the size of your phylogeny - for example, you might not have enough replication at the phylum level (in which case you could fit the model ~ temp*(fac phylum) (1 temp | phylum:(genus/family/species)), i.e. pull out the phylum effects as fixed effects).
• This is assuming that you're willing to assume that the effects of fac, and its interaction with temp, do not vary across the phylogeny ...