Classic linear model - model selection

I have a classic linear model, with 5 possible regressors. They are uncorrelated with one another, and have quite low correlation with the response. I have arrived at a model where 3 of the regressors have significant coefficients for their t statistic (p<0.05). Adding either or both of the remaining 2 variables gives p values >0.05 for the t statistic, for the added variables. This leads me to believe the 3 variable model is "best".

However, using the anova(a,b) command in R where a is the 3 variable model and b is the full model, the p value for the F statistic is < 0.05, which tells me to prefer the full model over the 3 variable model. How can I reconcile these apparent contradictions ?

Thanks PS Edit: Some further background. This is homework so I won't post details, but we are not given details of what the regressors represent - they are just numbered 1 to 5. We are asked to "derive an appropriate model, giving justification".

• An appropriate model might be taken to mean a model that effectively uses all pre-specified variables (accounting for nonlinearity, etc.). I hope your instructor understands that stepwise variable selection is invalid. Few do. – Frank Harrell Oct 17 '11 at 13:56
• Hi again and thanks. Sorry for all the back and forth. The instructions also say "There is not necessarily one "best" model, and you do not necessarily have to include all predictors". Also, there is no collinearity or nonlinearlity. Actually, all the 5 predictors are generated by independent simulations from a normal distribution. Consequently, the correlations between the predictors and the response are also small (the largest is less than 0.1). Frankly my intuition says that the "best" model may just be the sample mean (adjusted r squared is less than 0.03) – P Sellaz Oct 17 '11 at 14:02
• @P Sellaz: given that this is homework using simulated data, your intuition might serve you well here. Write up a well-reasoned explanation for your intuition. – Zach Oct 17 '11 at 16:00
• You can't go by the $R^2$ in general as how large is large is context dependent. But depending on exactly how the simulation was supposed to be executed, you are right that the overall mean may be what's needed. – Frank Harrell Oct 17 '11 at 16:01
• In general it is correct that one does not have to include all predictors to do a good job. But the data are incapable of telling you which predictors to use. – Frank Harrell Oct 17 '11 at 16:01

The problem began when you sought a reduced model and used the data rather than subject matter knowledge to pick the predictors. Stepwise variable selection without simultaneous shinkage to penalize for variable selection, though often used, is an invalid approach. Much has been written about this. There is no reason to trust that the 3-variable model is "best" and there is no reason not to use the original list of pre-specified predictors. P-values computed after using P-values to select variables is not valid. This has been called "double dipping" in the functional imaging literature.

Here is an analogy. Suppose one is interested in comparing 6 treatments, but uses pairwise t-tests to pick which treatments are "different", resulting in a reduced set of 4 treatments. The analyst then tests for an overall difference with 3 degrees of freedom. This F test will have inflated type I error. The original F test with 5 d.f. is quite valid.

• Thanks for your reply. I have added an edit the original question. I hope that is OK. Any further advice would be most welcome. – P Sellaz Oct 17 '11 at 13:37

One answer would be "this cannot be done without subject matter knowledge". Unfortunately, that would likely get you an F on your assignment. Unless I was your professor. Then it would get an A.

But, given your statement that $R^2$ is 0.03 and there are low correlations among all variables, I'm puzzled that any model is significant at all. What is N? I'm guessing it's very large.

Then there's

all the 5 predictors are generated by independent simulations from a normal distribution.

Well, if you KNOW this (that is, your instructor told you) and if by "independent" you mean "not related to the DV" then you know that the best model is one with no predictors, and your intuition is correct.

• Hi Peter, and thanks. N is 900. The data were all produced by simulation. I KNOW this because we had to do the simulatons ourselves. They are supposed to represent real data, as far as this homework is concerned. 100 simulations were conducted, and the 5 with the largest correlations to the response (also simulated but only once) were chosen as the candidate regressors. – P Sellaz Oct 17 '11 at 15:13
• Just be certain that you were to simulate no connection between any X and Y. Then as others have said a regression model is irrelevant and the overall mean is sufficient. – Frank Harrell Oct 17 '11 at 16:00
• Yes, they are completely independent. We chose the data with the largest 5 correlations as the candidate regressors, from which we have to "derive an appropriate model, giving justification" but we "do not necessarily have to include all 5 predictors". – P Sellaz Oct 17 '11 at 16:06
• It sounds like your professor is either a) Completely confused or b) doing something quite interesting. Hard to tell which. If he/she intended this to show the sort of thing that @FrankHarrell and I and others have been pointing out, then good! (that would be b). OTOH, if he/she is intending this to be a "real" regression, then uh-oh it's a). – Peter Flom - Reinstate Monica Oct 17 '11 at 18:12
• I'll let you know which it is when the papers are marked :) – P Sellaz Oct 17 '11 at 18:40

You might try doing cross validation. Choose a subset of your sample, find the "best" model for that subset using F or t tests, then apply it to the full data set (full cross validation can get more complicated than this, but this would be a good start). This helps to alleviate some of the stepwise testing problems.

See A Note on Screening Regression Equations by David Freedman for a cute little simulation of this idea.

I really like the method used in the caret package: recursive feature elimination. You can read more about it in the vignette, but here's the basic process:

The basic idea is to use a criteria (such as t statistics) to eliminate unimportant variables and see how that improves the predictive accuracy of the model. You wrap the entire thing in a resampling loop, such as cross-validation. Here is an example, using a linear model to rank variables in a manner similar to what you've described:

#Setup
set.seed(1)
p1 <- rnorm(50)
p2 <- rnorm(50)
p3 <- rnorm(50)
p4 <- rnorm(50)
p5 <- rnorm(50)
y <- 4*rnorm(50)+p1+p2-p5

#Select Variables
require(caret)
X <- data.frame(p1,p2,p3,p4,p5)
RFE <- rfe(X,y, sizes = seq(1,5), rfeControl = rfeControl(
functions = lmFuncs,
method = "repeatedcv")
)
RFE
plot(RFE)

#Fit linear model and compare
fmla <- as.formula(paste("y ~ ", paste(RFE\$optVariables, collapse= "+")))
fullmodel <- lm(y~p1+p2+p3+p4+p5,data.frame(y,p1,p2,p3,p4,p5))
reducedmodel <- lm(fmla,data.frame(y,p1,p2,p3,p4,p5))
summary(fullmodel)
summary(reducedmodel)


In this example, the algorythm detects that there are 3 "important" variables, but it only gets 2 of them.