# AICc is picking overly complex models - something stricter?

I'd like to know if there are stricter alternatives to automated model selection than AICc / AIC / BIC.

We have approximately ten thousand curves, and for each we'd like to find the most parsimonious Gaussian mixture model. Our maximum complexity is a mixture of 5 Gaussians. The procedure so far is to fit these 5 models:

gauss1 = a1*exp(-((x-b1)/c1)^2)
gauss2 = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)
...
gauss5 = a1*exp(-((x-b1)/c1)^2) + ... + a5*exp(-((x-b5)/c5)^2),


and then calculate AICc for each and choose the model with the lowest AICc value. We repeat this for each of the ten thousand curves.

To verify whether it's working, I generated some simulation data which looks a lot like our real data. The simulated data is generated by models gauss1, gauss2, ..., gauss5 + noise. To include noise, I added a random number to each point (between 0 and 10% of the maximum value of the simulated data, uniform distribution). This way I know what the model complexity should be.

Unfortunately it looks like AICc is choosing overly complex models. An example fit for simulated data is below. The complexity of the generating model is 2 Gaussians (gauss2) and the best fitted model, as chosen by AICc, was gauss5.

Using 50 simulated curves, here is the model complexity breakdown. The simulated curves were chosen to be biased toward a simpler model. However the complexity of the fitted models, as chosen by AICc, is biased toward the more complex. For example from these 50 curves, AICc never chose gauss1, even though 19 out of 50 were generated by gauss1. Similar results are seen for AICc, AIC, and BIC.

Choosing between the 5 models with AICc (or AIC or BIC) does not look like it's working for us. Are there other accepted ways to pick between a small number of models? I don't think I can use LASSO / ridge regression since I believe that's only for linear regression. I also don't want to build something ad hoc.

• Would be good to know how you fit the models and how did you computed AICc, – AlefSin May 24 '16 at 20:27
• in particular: are you assuming Gaussian noise when your noise is actually uniform? (out of curiosity: is BIC more conservative than AICc, even if not as conservative as you'd like?) – Ben Bolker May 24 '16 at 20:35
• @AlefSin Models were fit using 'fit' in Matlab along with 'fitOptions' to specify the model. AIC is calculated as N * log(MSE) + 2*k, BIC = N * log(MSE) + k*log(N), and AICc = AIC + 2k(k+1)/(N-k-1). – R Greg Stacey May 24 '16 at 20:41
• The worst problem with information measures is likely just what you are experiencing. Why not enforce your own regularization by limiting the number of gaussians. – mandata May 24 '16 at 20:58
• @BenBolker Yes, you're right. For the simulation I'm adding uniformly distributed noise, but the fitting assumes gaussian noise. I'll try adding normally distributed noise... Re: BIC, it doesn't seem to be any more conservative than AIC or AICc. (I just checked for the same simulation data, and the distribution in the second figure looks identical.) – R Greg Stacey May 24 '16 at 21:01

You could try , fitting all five models for each fold and picking the model that has the lowest MSE (or whatever other error measure you are interested in) on the holdout folds. That is, for each curve, label all dots at random with labels "1", "2", ..., "10" (for ten-fold cross validation). Fit a one-Gaussian model to all points not labeled "1", and using this fitted model, predict the $y$ component of the holdout data (i.e., the points labeled "1"). Record the Mean Squared Error. Do the same for two-Gaussian, three-Gaussian etc. models. Then repeat the process for points labeled "2". Finally, summarize the MSEs per mixture model over all labels, and you will get grand MSEs per mixture model. Pick the model with the lowest overall MSE.

That said, maybe it's not all that problematic if you have "too complicated" models? The problem with too complex models is usually that they are too variable, but looking at your example with two Gaussians in the true DGP and five in the AICc-fitted model, it's not like the five-Gaussian fit is obviously problematic. This will depend on what you actually want to do with your data. In the particular example you show, it looks to me like (say) an interpolation based on five Gaussians wouldn't be overly dramatically different from one based on two Gaussians.

Edit: Alternatively, you could do a Bayesian approach. Assign prior probabilities to the number of mixture components, which are higher for one or two than for five Gaussians. Next, assign sufficiently uninformative priors to the actual parameters of these $k$ Gaussians. Derive posterior distributions and pick the number of Gaussians with the highest posterior probability, which would be a kind of a maximum a posteriori (MAP) approach. (Or use the full posterior distribution, which would be a kind of a "mixture of mixtures", since it would mix one-Gaussian, two-Gaussian and so forth mixtures.) This approach has the advantage that you can fine-tune the (non-)complexity you want, by changing the prior probabilities you assign to the numbers of mixture components.

• Thanks! The problem might not be AICc per se, and instead that we'd like something stricter. At least for our purposes, 2 Gaussians in the top figure is "more correct" than 5 Gaussians. Ideally there would be a tuneable solution with a parameter that we could adjust to make the selection as strict as we like. I can imagine how to do this, but it would be ad hoc... – R Greg Stacey May 24 '16 at 20:48
• If we used cross-validation, what would it look like here? Would the training and testing data be subsets of a single curve? – R Greg Stacey May 24 '16 at 20:53
• I edited my post to include a Bayesian approach, and to clarify the cross-validation idea. Yes, you'd do this independently for each curve (unless you had external information that told you that all curves had the same underlying mixture DGP, but that doesn't sound like what you are describing). The result would be one mixture model for each curve. – Stephan Kolassa May 24 '16 at 21:01
• Hm. Thank you for accepting my answer, but often questions with accepted answers don't get much attention any more. Please feel free to un-accept and see whether something better comes along - if you still want to accept in a day or two, be my guest ;-) – Stephan Kolassa May 24 '16 at 21:16
• Ah, thank you! I did just that, but only because you suggested... =) Thanks again. – R Greg Stacey May 24 '16 at 21:19

K fold cross validation using the one standard deviation rule. With the sd rule, the model with the lowest error is not selected. Rather, the least complex model that is within one sd of the model with the lowest error is select as the best model. Please see the r package bestglm.

• Thanks, that looks very helpful. Have you used e.g. BICq, BICg? What method in bestglm does your answer correspond to (just to save myself some reading)? – R Greg Stacey May 26 '16 at 20:57
• I used those quantities in bestglm, all of which picked an overly complex model. After reading about k-fold cross validation and the bestglm vignette, I decided it was best for what I was trying to accomplish. – D_Williams May 26 '16 at 21:30