# Issues with regression and prediction

I am having a lot of fun with regression analysis at the moment, and by fun I mean bashing myself repeatedly over the head. I have a set of 200 data points, by filtering on a property of interest, I end up with 153 points of use.

I initially used these 153 points to generate a linear regression, with an excellent R${^2}$ and a plot of fitted vs actual variables of almost a perfect diagonal. Great! However, it was suggested that this might only be an internally predictive model (which as I understand it means the model fits the data, rather than the opposite). So, I then tried this: I randomly selected a sample of 100 of the 153 results, and built the same model, it still gave a relatively good fit. I then used the predict function in R to try to predict the outcome of the other 53 records. It did not go well. What I got was one of 2 things.

1. the predictions made no sense at all, not even on the same scale as the actual values.
2. most of the predictions made sense (although weren't very accurate) and one or two, were on an entirely different scale (orders of magnitude larger, or smaller).

Since the model I am fitting has time as the response variable, it was suggested I use a Gamma fit regression instead of a plain old linear regression. I tried this and ended up essentially with the result.

So, am I using R correctly, was Gamma a good choice for this? I'm pretty sure my data is good (non biased) so if I am unable to predict, despite the good model - does this mean my model is useless? I've been working on this for some weeks now, and it would be great if I could salvage something.

The R commands I have used:

modelSet<-sample(1:nrow(myData),100)
modelData<-myData[modelSet,]
predictData<-myData[-modelSet,]

fit<-lm("time~(x1+x2+x3+x4+x5+x6)^3", data=modelData)
pred<-predict(fit, predictData)
plot(predictData$time, pred) <- gives a really not useful plot fit2<-glm("time~(x1+x2+x3+x4+x5+x6)^3", data=modelData, family=Gamma) # tried with link=log too pred2<-predict(fit2, predictData) plot(predictData$time, pred2) <- gives an even less useful plot

• You have a badly overfitted function. The advice you got about having a hold-out subsample and trying out-of-sample prediction has saved you from making a serious error. The Gamma glm advice is probably sensible, but it's no good to you until you're clear on how to go about avoiding overfitting. Consider using cross-validation, for example. – Glen_b Jun 12 '14 at 9:48
• You might benefit from chapter 7 of Elements of Statistical Learning (you'll see the 10th printing of the book is downloadable there) – Glen_b Jun 12 '14 at 9:56

Hidden away in the R code is the information that you're trying to estimate 42 regression coefficients from 153 observations, & doubtless over-fitting. If that model, which includes all two-way & three-way interactions between six predictors, is of special interest you need to collect more observations to fit it well; otherwise fit one more appropriately sized for the number of observations you have—perhaps six for linear terms & just a few likely interactions or non-linear terms.

Rules of thumb say that in most situations where you're wanting to fit a regression model to observational data you'll need at least 10 to 20 observations for each estimated coefficient in addition to the intercept to avoid badly over-fitting it. The hand-outs for @Frank Harrell's Regression Modelling Strategies course explain how to use le Cessie & van Houwelingen's heuristic shrinkage estimator to help decide how many coefficients you can sensibly estimate in a particular case, when the model you came up with at first is over-fitted (§4.7.7, "How Much Data Reduction Is Necessary?").

Your way of checking the predictive ability of the model is based on a sound idea, & seems to have rightly shown up a problem in this case; but the results are going to vary a lot depending on which 53 observations you happen to exclude. Cross-validation splits the sample randomly many times & averages the out-of-sample fit metric, to give a more stable estimate. Note that when you're doing ordinary least-squares regression, the predicted residual sum of squares (PRESS) can be got analytically. In R press <- sum((residuals(fit)/(1 - lm.influence(fit)\$hat))^2).

As @Glen_b says, think about over-fitting first, then read his answer here on gamma GLMs. Fitting a log-normal model would be more straightforward than a gamma GLM with a log link, & I'd guess with so few observations there might be little to choose between them.

• Thanks for the answer, getting more data will be difficult, each trial takes upto 1 hour to run, and the trials being run are from a set of already approved trials. Is cross validation also called bootstrapping? It sounds very much like the next step I was advised to do, perform the model/predict sample 1000 times, and plot the error. – Zack Newsham Jun 12 '14 at 14:56
• Bootstrap validation often gets lumped in with cross-validation as they have the same goals (& give similar results). I try to remember to distinguish them. Fit a model to each bootstrap sample, calculate your fit metric; score out the actual sample with that model, calculate your fit metric: the average difference gives an estimate of optimism, which you subtract from the fit metric of the original sample to give an idea how the model will perform on new data. Made very easy using the rms package. – Scortchi Jun 12 '14 at 15:30
• But don't expect a happy result in this case, just a more accurate measure of how badly the model performs - it's already clear it's far too big. – Scortchi Jun 12 '14 at 16:19