1
$\begingroup$

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 ...
head(df_reshaped)
## 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

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Dave
    Aug 25, 2023 at 2:36
  • 3
    $\begingroup$ Have you thought about using Dirichlet regression? $\endgroup$
    – dimitriy
    Aug 25, 2023 at 3:53
  • 1
    $\begingroup$ @dimitriy - this was an excellent suggestion, and provides the answer I was after. Thank you! $\endgroup$
    – nullglob
    Aug 28, 2023 at 3:17
  • 1
    $\begingroup$ @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). $\endgroup$
    – nullglob
    Aug 28, 2023 at 6:23
  • 1
    $\begingroup$ It's better to add answers as anwers, so that the question can be counted as answered. You can answer your question btw @nullglob $\endgroup$
    – Firebug
    Aug 28, 2023 at 7:09

1 Answer 1

2
$\begingroup$

You can use maximum likelihood Dirichlet regression to model outcomes that are proportions in [0,1] that add up to one and covariates.

You can also try a quasi-maximum likelihood fractional multinomial logit model.

There are R and Stata packages for both, and maybe others that I am unaware of. I am unsure if either can accommodate mixed effects or zero-/one-inflated response variables off the shelf. I suggest that you elaborate more on that aspect of the problem in a new question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.