Is it possible to fit a Dirichlet Regression to changing response variables?

Question

A toy problem to illustrate my issue:

We put a 100 people in a room with 10 candy bars. Each bar is different and has a different brand, flavor, size, color, etc. We ask each person in the group to choose the candy bar that looks the most appealing.

• 100 people
• 10 candy bars
• each person chooses 1 bar

At the end of the experiment we record each person's candy bar choice and figure out what percent of the group choose each candy bar.

obs_1     
Bar 1: 10%
Bar 2: 25%
Bar 3: 14%
Bar 4: 1% 
Bar 5: 2% 
Bar 6: 40%
Bar 7: 2% 
Bar 8: 3% 
Bar 9: 1% 
Bar 10: 2%

Let's say we ran this experiment 10 times with different groups and for the 11th experiment, we wanted to try to forecast what percent of people would choose each candy bar.

Because the sum of all candy choice percentages will be 1 (sum to unity) every time and there is an element of randomness, this seems like a perfect problem to model with a Dirichlet distribution. I'm hopeful that we can forecast choices in future experiments with a Dirichlet Regression (DirichletReg in R) or possibly using STAN.

Question 1: I'm at a loss as to how to format my data for a Dirichlet Regression. In my real data the "candy bars" are constantly being replaced. In one observation I'll have bars 1-10, but in the next observation bars 1 and 4 are replaced by bars 11 and 12 and might look like this:

obs_2      
Bar 11: 9% 
Bar 2: 13% 
Bar 3: 26% 
Bar 12: 22%
Bar 5: 1%  
Bar 6: 20% 
Bar 7: 2%  
Bar 8: 3%  
Bar 9: 1%  
Bar 10: 2% 

My full dataframe looks like this:

+--------------+------------+----------------+-----+------------+----------------+
| obs_number   | bar_1_size | bar_1_chosen_% | ... | bar_15_size| bar_15_chosen_%|
+--------------+------------+----------------+-----+------------+----------------+
| 1            | 15.5       | 10.0           | ... | NULL       | NULL           |
+--------------+------------+----------------+-----+------------+----------------+
| 2            | NULL       | NULL           | ... | NULL       | NULL           |
+--------------+------------+----------------+-----+------------+----------------+
| ...          | ...        | ...            | ... | ...        | ...            |
+--------------+------------+----------------+-----+------------+----------------+
| 11           | 15.5       | NULL           | ... | 33.6       | NULL           |
+--------------+------------+----------------+-----+------------+----------------+

There's a column for every bar's size and chosen_percentage even if it wasn't an option in an experiment(like Bar 15 which wasn't an option in experiment 2, but was an option in experiment 11).

How would I properly format this data for DirichletReg?

Question 2: This is more advanced, but how would I format the regression formula? I'm really lost on this one, I'm not sure if I should use the "common" or "alternative" models or how I would format the formula.

I know this is a long question, but I appreciate you taking the time to read it and any help you all can offer. Cheers!

Answer 1

Disclaimer/ Important Note:

The example below assumes that the ratio of probabilities to select certain bars is completely independent from the specific selection available.

This is likely not the case. For instance the decoy effect is an obvious example that shows that relative probabilities to choose a specific bar may drastically change (even inverse) depending on the options that are available to the participants/consumer.

If you wish to capture these effects then you will need to make the model much more flexible. But then the question is more about how to model preference (the deterministic part of the underlying probabilities) rather than modelling the random behaviour.

The example below is an example how you could proceed once you have a model. It shows how you can use the data to estimate the parameters by finding the parameters that maximize the likelihood function.

But note that to obtain a reasonable idea about the model to use you should first start by exploring/plotting your data, see what it is and how it behaves, before you apply a statistic model (such as assuming a Dirichlet distribution).

---

Logistic regression

You might model your system in terms of the Gibbs Measure (that means a multinomial distributed variable, and not a Dirichlet distributed variable).

Let $P(Y_i = j)$ be the probability that a person $i$ selects candybar $j$ with properties $X_{j}$:

$$P(Y_i = j) = \frac{e^{-\beta X_j}}{\sum_{\forall j} e^{-\beta X_j}}$$

