From smckinney at bccrc.ca  Tue Jun 25 21:23:22 2013
From: smckinney at bccrc.ca (Steven McKinney)
Date: Tue, 25 Jun 2013 12:23:22 -0700
Subject: [Rd] F statistic in add1.lm vs add1.glm
In-Reply-To: <E66794E69CFDE04D9A70842786030B931C3072F8@PA-MBX01.na.tibco.com>
References: <E66794E69CFDE04D9A70842786030B931C3072F8@PA-MBX01.na.tibco.com>
Message-ID: <DCE81E14EB74504B971DAD4D2DB0356B0CB926DF21@crcmail4.BCCRC.CA>



> 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



[Rd] F statistic in add1.lm vs add1.glm
Steven McKinney