R - multinomial logistic regression with relative frequencies as response variable

Question

My colleagues observed in an experiment involving categorical and continuous independent variables, how the species composition changes. Approximately equal numbers of microbes were used in the experiments, and at the end of the trials, the surviving species were genetically determined. Only relative frequencies (not count data) of the individuals of each species could be estimated based on quantitative genetic measures.

However, this also means that my response variable is not binary coded (0/1), and probabilities for survival are already available.

Now my question: Can I still perform a multinomial logistic regression in this case? How should I do it in R? Or: Is there already an implementation for this in R? Thank you for suggestions, ideas, and assistance!

Additional background info

A community of bacteria has been grown together in medium. For every experiment, the community has been initialized with the same composition, i.e. every culture started had the same number of individuals belonging to one of the 12 species contained in the population.

Moreover, we tested whether the combination of the supernatants of all species together either increased or decreased the growth of any of these 12 species when cultured separately. This effect of supernatants has been expressed in either a positive or a negative continuous value.

Furthermore, it could be that the individuals belonging to the same species may build microcolonies/clusters inside the liquid medium. Thus, we either slowly shaked the medium with the bacteria composition or didn't shake it.

In the end we were interested which species surveyed finally in each experimental group. Thus, my colleagues performed qPCRs to quantitatively determine how many characteristic DNA sequences for every bacterium could be found. We used this as a proxy for the number of cells in the medium, and thus calculated relative frequencies as proxies for the species composition.

Since I have been asked to analyze this setup, I'd like to know how I can evaluate whether the final species composition is impacted by which explanatory variable (either numeric or categorical).

Thanks in advance for any suggestions on how to analyze such a kind of multinomial model (species composition = starting composition + supernatant effect + treatment).

Answer 1

The regular multinomial regression model is usable as long as the model-fitting procedure accepts fractional responses and robust standard errors are used. This is known as the fractional multinomial regression model based on a quasi-likelihood function. See their tiny difference in likelihood function construction in "Relationship between the likelihood functions used in fmlogit and mlogit" https://maartenbuis.nl/software/likelihoodFmlogit.pdf.

Since qPCRs counted characteristic DNA sequences for every bacterium among 12 species, you do have count data. When the count is large enough, the process can be represented by a normal distribution instead of Poisson. The challenge here is that the counts of 12 species cultured in the same dish is correlated, positively if they thrive from the same nutrients or excrete activators and negatively if they compete for nutrients or excrete inhibitors. Here seemingly unrelated regression (SUR) is helpful in fitting cross-equation correlation of the error term. This is done in R by systemfit::systemfit() and will generate a 12-by-12 covariance matrix of the error term across species-specific equations. See Henningsen, A., & Hamann, J. D. (2007). systemfit: A package for estimating systems of simultaneous equations in R. Journal of Statistical Software, 23(4). https://doi.org/10.18637/jss.v023.i04.

We could also consider a generalized linear model with random intercepts and maybe also random slopes clustered by dish to represent the within-dish correlation across the 12 species. In R, we can use nlme::lme(count ~ supernatant * shake * species, random = ~ 1 | dish, weights = ~ varIdent(form = ~ 1 | species)), lme4::glmer(count ~ supernatant * shake * species + (1 | dish), family = poisson), and glmmTMB::glmmTMB(count ~ supernatant * shake * species + (1 | dish), family = poisson) and several families for generalized Poisson and negative binomial models.

If the response is better represented by proportions of 12 species in each dish, several options are outlined in slides by Buis (2010) "Analyzing Proportions" https://maartenbuis.nl/presentations/berlin10.pdf.

• We can break the 12 proportions into 11 relative proportions between each species and a reference species. Then we can estimate 11 models of a fractional response ranging between zero and one. For binary proportion models, see my answer https://stats.stackexchange.com/a/637309/284766 that includes R tutorials. The strategies include ordinal regression, beta regression (potentially with zero and one inflation), and binary regression with proportional responses and robust standard errors.
• The implementation of Stata fmlogit https://maartenbuis.nl/software/fmlogit.html for fractional multinomial logit in R is at https://github.com/f1kidd/fmlogit with technical documentation https://github.com/f1kidd/fmlogit/blob/master/Documentation/fmlogit_docs.pdf. It is an extension of the quasi-likelihood approach in binary logit models. Fine if not all species are present in a dish or any species has an observed proportion of zero or one. If you concern that some species are similar and share unobserved traits, you could consider nested logit models and see See https://cran.r-project.org/web/packages/mlogit/vignettes/e2nlogit.html. If you concern that unobserved errors are clustered by time or location, you can consider random-effect models and see mixed-effect models https://cran.r-project.org/web/packages/mlogit/vignettes/e3mxlogit.html. Unfortunately, I do not think the package mlogit accepts fractional responses.
• To model multiple subcompositions together, see the textbook van den Boogaart, K. G. G., & Tolosana-Delgado, R. (2013). Analyzing compositional data with R. Springer. https://link.springer.com/book/10.1007/978-3-642-36809-7. It uses OLS linear regression but transforms the response variable using an isometric log-ratio transformation (ilr) into orthonormal coordinates. Chapter 7 discusses what to do with zeroes and missingness. To understand the transformation procedure, see "Compositional Data Analysis in a Nutshell" by Tolosana-Delgado, http://www.sediment.uni-goettingen.de/staff/tolosana/extra/CoDaNutshell.pdf. Section 6.4.2 shows how to fit multinomial regression using nnet::multinom(), which appears to allow fractional responses. You could also test VGAM::vgam(, family = multinomial()) for nonparametric smoothing curves as predictors and let me know if it allows fractional responses.
• A Dirichlet regression specifies a full continuous multinomial distribution and likelihood function, which is an extension of the beta distribution to multiple proportions. See Maier, M. J. (2014). DirichletReg: Dirichlet regression for compositional data in R (Research Report Series 125; Institute for Statistics and Mathematics). Wirtschaftsuniversität Wien. https://doi.org/10.57938/ad3142d3-2fcd-4c37-aec6-8e0bd7d077e1 with tutorials https://www.psychoco.org/2011/slides/Maier_hdt.pdf and https://cran.r-project.org/web/packages/DirichletReg/vignettes/DirichletReg-vig.html.

Controlling the starting composition as predictors is a good practice in causal inference, instead of taking the difference in composition as the response. How to represent the compositional predictors may need trials. Typical ways of incorporating 11 proportions may have severe multicollinearity that inflates the standard errors a lot. The book "Analyzing compositional data with R" instructs to use the ilr transformation, but this may be difficult to interpret. Using additive log-ratio transformation of compositional predictors seems okay.