This question is the first part of a larger question that is continued here. I thought that it could be easier to split it into two questions to generate better answer to both and to help further people encounter them.
So, my problem. I have a series of observations of two variables, let's say $x$ and $y$ and I expect that $y$ is some polynomial function of $x$ (for ex $y = \alpha x + \beta x^3$). However, I'm not sure of how many polynomial degrees should be used, and if some of the smaller can be removed. From the literature I expect that $x^3$ should be the maximum needed.
My overall strategy was to fit a model with all polynomial up to $x^4$ (as in $y = a + bx + cx^2 + dx^3 + fx^4$) and compute it's AIC. Then I tried to remove each individual component in the previous equation and compute the AIC of each fit. I select the one with lower AIC and if it's lower than the AIC with all component I pick this as my new model and restart the process. My initial question here is, thus this strategy seems reasonable?
Now, going back to the problem in the title, I've been an huge fan of robust methods for fitting lately. Specifically, with the help of R's robustbase package I've changed all my fitting for robust methods and found and overall increase in performance. However, I'm not sure how to do it here, or if it's even correct. How one computes an AIC within the robust framework? What happens when one model find's $0$ outliers and another $10$? Intuitively, the fitting of the model with $10$ outliers can have a better score (lower sum square of residuals and higher likelihood) but that doesn't mean it's a better model!