Bill
I think you are correct, there's something funny in add1, but is it just
degrees of freedom? Example below..
On Tue, May 14, 2013 at 2:23 PM, William Dunlap <wdun...@tibco.com> wrote:
> Shouldn't the F statistic (and p value) for the x2 term in the following
> calls
> to anova() and add1() be the same? I think anova() gets it right and
> add1()
> does not.
>
> > d <- data.frame(y=1:10, x1=log(1:10), x2=replace(1/(1:10), 2:3, NA))
> > anova(lm(y ~ x1 + x2, data=d))
> Analysis of Variance Table
>
> Response: y
> Df Sum Sq Mean Sq F value Pr(>F)
> x1 1 52.905613 52.905613 1108.61455 4.5937e-07 ***
> x2 1 6.355775 6.355775 133.18256 8.5678e-05 ***
> Residuals 5 0.238611 0.047722
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> > add1(lm(y ~ x1, data=d), y ~ x1 + x2, test="F")
> Single term additions
>
> Model:
> y ~ x1
> Df Sum of Sq RSS AIC F value Pr(>F)
> <none> 6.5943869 2.4542182
> x2 1 6.3557755 0.2386114 -22.0988844 186.45559 2.6604e-06 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> Warning message:
> In add1.lm(lm(y ~ x1, data = d), y ~ x1 + x2, test = "F") :
> using the 8/10 rows from a combined fit
>
> It looks like add1 is using 7 instead of 5 for the denominator degrees of
> freedom,
> 7 being the value in the original fit, before the 2 rows containing NA's
> in x2
> were omitted.
>
> > (6.355775/1) / (0.238611/5)
> [1] 133.1827745
> > (6.355775/1) / (0.238611/7)
> [1] 186.4558843
>
>
Here's another view of the same problem.
> m1 <- lm(y ~ x1, data = d)
> m2 <- lm(y ~ x1 + x2, data = d)
> anova(m2, m1, test = "F")
Error in anova.lmlist(object, ...) :
models were not all fitted to the same size of dataset
# It is necessary to refit m1 on the smaller dataset, add1 claims it is
doing it.
> m3 <- lm(y ~ x1, data = model.frame(m2))
> anova(m2, m3, test = "F")
Analysis of Variance Table
Model 1: y ~ x1 + x2
Model 2: y ~ x1
Res.Df RSS Df Sum of Sq F Pr(>F)
1 5 0.239
2 6 6.594 -1 -6.356 133.2 8.57e-05 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
## Now compare add1 against m1 and m3
> add1(m1, y ~ x1 + x2, test = "F")
Single term additions
Model:
y ~ x1
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 6.594 2.454
x2 1 6.356 0.239 -22.099 186.5 2.66e-06 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
Warning message:
In add1.lm(m1, y ~ x1 + x2, test = "F") :
using the 8/10 rows from a combined fit
## That F value is just wrong, isn't it?
> add1(m3, y ~ x1 + x2, test = "F")
Single term additions
Model:
y ~ x1
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 6.594 2.454
x2 1 6.356 0.239 -22.099 133.2 8.57e-05 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
It seems to me that those last 2 uses of add1 ought to return the same
thing, shouldn't they?
pj
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
>
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
--
Paul E. Johnson
Professor, Political Science Assoc. Director
1541 Lilac Lane, Room 504 Center for Research Methods
University of Kansas University of Kansas
http://pj.freefaculty.org http://quant.ku.edu
[[alternative HTML version deleted]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
