Issues with regression and prediction

Question

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

Answer 1

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.