Here the term $e^{-\beta X_j}$ relates to logistic regression (you might wish to model it differently, more complex, as mentioned before) and the term $\sum_{\forall j} e^{-\beta X_j}$ is a normalization constant.

---

Example computation, maximizing log likelihood

This normalization constant is problematic in your situation because the available bars $j$ are not constant every time. I do not know of a package that allows you compute this (it is also not a typical situation). The example below shows how you might do it (manually) by finding the parameters that optimize the likelihood function:

set.seed(1)
# sample size
n <- 1000

# generate theoretic probabilities and properties of bar choice
bar_len <-   c(0.5,1,1.5,2,0.5,1,1.5,2)
bar_sugar <- c(1,1,1,1,2,2,2,2)
p_bars = c(exp(bar_len+bar_sugar)/sum(exp(bar_len+bar_sugar)),0)

# some table for sampling of the bars with different type and number of bars
barsamples <- t(sapply(1:n, FUN = function(x) {r <- 1+rbinom(1,7,0.5);  c(sample(1:8,r),rep(9,8-r))}))
barsamples

# simulating the choice of participants 
barprobs <- matrix(p_bars[barsamples],n)
choice <- sapply(1:n, FUN = function(x) sample(barsamples[x,],1,prob=barprobs[x,]))

#
# estimating model parameters by maximizing likelihood
#

# likelihood function 
# the first argument 'par' will be optimized by the call to the function 'optim'
#
loglike <- function(par, choi = choice, bsamp = barsamples, blen=bar_len, bsug=bar_sugar) {
  # compute theoretic probabilities
  p <- c(exp(par[1]*blen + par[2]*bsug),0)
  # compute probailities
  p_choice <- p[choi]/rowSums(matrix(p[bsamp],n))
  -sum(log(p_choice))
}

optim(c(0,0),loglike)

which in this case has model parameters $1,1$ and will return as result $0.87, 0.91$.

I imagine that in a more typical situation you would not use a linear model as above but instead some sort of neural network that allows to capture more various features than the fixed presupposed (limited) linear model.

---

Dirichlet regression

Dirichlet regression models your variables as a Dirichlet distributed variable. But a multinomial distribution feels more natural to me (probabilities for counts). Your raw data is categorical, and has values 0 or 1. This is not what a Dirichlet distribution describes (it describes a continuous distribution).

Possibly, you wish to model the distribution like a Dirichlet-multinomial distribution.

If you insist on using a Dirichlet distribution then you would need to have outcome variables that describe fractions. E.g. you have several times a group of 100 people and you observe the fractions of bars in those groups. Here you will assume that for a single group the distribution will be multinomial, but the parameters that describe that multinomial distribution are themselves variable (not every group of 100 people is the same) and distributed according to a Dirichlet distribution.

In any case, because you are dealing with a result vector that is not continuous you may not be able to use standard R packages like DirichletReg. Such packages model the outcome variable $Y_i$ for each sample $i$ as a ratio of different fractions ($Y_i$ is a vector) as a function $\beta X_i$, but here $X_i$ are properties of the sample $i$ (describing the group, e.g. age, gender, etc.).

---

Hack example

In those standard packages (which model $Y_i$ as a single function of properties $X_i$ that does not allow different $Y_i$) there is no way to make the vector $Y_i$ with different sizes for different samples/observations $i$.

You might hack it by adding some variables in the matrix $X_i$ that describes whether or not the particular option was present or not, but it is not natural/exact/realistic/precise. (I show this to illustrate the idea, I am not suggesting that you should do this)

Example with multinom package

# make predictor variable 'choicematrix' that is defined by whether or not the 
# particular bar was a choice or not.
sel <- 1-sapply(1:8, FUN <- function(x) rowSums(barsamples == x))
choicematrix <- matrix(rep(0,8*n),n)
for (i in 1:1000) {
  choicematrix[i,choice[i]] <- 1
}

# modeling as a multinomial distribution with 'multinom'
mod <- multinom(choicematrix ~ 1+sel)

