3
votes
$\begingroup$

How do I test for Lack Of Fit (F-test) using R? I've seen a similar question, but that was for SPSS and it was just said that is can be easily done in R, but not how.

I know in simple linear regression I would use anova(fm1,fm2), fm1 being my model, fm2 being the same model with x as a factor (if there are several y for x). How do I do it in multiple linear regression?

$\endgroup$
2
  • 2
    $\begingroup$ Can you say which CV question you are referring to? Note that an ANOVA is a multiple regression model, just one w/ only categorical covariates. In R, the code / formulation should be identical for performing a nested model test b/t 2 ANOVA's & 2 MR's. $\endgroup$
    – gung - Reinstate Monica
    Commented Oct 21, 2012 at 15:51
  • 1
    $\begingroup$ stats.stackexchange.com/questions/4762/… is the one I was talking about $\endgroup$
    – lisa
    Commented Oct 21, 2012 at 16:33

2 Answers 2

Reset to default
2
votes
$\begingroup$

As @gung says in the comment, your question title and text conflict. The F-test for joint significance of all parameters in a model is on a single model fit; it is displayed each time you do summary().

Comparisons of models is a whole different ball game -- as the models need to be nested for inference to be valid.

The lmtest adds a number of common econometrics tests for linear models. As an illustration, here is the beginning of examples(lrtest) for using a likelihood-ratio test to compare two nested models:

R>      ## with data from Greene (1993):
R>      data("USDistLag")
R>      usdl <- na.contiguous(cbind(USDistLag, lag(USDistLag, k = -1)))
R>      colnames(usdl) <- c("con", "gnp", "con1", "gnp1")
R>      fm1 <- lm(con ~ gnp + gnp1, data = usdl)
R>      fm2 <- lm(con ~ gnp + con1 + gnp1, data = usdl)
R>      lrtest(fm2, fm1)
Likelihood ratio test

Model 1: con ~ gnp + con1 + gnp1
Model 2: con ~ gnp + gnp1
  #Df LogLik Df Chisq Pr(>Chisq)    
1   5 -56.07                        
2   4 -65.87 -1 19.61   9.52e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
R>
$\endgroup$
7
  • 1
    $\begingroup$ What I want to do is test a model for Lack Of Fit. I do not want to compare it, I just read that it is done this way. That seems to be not true. So the question is: how do I test my multiple regression model for lack of fit? $\endgroup$
    – lisa
    Commented Oct 21, 2012 at 16:35
  • $\begingroup$ You can look at summary(lm(y ~ . , data=X)) which prints the F-test for you. $\endgroup$
    – Dirk Eddelbuettel
    Commented Oct 21, 2012 at 16:36
  • $\begingroup$ but that is f-test for different variables, not for the whole model. am I wrong? $\endgroup$
    – lisa
    Commented Oct 21, 2012 at 16:38
  • $\begingroup$ Well, I fear you are wrong. The F-test corresponds to testing all variables at the same time, and as such can be seen as a test of the whole model. You may want to check the documentation of the lmtest package as well as the references it offers. $\endgroup$
    – Dirk Eddelbuettel
    Commented Oct 21, 2012 at 17:50
  • 5
    $\begingroup$ @DirkEddelbuettel The F-test is not the lack of fit test ! $\endgroup$
    – Stéphane Laurent
    Commented Oct 21, 2012 at 19:37
0
votes
$\begingroup$

You can perform the lack-of-fit test with the alr3 package.

> library(alr3)
> x1 <- c(1,1,1,2,3,3,4,4,4,4)
> x2 <- c(1,2,3,1,2,3,1,2,3,3)
> y <- rnorm(10, x1+x2)
> fit <- lm(y ~ x1+x2)
> pureErrorAnova(fit)
Analysis of Variance Table

Response: y
             Df  Sum Sq Mean Sq F value  Pr(>F)  
x1            1  8.6412  8.6412  53.857 0.08622 .
x2            1 11.9019 11.9019  74.180 0.07359 .
Residuals     7 10.8198  1.5457                  
 Lack of fit  6 10.6593  1.7766  11.073 0.22608  
 Pure Error   1  0.1604  0.1604                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note the equivalence here:

> xx1 <- factor(x1)
> xx2 <- factor(x2)
> fit1 <- lm(y ~ xx1*xx2)
> anova(fit, fit1)
Analysis of Variance Table

Model 1: y ~ x1 + x2
Model 2: y ~ xx1 * xx2
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1      7 10.8198                           
2      1  0.1604  6    10.659 11.073 0.2261
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.