Is it possible to fit a Dirichlet Regression to changing response variables? 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!
 A: 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).