# getting intercept term
intercept <- c(0,coef(mod)[,1])
# normalizing
intercept <- exp(intercept)/sum(exp(intercept))

# compare
plot(p_bars[-9],intercept, xlim=c(0,1),ylim=c(0,1))
lines(c(0,1),c(0,1))

This models the probability for each bar to be selected as a function of whether or not certain bars are in the sample. (so the properties $X_i$ refer to the test conditions defining which bars where not in the sample). The intercept will refer to the solution giving a probability to each bar (but not as some function of the properties of that bar).

Answer 2

Bayesian multinomial regression

I would recommend (as stated in the comments) Multinomial Regression, whose conjugate prior is the Dirichlet Distribution, so your posterior will be a Dirichlet Distribution.

In this case, there is no difference between running the experiment 10 times with N different candidates each time, or thinking of this as "1 experiment with 10N" candidates.

If your prior is $Dirichlet(\alpha_{1}, \alpha_{2}, \ldots, \alpha_{10})$, your posterior will be $Dirichlet(\alpha_{1}+n_{1}, \alpha_{2}+n_{2}, \ldots, \alpha_{N}+n_{10})$

In which $n_{1}$ is the number of people who chose candy bar 1 as their favourite, and so forth.

Your prior should probably be (unless you have information that you did not state), that all candy bars are equally popular, and thus $\alpha_{1}=\alpha_{2}=\ldots=\alpha_{10}$, but what value these take, will govern how restrictive you want to be, i.e. how much you wish to penalise as implausible, posteriors in which some candy bars are much more popular than others.

If you set all $\alpha_{i}$ to 1, that's the equivalent of a uniform prior.

In terms of how to prepare your data, for each candy bar, you only need to calculate how many people chose it as their favourite, and then your posterior is analytically known, as above.

