Beta regression for interaction terms without full rank? (R)

Question

I have a series of beta regressions I have performed on the effects of tissue type, year, and age class on the eccentricity of ellipses, which varies between 0 and 1. I have created interaction models and a global model for all combinations of these independent variables except age and year, because not all age classes exist in each year (I think this is what is meant by "not full rank", but here's the error message: "Error in optim(par = start, fn = loglikfun, gr = if (gradient) gradfun else NULL, : non-finite value supplied by optim"). Specifically, I only have the HY class in 2021 and the JU class in 2019. That being said, based on how the distribution of my data breaks down by independent variable, I suspect that an age x year model would be important for capturing the skew in my underlying data (see below images compared to the black line representing uncategorized data density in the last image).

When I plot the fitted distributions for the best and worst models (see image below), none of them appear to capture the left skew very well, and I believe a betareg(eccentricity ~ age*year, data = data) model would, but I can't run that model. So I have two questions: 1) Is there any way to address the lack of full rank and still be able to run the age x year betareg, and if not, 2) should I be choosing a different modeling approach to better capture the wonkiness of the underlying distribution? Generalized Additive Models have been suggested to me...

Answer 1

Answer: Based on your description, I think that your assessment is correct: betareg() cannot estimate all coefficients because there are some factor combinations that do not occur. This is currently not captured by betareg() and thus the user needs to handle it. One way of doing so would be to set up the interaction manually:

data <- transform(data, year_age = factor(interaction(year, age)))

The interaction(year, age) sets up all possible combinations and then the application of factor() reduces it to those combinations that actually occur in data. Afterwards you can do:

betareg(eccentricity ~ year_age, data = data)

Additional remark: The empirical distributions you show do not look like beta distributions but rather like a mixture of different beta distributions with different mean and precision parameters. Ideally, these different means/precisions would be explained by the observed covariates. I would encourage you to explore more model specifications, in particular including covariate effects for the precision parameter as well. If this does not suffice, possibly a finite mixture model of beta regressions might help. See: Grün, Kosmidis, Zeileis (2012). "Extended Beta Regression in R: Shaken, Stirred, Mixed, and Partitioned." Journal of Statistical Software, 48(11), 1-25. doi:10.18637/jss.v048.i11.