I would like to fit a multivariate regression model of an ordinal random vector using ordinal variables as covariates. I am wondering if there is anything implemented in R or other software that could do that. Or, alternatively using a Bayesian approach.
I would use a Multinomial Generalized Linear Model for this, or rather a collection of them, one for each response. There do also exist Multivariate GLMs, but I don't know of an R package to do that and I also have never seen someone do Multivariate Multinomial GLM, but it should be possible. That said, the benefit of Multivariate GLM is usually not huge over just fitting a separate GLM for each response. Starting from there, then the important point is that you need a Multinomial GLM.
Multinomial GLM is suitible for an ordinal response variable. You set this up effectively by doing an ordered sequence of binary responses, so you have
$g(\mu_m) = X\beta + c_0 +c_1 +\cdots + c_m$
and where $\mu_m = P(Y<m)$
Now you have a choice of link function - probit or logit are the main candidates, and lead to the ordered probit and order logit models. Logit is "perfect" if you fulfil the proportional odds assumption, which basically says that your ordinal levels are evenly spaced. That is, the model will be the same between any pairs of responses equally spaced. So if your responses are A to E, then the proportional odds between A and B will be the same as between D and E. Probit, as usual, is particularly suited when you have "hidden gaussianity", in this case you are going to be thinking about a continuous variable that you are observing as ordinal (exam grades might be a good example).
glm function in R doesn't do this for you, but there is a similar command in the
MASS package called
polr. It works basically the same as
glm: see here for more info: http://stat.ethz.ch/R-manual/R-patched/library/MASS/html/polr.html
As to the issue of ordinal inputs, the usual method is to map them to "sensible numbers" unless the number of levels is low, in which case factors are probably easier. There is, however, an excellent answer on this subject here: https://stats.stackexchange.com/a/5398/19879
Also you can fit an Ordinal Logit model or a Beyesian Ordinal Probit model (amongst others) using the package Zelig. It's a very well documented package. Type
demo(Zelig::oprobit.bayes) to see a demo of each model.