Also, if you want to know the posterior for any individual candy bar (i.e. what fraction of the world (call this q) is likely to prefer candy bar k), that is given by the beta distribution, $\beta(q| \alpha_{k} + n_{k}, \sum_{j\neq k} \alpha_{j} + n_{j}$

Bayesian multinomial regression when not sampling all candidates

Edit (as I had originally misread your post/you added more information, see comments): I would still propose a solution which involves Multinomial Regression, leading to a Dirichlet Posterior, which ignores the features of the candy bars (e.g. size, flavour). What underpins all of this, is that if you knew the probability of a member of the general population preferring each of candy bars, assuming they can choose from the full selection, you need to make an assumption of how their choice would vary, if you removed some bars from the pool. If you make an assumption that their preferences are redistributed in some random/proportional way (see below), and not that people who prefer the largest bar are likely to prefer the second largest one if the largest is removed from the pool, then I suggest ignoring the features completely.

For simplicity, let us assume there are only 3 different bars for now.

We know that a random selected member of the population has a probability $q_{i}$ of saying they prefer the $i^{th}$ bar, such that $\sum_{i=1}^{3}q_{i}=1$. We now ask somebody chosen at random to tell us whether they prefer bar 1 or bar 2 (so we don't show them bar 3). In that case, I am going to assume that the probability they will prefer bar 1 is $p_{1}=\frac{q_{1}}{q_{1}+q_{2}}$ and the probability they will prefer bar 2 is given by $p_{2}=\frac{q_{2}}{q_{1}+q_{2}}$

Clearly, this is a simplifying assumption, as another response to this question states, the Decoy Effect makes this potentially a bit shakey, and as I stated above, candy bar 3 might be more similar to 2 than it is to 1. But this assumption might work well in scenarios in which all bars are sufficiently different and it allows for a tractable solution. In words, this solution assumes that if we remove bar 3 from the pool, the people who would have preferred it transfer their preferences to 1 and 2 randomly, but according to the ratio in which people prefer bars 1 and 2.

So given we perform one experiment only, and we want to determine $q_{1},q_{2},q_{3}$ from some data, we want to calculate $P(\underline{q}|D)$, which by Bayes is given by $\frac{P(D|\underline{q})P(\underline{q})}{P(D)}$. The only way this is different to more conventional multinomial regression, is that $P(D|\underline{q})$ is not just a categorical distribution.

In this case, if the first experiment resulted in $n_{11}$ people choosing bar 1 and $n_{12}$ people choosing bar 2 (the first index in the subscript indicating that this was experiment number 1, and the second denoting which bar it was), and $n_{13}=0$ because bar 3 was not an option in experiment 1, then the probability of seeing the data we saw, given some vector q, is given by (up to a constant) $\left(\frac{q_{1}}{q_{1}+q_{2}}\right)^{n_{11}}\left(\frac{q_{2}}{q_{1}+q_{2}}\right)^{n_{12}}$

Using a Dirichlet prior with parameter vector $\underline{\alpha}$ yields:

$P(q|D_{1}) \propto \left(\frac{q_{1}}{q_{1}+q_{2}}\right)^{n_{11}}\left(\frac{q_{2}}{q_{1}+q_{2}}\right)^{n_{12}} q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}$

(I've added a subscript to the D here, to denote this is the data from experiment 1). The proportionality coefficient is hard to calculate, and will have to be calculated numerically, but because there's more to come, let's just treat it as a number tbd later for now.

So now, let's say in experiment 2, we only let them choose between candy bars 2 and 3. Then $P(D_{2}|q)\propto \left(\frac{q_{2}}{q_{2}+q_{3}}\right)^{n_{22}}\left(\frac{q_{3}}{q_{2}+q_{3}}\right)^{n_{23}}$, where $n_{22}$ and $n_{23}$ are the numbers of people who chose bars 2 and 3 in experiment 2 respectively.

So now if you want to know $P(q|D_{2})$ (which technically is $P(q|D_{2}, D_{1})$), you can do the same trick, but now, you use $P(q|D_{1})$ where you would have used your Dirichlet prior $P(q)$ previously.

Explicitly: $P(\underline{q}|D_{1},D_{2})\propto \left(\frac{q_{2}}{q_{2}+q_{3}}\right)^{n_{22}}\left(\frac{q_{3}}{q_{2}+q_{3}}\right)^{n_{23}}\left(\frac{q_{1}}{q_{1}+q_{2}}\right)^{n_{11}}\left(\frac{q_{2}}{q_{1}+q_{2}}\right)^{n_{12}} q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}$

I hope you can see a pattern developing, and how this process can be continued over many experiments with any combination of available candy bars.

There is one final step, which is to turn the constant of proportionality into an equality. I'll outline here for this case, where there were only 2 experiments.

We know that $\int d\underline{q} P(\underline{q}|D_{1},D_{2})=1$, and thus:

$P(\underline{q}|D_{1},D_{2})= \frac{\left(\frac{q_{2}}{q_{2}+q_{3}}\right)^{n_{22}}\left(\frac{q_{3}}{q_{2}+q_{3}}\right)^{n_{23}}\left(\frac{q_{1}}{q_{1}+q_{2}}\right)^{n_{11}}\left(\frac{q_{2}}{q_{1}+q_{2}}\right)^{n_{12}} q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}}{\int d\underline{q} \left(\frac{q_{2}}{q_{2}+q_{3}}\right)^{n_{22}}\left(\frac{q_{3}}{q_{2}+q_{3}}\right)^{n_{23}}\left(\frac{q_{1}}{q_{1}+q_{2}}\right)^{n_{11}}\left(\frac{q_{2}}{q_{1}+q_{2}}\right)^{n_{12}} q_{1}^{\alpha_{1}-1}q_{2}^{\alpha_{2}-1}q_{3}^{\alpha_{3}-1}}$

Note that this integral is over a complex multidimensional surface on which $q_{1}+q_{2}+q_{3}=1$ and $0 \leq q_{i} \leq 1, \forall i$

I suspect this integral cannot be done analytically, so you'll need to solve numerically. I'm no expert on which Monte Carlo sampling method you'll want to use, perhaps somebody else can suggest whether this is best suited to HCMC or other. Whether you need to normalise will however depend on what you want to do with it.

Similarly, you'll need to do more numerical integrals if you want to calculate $\langle q_{i} \rangle$, i.e. the expected fraction of the general population who will prefer bar i.