Beta regression - simultaneously tracking all components so that proportions sum to 1.0

I am trying to fit a beta regression model in R with the outcome variable representing proportions, so that the model accounts for all of the constituents of the whole fraction. I am unclear how to keep track of multiple constituents at once. I have looked at a lot of different examples (using, for example, the betareg or glmmTMB packages), but couldn't quite find something that deals with the issue. The examples I have just focus on a single constituent among the proportions. It is also relevant that it should account for the fact that they need to sum to 1.0, although I realise that this over-specifies the problem, because the final proportion is determined by the sum of the others.

Here's a synthetic dataset to try to illustrate what I am trying to do here:

library(dplyr)
library(betareg)
library(glmmTMB)
set.seed(15)

df <- expand.grid(a = 1:5,b = 1:6,c=1:7)
n <- nrow(df)
## constituents of the sum
df$$A <- sin(df$$a) + 1 + runif(n)
df$$B <- cos(df$$b) + 1 + runif(n)
df$$C <- log(df$$c) + runif(n)
## proportions of the sum
df$$pA <- with(df,A/(A+B+C)) df$$pB <- with(df,B/(A+B+C))
df$$pC <- with(df,C/(A+B+C)) ## factor levels to represent which constituent was being sampled LETS <- LETTERS[1:3] df_reshaped <- rbind(cbind(rename(df[,c(1:3,7)], proportion = pA), component = factor(rep('A',n),levels = LETS)), cbind(rename(df[,c(1:3,8)], proportion = pB), component = factor(rep('B',n),levels = LETS)), cbind(rename(df[,c(1:3,9)], proportion = pC), component = factor(rep('C',n),levels = LETS))) ## represent the predictors as categorical variables df_reshaped$$a <- factor(df_reshaped$$a) df_reshaped$$b <- factor(df_reshaped$$b) df_reshaped$$c <- factor(df_reshaped$c) str(df_reshaped) ## output: ## 'data.frame': 630 obs. of 5 variables: ## $$a : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ... ##$$ b : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ... ## $$c : Factor w/ 7 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ... ##$$ proportion: num 0.459 0.489 0.42 0.282 0.119 ... ##$ component : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
## output:
##   a b c proportion component
## 1 1 1 1  0.4585355         A
## 2 2 1 1  0.4888762         A
## 3 3 1 1  0.4196652         A
## 4 4 1 1  0.2824253         A
## 5 5 1 1  0.1194853         A
## 6 1 2 1  0.6089280         A
tail(df_reshaped)
## output:
##     a b c proportion component
## 625 5 5 7  0.5939893         C
## 626 1 6 7  0.3491562         C
## 627 2 6 7  0.2814178         C
## 628 3 6 7  0.3306160         C
## 629 4 6 7  0.3385776         C
## 630 5 6 7  0.4702073         C

## this *does not* work (!!), but aims to capture the meaning of what I am after
mod1 <- betareg(proportion | component ~ a + b + c, data = df_reshaped)
## output:
## Error in betareg(proportion | component ~ a + b + c, data = df_reshaped) :
##   empty model

## this *does* run, but I don't think it represents what I am trying to achieve here
mod2 <- glmmTMB(proportion ~ a + b + c + (1 | component), data = df_reshaped, beta_family())


Neither of the model specifications above really capture what I am aiming to achieve.

Edit: Thanks to an insightful comment from @dimitriy, I was able to find the representation I was after. This uses Dirichlet regression, and was achieved using the DirichletReg package.

library (DirichletReg)

## create 'factor' objects for these unordered categorical variables
df$$fa <- factor(df$$a)
df$$fb <- factor(df$$b)
df$$fc <- factor(df$$c)
df\$Y <- DR_data(df[,c('pA','pB','pC')])

res1 <- DirichReg(Y ~ fa + fb + fc, data = df)
summary(res1)

## Just as a sanity-check, extract the coefficients, and compare them
## for what we might expect based on the construction of the data.
x1 = rowMeans(sapply(coef(res1),'[',2:5));      ## extract components for the 'A' factor
x2 = sin(2:5) - sin(1)                          ## contributions from the 'A' factor
y1 = rowMeans(sapply(coef(res1),'[',1:5 + 5));  ## extract components for the 'B' factor
y2 = cos(2:6) - cos(1)                          ## contributions from the 'B' factor
z1 = rowMeans(sapply(coef(res1),'[',1:6 + 10)); ## extract components for the 'C' factor
z2 = log(2:7)                                   ## contributions from the 'C' factor
cor(x1,x2) ## comparison for the 'A' factor (only 4 coefficients to account for)
## [1] 0.9774279
cor(y1,y2) ## comparison for the 'B' factor (5 coefficients to consider)
## [1] 0.9472137
cor(z1,z2) ## comparison for the 'C' factor (5 coefficients to consider)
## [1] 0.8884948


However, there is a trade-off. There don't appear to be many options for performing Dirichlet regression in R (there's also the zoid package - maybe others too). According to [Douma & Weedon, 2019][1], at the time when that review was conducted there were no R packages that could handle mixed-effects or zero-/one-inflated response variables.

[1]: Douma & Weedon (2019): Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression. Methods in Ecology and Evolution. Vol. 10, Iss. 9, pp. 1412-1430. https://doi.org/10.1111/2041-210X.13234

• How do you wind up with the proportion values? If you can model the discrete events (maybe you can’t, but maybe you can), a multinomial accomplishes your goal of having the sum be one.
– Dave
Aug 25, 2023 at 2:36
• Have you thought about using Dirichlet regression? Aug 25, 2023 at 3:53
• @dimitriy - this was an excellent suggestion, and provides the answer I was after. Thank you! Aug 28, 2023 at 3:17
• @dimitriy - I'm satisfied with your suggestion. I can accept it if you can post it as an Answer (preferably in a more expanded form). Aug 28, 2023 at 6:23
• It's better to add answers as anwers, so that the question can be counted as answered. You can answer your question btw @nullglob Aug 28, 2023 at 7:09