How can I predict a continuous variable using only ordinal covariates? I have a quite large data set (approximately 1500 individuals, with very few missing values). My goal would be to predict age (thus, a positive, continuous outcome) using approximately 10 ordinal variables representing biological/developmental indicators. Each ordinal variable has exactly 8 stages. What method would be best suited for such a question?
In particular:

*

*I think there are few methods that natively handle ordinal variables. So, should I treat the ordinal predictors as numerical, or as nominal?

*Is there some way of using simply a linear regression here? (Maybe in combination with some penalization method?)

*Ideally, I would like to get not only a point estimate of age, but also a prediction interval.

I've already tried random forests, which perform quite correctly, but I wonder about possibly better alternatives in this use case.
 A: A problem with treating ordinal variables as numeric/continuous is that it assumes, often incorrectly, that predictor categories are equidistant with respect to their effect on the response variable. So, a change from category 1 to category 2 has the equivalent effect on the response variable as a change from category 2 to 3, etc.
Treating ordinal predictors as categorical/unordered can also be problematic for several reasons. Firstly, the information contained in the natural ordering of the predictor is lost. Secondly, because $n - 1$ coefficients are estimated for each predictor, where $n$ is the number of categories, the estimates can be noisy and the model prone to overfitting, especially when data are sparse within some categories.
Monotonic effects
A compromise position between these two extremes is to model ordinal predictors as having monotonic effects (assuming that predictor effects can be assumed to be truly monotonic and not, for instance, parabolic).
In short, the method involves estimating a coefficient $b$, which expresses, similar to a typical regression coefficient, the expected average difference in the response variable between two adjacent categories of the ordinal predictor. In addition, a simplex vector $\zeta$ is estimated which describes the  normalized distances between consecutive predictor categories, thus providing the shape of the monotonic effect.
A software implementation is available in R package brms. The vignette provides detail and some useful extensions (monotonic interactions, random effects, etc). The original paper by Bürkner & Charpentier is also very good.
