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.

• I found the code in this excellent primer on GAM in R - www3.nd.edu/~mclark19/learn/GAMS.pdf The code: library(car) scatterplotMatrix(mydata[,c(1,1:8)],pch=19,cex=.5,reg.line=F, lwd.smooth=1.25, spread=F,ellipse=T, col=c('gray60','#2957FF','#FF8000'), col.axis='gray50') Jun 12, 2013 at 14:05
• Can you share your dataset? Also, +1 for that scatterplot matrix. It's excellent.
– Zach
Jun 12, 2013 at 14:10
• That's too bad, it looks like an interesting dataset. My suggestion would be to try some other machine learning algorithms, e.g. a random forest.
– Zach
Jun 12, 2013 at 14:24
• Random forests are based on decision trees. Check out the randomForest package in R, and the random forest wikipedia page: en.wikipedia.org/wiki/Random_forest
– Zach
Jun 12, 2013 at 14:47
• "The dependent variable...is measured as the sum of offensive contributions for every play the player is involved in over the course of his career." This jumps out to me. A serious confounder here is how long a player has been playing, both in that [a] longer playtime means more time to "collect" oWAR [b] better players will probably play for longer periods of time. Jun 13, 2013 at 14:00

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.

• By levels, do you mean binning the Y variable into 100-200 buckets? If so, is there any preferred method for choosing the bin size? Should they be equally sized? Jun 12, 2013 at 14:20
• Just do the binning temporarily unless we have the continuous solution. You can bin into 100 percentiles, e.g. require(Hmisc); cut2(y, g=100, levels.mean=TRUE) Jun 12, 2013 at 22:45
• When you say a new version of rms will be released soon, do you have any idea when that might be? Jun 13, 2013 at 14:30
• If you use linux I can give it to you now, otherwise expect 2 weeks. Jun 13, 2013 at 15:10
• I do not use Linux so I guess I'll have to wait. Please let me know when it is available. Jun 14, 2013 at 12:59

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.

• I'm not sure I fully understand. If the dependent variable is seasonal WAR, what are the independent variables? An identical minor league stat line for each player? Are we essentially saying that minor league stat line A can lead to major league WAR B, C, D and E? Jun 14, 2013 at 12:46
• Also, since posting this model, I have followed up with an additional post you may want to check out here: stats.stackexchange.com/questions/61711/… Jun 14, 2013 at 12:47