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 [EMAIL PROTECTED] 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 > > ______________________________________________ > [EMAIL PROTECTED] mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-devel > > -- Brian D. Ripley, [EMAIL PROTECTED] 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 ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-devel
