[Rd] Bug in update()? (PR#6902)

[Prof Brian Ripley]

Berwin,

It's fortuitous that this works at all.

The problem is not in update but in glm itself: take a look at fm$formula
after the first fit.  fm$terms has the right formula, but formula() 
extracts the $formula component first.

I am not currently sure what the right fix is, so will not try to fix this 
in R-patched/1.9.1.

Brian

On Fri, 21 May 2004 berwin at maths.uwa.edu.au wrote:

> Dear all,
> 
> I noticed the following while playing around with fitting log-linear
> models to contingency tables using R 1.8.1, but the problem also
> exists under R 1.9.0.
> 
> A reproducible example uses the following contingency table:
> 
> 
> > library(MASS)
> > data(quine)
> > tmp <- with(quine, expand.grid(Eth=levels(Eth), Sex=levels(Sex), 
> +                    Lrn=levels(Lrn), Age=levels(Age)))
> > n <- nrow(quine)
> > quine2 <- with(quine,
> +       data.frame(tmp,
> +       Count=as.vector(tapply(rep(1,n), list(Eth, Sex, Lrn, Age), sum))))
> 
> First fit a saturated model and see which term we can drop:
>     > fm <- glm(Count ~ .^4, quine2, family=poisson)
>     > drop1(fm, test="Chisq")
>     Single term deletions
>     
>     Model:
>     Count ~ .^4
>                     Df  Deviance     AIC     LRT Pr(Chi)
>     <none>             6.661e-16 148.916                
>     Eth:Sex:Lrn:Age  2     0.422 145.338   0.422    0.81
> 
> Drop the four-way interaction and check which term we can drop next:
>     > fm <- update(fm, . ~ . - Eth:Sex:Lrn:Age)
>     > drop1(fm, test="Chisq")
>     Single term deletions
>     
>     Model:
>     Count ~ .
>            Df Deviance     AIC     LRT  Pr(Chi)   
>     <none>      35.893 142.810                    
>     Eth     1   36.332 141.248   0.439 0.507811   
>     Sex     1   37.238 142.154   1.345 0.246237   
>     Lrn     1   37.392 142.308   1.499 0.220842   
>     Age     3   49.818 150.735  13.925 0.003009 **
>     ---
>     Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
> 
> All of a sudden there are only the main effects left and all
> interaction terms have been removed.  I expected to see at this point
> a test for all three-way interactions.  Since I started of with the
> saturated model and used update to remove the four-way interaction,
> all two-way and three-way interactions should still be in the model.
> Or am I missing some subtlety of the "^" and "." operators in model
> formulae? 
> 
> On the other hand, if we specify the saturated model as follows:
> 
>     > fm <- glm(Count ~ Lrn*Eth*Sex*Age, quine2, family=poisson)
> 
> then the above steps have the desired effect:
> 
>     > drop1(fm, test="Chisq")
>     Single term deletions
>     
>     Model:
>     Count ~ Lrn * Eth * Sex * Age
>                     Df  Deviance     AIC     LRT Pr(Chi)
>     <none>             4.219e-15 148.916                
>     Lrn:Eth:Sex:Age  2     0.422 145.338   0.422    0.81
>     > fm <- update(fm, . ~ . - Eth:Sex:Lrn:Age)
>     > drop1(fm, test="Chisq")
>     Single term deletions
>     
>     Model:
>     Count ~ Lrn + Eth + Sex + Age + Lrn:Eth + Lrn:Sex + Eth:Sex + 
>         Lrn:Age + Eth:Age + Sex:Age + Lrn:Eth:Sex + Lrn:Eth:Age + 
>         Lrn:Sex:Age + Eth:Sex:Age
>                 Df Deviance     AIC     LRT  Pr(Chi)   
>     <none>            0.422 145.338                    
>     Lrn:Eth:Sex  1    0.493 143.409   0.071 0.789909   
>     Lrn:Eth:Age  2    0.629 141.545   0.207 0.901471   
>     Lrn:Sex:Age  2   12.728 153.644  12.306 0.002127 **
>     Eth:Sex:Age  3    1.019 139.935   0.597 0.897103   
>     ---
>     Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
> 
> The update command only removed the four-way interaction but left all
> the two-way and three-way interaction in the model.  
> 
> I looked at update.formula to see what the problem might be, but after
> reading the help file of .Internal I stopped looking around.... :)
> 
> Cheers,
> 
>         Berwin
> 
> --please do not edit the information below--
> 
> Version:
>  platform = i686-pc-linux-gnu
>  arch = i686
>  os = linux-gnu
>  system = i686, linux-gnu
>  status = 
>  major = 1
>  minor = 9.0
>  year = 2004
>  month = 04
>  day = 12
>  language = R
> 
> Search Path:
>  .GlobalEnv, package:MASS, package:methods, package:stats, package:graphics, package:utils, Autoloads, package:base
> 
> ______________________________________________
> R-devel at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
> 
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595