model.matrix (via predict) (PR#2206)

ripley@stats.ox.ac.uk ripley@stats.ox.ac.uk
Thu, 24 Oct 2002 10:10:33 +0100 (BST)


I have a clue as to what is going on, given that it changed between 1.5.1
and 1.6.0.  I think it may be related to (NEWS)

        model.matrix() was sometimes incorrectly determining the first
        factor in a formula without an intercept.

It seems a marginal bug, one not worth spending much time on.  Anyone who
has read the internals of model.matrix.default will appreciate my
reluctance.

Brian


On Thu, 24 Oct 2002 Bill.Venables@cmis.csiro.au wrote:

> Glenn,
>
> I suspect that I(i100*j0) will not work, actually.  It certainly won't work
> if either of them are factors, since multiplication, (i.e. the arithmethic
> "*" as opposed to the linear model formula "*" operator) is explicitly
> inhibited for factors.
>
> Even if one is confounded within the other, there is no reason to omit the
> main effect.  Both glm and aov will detect that and silently omit the
> confounded term.  You could also use a/b, of course, and my guess is that
> works as well.
>
> Cheers,
> Bill.
>
> -----Original Message-----
> From: Stone, Glenn (CMIS, North Ryde)
> Sent: Thursday, October 24, 2002 4:53 PM
> To: Venables, Bill (CMIS, Cleveland); 'r-devel@stat.math.ethz.ch '
> Subject: RE: model.matrix (via predict) (PR#2206)
>
>
>
> In my actual application (this is a dummy) i100 is confounded with ihalf, so
> I don't fit its main effect  I guess I(i100*j0) would work?
> -----Original Message-----
> From: Venables, Bill (CMIS, Cleveland)
> To: Stone, Glenn (CMIS, North Ryde); r-devel@stat.math.ethz.ch
> Sent: 24/10/2002 4:29 PM
> Subject: RE: model.matrix (via predict) (PR#2206)
>
> I'm not sure what Glenn's "obvious" workarounds are, but if you fit the
> model in the equivalent and more conventional form:
>
> gg <- glm(clms ~ logj + i100*j0 + ihalf*jmin6, data=dfr, family=poisson)
>
> The predict(gg, newdata = dfr) step now works fine (for me, Windows
> 2000, R 1.6.0).
>
> I present this in the hope someone can use it as a clue to find out
> what's going on.  Sorry, I don't have time to do it myself.
>
> Bill Venables.
>
>
> -----Original Message-----
> From: Glenn.Stone@csiro.au [mailto:Glenn.Stone@csiro.au]
> Sent: Thursday, October 24, 2002 2:33 PM
> To: r-devel@stat.math.ethz.ch
> Cc: R-bugs@biostat.ku.dk
> Subject: model.matrix (via predict) (PR#2206)
>
>
> Full_Name: Glenn Stone
> Version: 1.5.1 and 1.6.0
> OS: win2000
> Submission from: (NULL) (168.140.227.9)
>
>
> The following code produces incorrect fitted values in version 1.5.1 and
> an
> error in 1.6.0
>
>    Error in "contrasts<-"(*tmp*, value = "contr.treatment") :
>            contrasts apply only to factors
>    In addition: Warning message:
>    variable ihalf is not a factor in: model.frame.default(object, data,
> xlev =
> xlev)
>
> The problem in 1.5.1 seems to be model.matrix return inconsistently
> ordered
> interaction terms. There are abviously several work arounds.
>
>
> ### Start of example
> dfr <- expand.grid(i=0:154,j=0:154)
> dfr <- dfr[with(dfr,i+j<=154),]
>
> dfr$logj <- with(dfr,log(j+0.5))
> dfr$jmin6 <- with(dfr, pmin(j,6))
> dfr$ihalf <- with(dfr, factor(floor(i/6)))
> dfr$i100 <- with(dfr, ifelse(i<100,1,0))
> dfr$j0 <- with(dfr, ifelse(j==0,1,0))
>
> Ai <- runif(length(levels(dfr$ihalf)))
> Ai <- Ai-Ai[1]
> Bi <- runif(length(levels(dfr$ihalf)))
> Bi <- (Bi-Bi[1])/10
>
> eta <- with(dfr, -5 -2.5*logj - 3*i100*j0
>             + 0.01*jmin6+Ai[as.numeric(ihalf)]
>             +Bi[as.numeric(ihalf)]*jmin6)
>
> dfr$clms <- sapply( eta , function(x) rpois(1,1e6*exp(x)))
>
> gg <- glm(clms~ logj + i100:j0 + jmin6 + ihalf + ihalf:jmin6,
>           data=dfr,
>           family=poisson)
>
> plot(predict(gg),predict(gg,newdata=dfr))
> ### End of example
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
> -.-.-.-
> r-devel mailing list -- Read
> http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To:
> r-devel-request@stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
> _._._._
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>

-- 
Brian D. Ripley,                  ripley@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 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._