Help me fit this non-linear multiple regression that has defied all previous efforts EDIT: Since making this post, I have followed up with an additional post here.
Summary of the text below: I am working on a model and have tried linear regression, Box Cox transformations and GAM but have not made much progress
Using R, I am currently working on a model to predict the success of minor league baseball players at the major league (MLB) level. The dependent variable, offensive career wins above replacement (oWAR), is a proxy for success at the MLB level and is measured as the sum of offensive contributions for every play the player is involved in over the course of his career (details here - http://www.fangraphs.com/library/misc/war/). The independent variables are z-scored minor league offensive variables for statistics that are thought to be important predictors of success at the major league level including age (players with more success at a younger age tend to be better prospects), strike out rate [SOPct], walk rate [BBrate] and adjusted production (a global measure of offensive production). Additionally, since there are multiple levels of the minor leagues, I have included dummy variables for the minor league level of play (Double A, High A, Low A, Rookie and Short Season with Triple A [the highest level before the major leagues] as the reference variable]). Note: I have re-scaled WAR to to be a variable that goes from 0 to 1. 
The variable scatterplot is as follows:

For reference, the dependent variable, oWAR, has the following plot:

I started with a linear regression oWAR = B1zAge + B2zSOPct + B3zBBPct + B4zAdjProd + B5DoubleA + B6HighA + B7LowA + B8Rookie + B9ShortSeason and obtained the following diagnostics plots: 

There are clear problems with a lack of unbiasedness of the residuals and a lack of random variation. Additionally, the residuals are not normal. The results of the regression are shown below:

Following the advice in a previous thread, I tried a Box-Cox transformation with no success. Next, I tried a GAM with a log link and received these plots:

Original

New Diagnostic Plot

It looks like the splines helped fit the data but the diagnostic plots still show a poor fit. EDIT: I thought I was looking at the residuals vs fitted values originally but I was incorrect. The plot that was originally shown is marked as Original (above) and the plot I uploaded afterwards is marked as New Diagnostic Plot (also above)

The $R^2$ of the model has increased
but the results produced by the command gam.check(myregression, k.rep = 1000) are not that promising. 

Can anyone suggest a next step for this model? I am happy to provide any other information that you think might be useful to understand the progress I've made thus far. Thanks for any help you can provide. 
 A: Very nice work.  I think this situation is a candidate for the proportional odds semiparametric ordinal logistic model.  The lrm function in the R rms package will fit the model.  For now you may want to round $Y$ to have only 100-200 levels.  Soon a new version of rms will be released with a new function orm that efficiently allows for thousands of intercepts in the model, i.e., allows $Y$ to be fully continuous [update: this appeared in 2014].  The proportional odds model $\beta$s are invariant to how $Y$ is transformed.  That means quantiles are invariant also.  When you want a predicted mean, $Y$ is assumed to be on the proper interval scale.
A: I think re-working the dependent variable and model can be fruitful here. Looking at your residuals from the lm(), it seems that the major issue is with the players with a high career WAR (which you defined as sum of all WAR). Notice that your highest predicted (scaled) WAR is 0.15 out of a maximum of 1! I think there's two things with this dependent variable that is exacerbating this issue:  


*

*Players that simply play for longer have more time to collect WAR

*Good players will tend to be kept around longer, and will thus have the opportunity to have that longer time to collect WAR


However in the context of prediction, including time played explicitly as a control (in any manner, whether as a weight, or as the denominator in calculating average career WAR) is counterproductive (also I suspect it's effect would also be non-linear). Thus I suggest modeling time somewhat less explicitly in a mixed model using lme4 or nlme. 
Your dependent variable would be seasonal WAR, and you would have a different number $j=m_i$ of seasons per player $i$. The model would have player as a random effect, and would be along the lines of:
$ sWAR_{ij} = \alpha + \sigma^2_i + \text{<other stuff>} + \varepsilon_{ij}$
With lme4, this would look something like
lmer(sWAR ~ <other stuff> + (1|Player), data=mydata)
You still may need to transform on $sWAR$, but I think this will help out on that feedback loop.
