model.matrix (via predict) (PR#2206)
Glenn.Stone@csiro.au
Glenn.Stone@csiro.au
Thu, 24 Oct 2002 16:52:34 +1000
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._