resettest in lmtest: unexpected behavior on perfect fit The function resettest in the R package lmtest runs the Ramsey RESET test.  A student today brought to my attention that this function behaves oddly when it is handed a regression which displays perfect fit:
library(lmtest)

x <- seq(1,100)
y <- x

reg1 <- lm(y~x)
summary(reg1)

# Should match but don't
summary(lm(y~x+I(x^2)))
resettest(reg1,power=2,type="regressor")

set.seed(12344321)
x <- rnorm(1000)
y <- x

reg2 <- lm(y~x)
summary(reg2)

# Should match but don't
summary(lm(y~x+I(x^2)))
resettest(reg2,power=2,type="regressor")

set.seed(12344321)
x <- seq(1,100)
y <- x + 0.001*rnorm(100)

reg3 <- lm(y~x)
summary(reg3)

# Should match and do, even though almost the same as reg1 example
summary(lm(y~x+I(x^2)))
resettest(reg3,power=2,type="regressor")

In this code sample, the p-value for the x^2 term should match the p-value for the resttest in each of the three pairs in the code.  What actually happens is that they fail to match in examples reg1 and reg2---two cases where the regression fits perfectly.  In example reg3---which is almost identical to example reg1, differing only by the addition of a tiny bit of noise---they do match.
This does not have a lot of practical significance since real world regressions basically never fit perfectly, but it is disconcerting if you are just playing around with the function to convince yourself that you understand what it does.
Does anyone understand why this seemingly undesirable behavior occurs or if I am doing something wrong in my use of resttest?  Maybe the test statistic ends up being calculated as (roundoff error)/(roundoff error)?
 A: This suffers from floating point errors. You can reproduce the result if you scale y like resettest does:
resettest(reg1,power=2,type="regressor")
#
#   RESET test
#
#data:  reg1
#RESET = 5.5247, df1 = 1, df2 = 97, p-value = 0.02078

summary(lm(scale(y)~x+I(x^2)))
#Coefficients:
#              Estimate Std. Error    t value Pr(>|t|)    
#(Intercept) -1.741e+00  7.418e-17 -2.347e+16   <2e-16 ***
#x            3.447e-02  3.390e-18  1.017e+16   <2e-16 ***
#I(x^2)      -7.644e-20  3.252e-20 -2.350e+00   0.0208 *  

The test statistics in resettest is calculated as 
reset <- (df2/df1) * ((res1 - res2)/res2)

where res1 and res2 are the sum of squared residuals of the original model (with scaled y) and the model with additional variable(s), i.e., really small numbers in the case of a perfect fit. They should be zero, but aren't due to floating point arithmetics. If you scale y those small numbers change from the 1e-26 magnitude to the 1e-30 magnitude and you get different results. Since the test statistics isn't defined if res2 is zero, resettest should at least give a warning if res2 gets too small.
In summary, you should tell your student that the software gives spurious p-values in case of a perfect fit due to floating point arithmetics. They should investigate if they have (something close to) a perfect fit before they do the test since the reset test isn't needed in such a case. You can't improve a perfect fit after all.
