Fitting a continuous result to categorical predictors semiparametrically Suppose one has a relatively large number of observations, each of which consists of a continuous result and a small number (2 or perhaps 3) of categorical variables, each of which has a large cardinality (up to 10k or so).  So, for instance, one might wish to predict income given a person's university and chosen industry (this is just an example).  Right, now I am essentially fitting a model that looks like
log(income) ~ 0 + university + industry 
in R notation, so each university and each industry has an associated "score", but there really is no reason to think the effect is additive, or that taking the log of income is the best transformation.  Still, it seems to work (somewhat).  
Another idea I would like to try out is to express log(income) as a sum of all possible polynomial powers (up to total degree 4 or so) of the associated scores - this would require a lot of regularization, I suppose.
However, I suppose what I would really like is just to fit something like
g(income) ~ f(score.industry,score.university) and simultaneously fit g,f,score.industry,score.university where the latter two variables are actually vectors.  Does anyone have any insight as to how this might be done?
 A: With up to 10,000 levels of each predictor, you are liable to have huge standard errors around each coefficient--plus, who is going to read the output containing all these results, or fashion a predictive equation out of them?  Instead of using a general linear model like you have, consider looking into data mining techniques like neural networks.  I think it would be more practical.
A: The proportional odds model is semiparametric and robust, and you can use quadratic (and other) penalization to deal with the large number of categories on the right hand side.  This is handled in the R rms package lrm function.  For reasons of computation time you might group income into 200 quantile groups so that you have only 199 intercepts in the model.  Someday I'll make the code more efficient so that no grouping is required at all.  JMP software makes use of patterned matrices so that the prop. odds model runs quickly no matter how many unique values there are in Y.
Quantile regression would be another possibility, although I don't know how to handle the large number of predictor categories.
