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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._