[Rd] F statistic in add1.lm vs add1.glm
Steven McKinney
smckinney at bccrc.ca
Tue Jun 25 21:23:22 CEST 2013
> stats:::add1.glm
function (object, scope, scale = 0, test = c("none", "Rao", "LRT",
"Chisq", "F"), x = NULL, k = 2, ...)
{
Fstat <- function(table, rdf) {
dev <- table$Deviance
df <- table$Df
diff <- pmax(0, (dev[1L] - dev)/df)
Fs <- (diff/df)/(dev/(rdf - df))
Is this where the double division is happening?
diff has df in the denominator, then Fs assignment sees
diff divided by df again. if df is 1, the double division
will go unnoticed.
Steven McKinney
Statistician
Molecular Oncology and Breast Cancer Program
British Columbia Cancer Research Centre
________________________________________
From: r-devel-bounces at r-project.org [r-devel-bounces at r-project.org] On Behalf Of William Dunlap [wdunlap at tibco.com]
Sent: June 25, 2013 12:01 PM
To: r-devel at r-project.org
Subject: [Rd] F statistic in add1.lm vs add1.glm
Should the F statistic be the same when using add1() on models created by lm and glm(family=gaussian)?
They are in the single-degree-of-freedom case but not in the multiple-degree-of-freedom case.
MASS:addterm shows the same discrepancy. It looks like the deviance (==residual sum of squares) gets
divided by the number of degrees of freedom for the term twice in add1.glm. Using anova() on the output
of lm and glm(family=gaussian) gives the save F statistic as add1.lm gives.
> # factor(carb) consumes 5 degrees of freedom, am 1, compare their F values.
> fit <- lm(mpg ~ hp, data=mtcars) ; add1(fit, ~ hp + factor(carb) + am, test="F")
Single term additions
Model:
mpg ~ hp
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 447.67 88.427
factor(carb) 5 103.54 344.13 90.009 1.5044 0.2241
am 1 202.24 245.44 71.194 23.8952 3.46e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> fit <- glm(mpg ~ hp, data=mtcars) ; add1(fit, ~ hp + factor(carb) + am, test="F")
Single term additions
Model:
mpg ~ hp
Df Deviance AIC F value Pr(>F)
<none> 447.67 181.24
factor(carb) 5 344.13 182.82 0.3009 0.9077
am 1 245.44 164.01 23.8952 3.46e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list