Proper approach for modeling multi-class probabilities (proportions/compositional data) I have a dataset in which the target values are counts of objects within each class, example (Xs are input features, Ys are values to be predicted):
example_id| X1, X2, X3, X4 | Y1,   Y2,   Y3,   Y4
----------+----------------+-----------------------
         1|  1,  2,  3,  4 |  1,    0,    5,   8
         2|  5,  6,  7,  8 |  3,   10,    0,   0
         3|  8,  9,  0,  1 |  1,    1,    1,   1

What I'm actually interested in predicting are the proportions of each classes within the sample, that is:
example_id| X1, X2, X3, X4 |   Y1,     Y2,     Y3,     Y4
----------+----------------+------------------------------
         1|  1,  2,  3,  4 | 1/14,   0/14,    5/14,   8/14
         2|  5,  6,  7,  8 | 3/13,  10/13,    0/13,   0/13
         3|  8,  9,  0,  1 |  1/4,    1/4,     1/4,    1/4

My initial idea was to transform the original data into the proportions (so model directly what I'm interested in) and then train the GBM-based model on it with a softmax-based objective that would allow passing in multi-dimensional and continuous input.
What I've found out is that neither of the common GBM libraries actually support that kind of setup, meaning it's either not very common or just inappropriate.
Clearly there are several ways I can try solving that problem, including:

*

*using neural nets instead of GBMs to train using the objective I initially wanted to use

*modeling the counts directly and then getting the proportions in the post-processing step

*use MultiCrossEntropy objective and then transform the scores so they sum up to 1, which seems quite dumb and very inappropriate

I can obviously try all the approaches and see which one works best, but I was wondering whether there's any recommended/theory-driven approach to this kind of problem that I failed to find?

If I understand the data correctly, what you have is a set of 4 mutually exclusive outcome classes. Under reasonable assumptions (in particular, the "independence of irrelevant alternatives") that can be handled by multinomial logistic regression, which can be implemented via a neural net as shown in the link. You then convert to probabilities after modeling.
You have to model counts rather than proportions, as a proportion based on 1000 observations is much more reliable than one based on 10 observations. Although the data format in the linked example above has a single outcome value for each data row, multinomial regression can handle aggregated data like you show. The multinom() function in the nnet package can accept an outcome that is "a matrix with K columns, which will be interpreted as counts for each of K classes" (from the help page). That seems to be just what you have. Then you use your Xs as the predictors for the regression, and ultimately represent the results in the probability scale.
Interfacing with other software like the emmeans package works better if you reformulate the data into a long form. Each data row is transformed into a number of rows equal to the number of outcome categories. The outcome is represented as levels of a categorical variable, and the number of counts with that outcome is provided to the weights argument of the function. There are a couple of worked-through examples on this page.
I'm not sure about how to implement multinomial outcomes with a gradient boosted machine, however.
