I am using the step() function to select a model using backward
elimination, with AIC as the selection criterion.  The full regression
model contains three predictors, plus all the second order terms and
two-way interactions.  The full model is fit via lm() using two different
model formulae.  One formula uses explicitly defined variables for the
second-order and interaction terms and the other formula uses the I(x^2)
and colon operators.  The fit generated by lm() is exactly the same for
both models, but when I pass these fitted models to the step() function, I
get two different results.  Apparently, step() does not recognize the three
main predictors unless the second order and interaction terms are
explicitly defined as separate variables.

I assigned this problem to my first-year graduate students, not realizing
that R would give two different answers.  Now I have to re-grade their
homework, but I would really like to give them a reasonable explanation for
the discrepancy.

The complete code is given below.

Could anyone shed some light on this mystery?

Thanks in advance,
Karen Keating
Kansas State University


# Exercise 9.13, Kutner, Nachtsheim, Neter & Li
temp<- scan()
49.0   45.0   36.0   45.0
55.0   30.0   28.0   40.0
85.0   11.0   16.0   42.0
32.0   30.0   46.0   40.0
26.0   39.0   76.0   43.0
28.0   42.0   78.0   27.0
95.0   17.0   24.0   36.0
26.0   63.0   80.0   42.0
74.0   25.0   12.0   52.0
37.0   32.0   27.0   35.0
31.0   37.0   37.0   55.0
49.0   29.0   34.0   47.0
38.0   26.0   32.0   28.0
41.0   38.0   45.0   30.0
12.0   38.0   99.0   26.0
44.0   25.0   38.0   47.0
29.0   27.0   51.0   44.0
40.0   37.0   32.0   54.0
31.0   34.0   40.0   36.0

dat<- matrix(temp,ncol=4,nrow=length(temp)/4,byrow=T)
colnames(dat)<-c('Y','X1','X2','X3')
dat <- data.frame(dat)
attach(dat)

# second order terms and interactions
X12<-X1*X2
X13<-X1*X3
X23<-X2*X3
X1sq <- X1^2
X2sq <- X2^2
X3sq <- X3^2

fit1 <- lm(Y~ X1sq  + X2sq  + X3sq  +X1+X2+X3+ X12 + X13 + X23 )
fit2 <- lm(Y~I(X1^2)+I(X2^2)+I(X3^2)+X1+X2+X3+X1:X2+X1:X3+X2:X3)
sum( abs(fit1$res - fit2$res) ) # 0, so fitted models are the same
dim(model.matrix(fit1))   # 19 x 10
dim(model.matrix(fit2))   # 19 x 10

dim(fit1$model)  # 19 x 10
dim(fit2$model)  # 19 x 7 -- could this cause the discrepancy?

back1 <- step(fit1,direction='backward')
back2 <- step(fit2,direction='backward')
# Note that 'back1' considers the three primary predictors X1, X2 and X3,
# while 'back2' does not.

	[[alternative HTML version deleted]]

