# How to calculate model error (like MSE) for a multivariate proportional response?

I have data where the response is multivariate and proportional (rows [observations] sum to 1). I am modelling this response using a Dirichlet regression via the DirichletReg R package where the predictors are the first $m$ principal components of a larger set of variables.

I could use AIC/BIC to select the model size (the value for $m$), but as this is a "prediction" exercise, I'd like to be able to compute some form of error for each of the $m$ models. If this were a simple linear regression I'd use the mean squared error or its square root.

In my data I have 4 proportions measured for each response and I can derive the predicted proportions for each training sample given the fitted model. What would a suitable measure of "error" be for a multivariate proportional response?

### Further details

The response variable is the proportion of each of 4 vegetation types around a lake, for a set of $n$ lakes. The predictors are the first $1, ..., m$ principal components of a data set of counts of pollen grains on $p$ species. The aim is to predict the proportions of the 4 vegetation types for samples where only the pollen counts are know (like the past 10,000 years). I'd like to compute an "error" statistic based on the lack of fit to the response given the model. Then I would use a $k$-fold CV to select the "best" model for predicting the proportions of the vegetation types.

• Could you use the (out-of-sample) predicted likelihood or some function thereof (e.g., -2 loglik or proper scoring rules)? The predictive likelihood is determined in full by the parameters you estimate and you could simply "plug in" the real values. A scoring rules provides a summary measures for the evaluation of probabilistic forecasts, by assigning a numerical score based on the predictive distribution P and on the event or value that materializes x. In your case you have the Dirichlet for P and might use the logarithmic score S(P,x)=log(px) (which is related to entropy btw.) – Momo Feb 19 '14 at 0:27
• @Momo Thanks for this suggestion. I'm not familiar with scoring rules so will have to do a bit of reading up on those. A colleague on Twitter suggested the out-of-sample predicted likelihood (or deviance), which I was able to use quite effectively in the CV to select the model size. This coincided with the same general level of complexity as indicated with BIC (which I presume is not surprising). – Gavin Simpson Feb 19 '14 at 14:54
• Perhaps what I didn't convey too well in the question is that whilst the out-of-sample likelihood is fine for selection, it is less useful for the end user in terms of the range of values around the estimated proportions for the predicted values. I thought of two ways to get at that; i) non-parametric bootstrap the selected model and for each model predict the new responses to derive a sampling distribution for each prediction, then form percentile CI, or ii) use the asymptotic distribution of the model parameters & simulate from their distribution to get a new model & predict using each one. – Gavin Simpson Feb 19 '14 at 14:57
• Sounds good to me. Re Scoring rules: They are in essence just convex functions of the predictive likelihood or its moments that have the property that they are minimal when the real value has been predicted. The squared error in a linear model is an example. They might not be what you need judging from the last comment as they are often not on a scale that is intuitive (there advantage lies in model selection not model communication imo). See. e.g., my answer here stats.stackexchange.com/questions/71720/… – Momo Feb 19 '14 at 15:02

Thumbs up for the fact that you are prepared to use cross validation. All too often researchers (including myself) just rely on model choice statistics such as AIC and BIC.

I think there is no simple answer as how to define a prediction error. Ultimately, it depends on your loss function. What is the cost associated with wrong prediction? Is that cost same for e.g. zero and 50 % frequencies?

Taking this principle seriously might lead into quite special error metrics. For example, you might say that the prediction is ok, if it reproduces the training proportions by ten percent accuracy, and otherwise it's worthless. Thus, you would obtain a binary metric. That might be quite close to how you really judge predictions in biological papers?

Anyway, rules such as square error, entropy or logarithmic error are ultimately arbitrary. They may have convenient properties in some parametric models, but they are not God's word.

• Thanks for this. I think if we take this idea much further, this (loss) is a key question to pose to the people interested in the proportions of past vegetation. At the moment I seem to be interested in 2 things; i) having some "error" to use in CV (for which I have now tried the out-of-sample log-lik, see comments on Q), and ii) a way to show range/uncertainty in predicted proportions to help the pollen experts interpret the predictions. For the latter I think a non-parametric bootstrap might be the easiest way to get sampling distribution for each predicted proportion & show percentile CI. – Gavin Simpson Feb 19 '14 at 15:09
• Non-parametric bootstrap is the easiest way, and probably sufficient. If you want to be a hipster, however, you can use the model estimates for prediction. I don't know what DirichletReg does, but I suspect that the estimates of parameters are based on a normal approximation - so sample from those normal distributions 1,000 time, and calculate the predicted proportions for each realization. – user40541 Feb 20 '14 at 14:40
• Which your strategy ii) above seems to be. It's a good question, whcich one is more accurate, non-parametric bootstrap or parametric simulation based on asymptotic distribution of estimators... One option would be to sample from the exact distribution of the predictions. In a Bayesian framework, these would be Dirichlet posteriors, but I don't know if the presence of covariates confounds this. Intuitively, I would say that this is not the thing you should worry about. Sure there are bigger issues with your model or data. – user40541 Feb 20 '14 at 14:47