Hello, I have a question about the formula and the user defined function: I can do following: ###Case 1:> clotting <- data.frame(+ u = c(5,10,15,20,30,40,60,80,100), + lot1 = c(118,58,42,35,27,25,21,19,18), + lot2 = c(69,35,26,21,18,16,13,12,12))> g1=glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma) > dc=clotting > dc$u=1 > predict(g1,dc)1 2 3 4 5 6 7 8 9 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 However, if I just simply wrap the poly as a user defined function ( in reality I would have my own more complex function) then I will get error: ###Case 2:> xpoly<-function(x,degree=1){poly(x,degree)} > g2=glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma) > predict(g2,dc)Error in poly(x, degree) : 'degree' must be less than number of unique points It seems that the predict always treat the user defined function in the formula with I(). My question is how can I get the results for Case2 same as case1? Anyone can have any idea about this? Thank you very much. Alex [[alternative HTML version deleted]]
Read about the 'makepredictcall' generic function. There is a method,
makepredictcall.poly(), for poly() that attaches the polynomial coefficients
used during the fitting procedure to the call to poly() that predict()
makes.
You ought to supply a similar method for your xpoly(), and xpoly() needs to
return an object of a a new class that will cause that method to be used.
E.g.,
xpoly <- function(x,degree=1,...){ ret <- poly(x,degree=degree,...);
class(ret) <- "xpoly" ; ret }
makepredictcall.xpoly <- function (var, call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(xpoly))) {
if (!is.null(tmp <- attr(var, "coefs"))) {
call$coefs <- tmp
}
}
}
call
}
g2 <- glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
predict(g2,dc)
# 1 2 3 4 5
#-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
#-0.01398928608
# 6 7 8 9
#-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
You can see the effects of makepredictcall() in the 'terms' component
of glm's output. The 'variables' attribute of it gives the original
function
calls and the 'predvars' attribute gives the calls to be used for
prediction:
> attr(g2$terms, "variables")
list(lot1, log(u), xpoly(u, 1))
> attr(g2$terms, "predvars")
list(lot1, log(u), xpoly(u, 1, coefs = list(alpha = 40, norm2 = c(1,
9, 8850))))
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Thu, Jul 16, 2015 at 12:35 PM, Kunshan Yin <yinkunshan at gmail.com>
wrote:
> Hello, I have a question about the formula and the user defined function:
>
> I can do following:
> ###Case 1:
> > clotting <- data.frame(
> + u = c(5,10,15,20,30,40,60,80,100),
> + lot1 = c(118,58,42,35,27,25,21,19,18),
> + lot2 = c(69,35,26,21,18,16,13,12,12))
> > g1=glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma)
> > dc=clotting
> > dc$u=1
> > predict(g1,dc)
> 1 2 3 4 5
> 6 7 8 9
> -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929
> -0.01398929 -0.01398929 -0.01398929
>
> However, if I just simply wrap the poly as a user defined function ( in
> reality I would have my own more complex function) then I will get error:
> ###Case 2:
> > xpoly<-function(x,degree=1){poly(x,degree)}
> > g2=glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
> > predict(g2,dc)
> Error in poly(x, degree) :
> 'degree' must be less than number of unique points
>
> It seems that the predict always treat the user defined function in the
> formula with I(). My question is how can I get the results for Case2 same
> as case1?
>
> Anyone can have any idea about this?
>
> Thank you very much.
>
> Alex
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
Thanks Bill for your quick reply.
I tried your solution and it did work for the simple user defined function
xploly. But when I try with other function, it gave me error again:
OPoly<-function(x,degree=1,weight=1){
weight=round(weight,0)# weight need to be integer
if(length(weight)!=length(x))weight=rep(1,length(x))
p=poly(4*(rep(x,weight)-mean(range(x)))/diff(range(x)),degree)
Z<-(t(t(p[cumsum(weight),])*sqrt(attr(p,"coefs")$norm2[-seq(2)]))[,degree])
class(Z)<-"OPoly";Z
}
##this OPoly is an FORTRAN orthogonal polynomial routine, it first maps the
x to range[-2,2] then do QR, then return the results with sqrt(norm2).
Comparing with poly, this transformation will make the model coefficients
within a similar range as other variables, the R poly routine will usually
give you a very large coefficients. I did not find such routine in R, so I
have to define this as user defined function.
#######
I also have following function as you suggested:
makepredictcall.OPoly<-function(var,call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(OPoly))) {
if (!is.null(tmp <- attr(var, "coefs"))) {
call$coefs <- tmp
}
}
}
call
}
But I still got error for following:
> g3=glm(lot1 ~ log(u) + OPoly(u,1), data = clotting, family = Gamma)
> predict(g3,dc)Error in poly(4 * (rep(x, weight) -
mean(range(x)))/diff(range(x)), degree) :
missing values are not allowed in 'poly'
I thought it might be due to the /diff(range(x) in the function. But even
I remove that part, it will still give me error. Any idea?
Many thanks in advance.
Alex
On Thu, Jul 16, 2015 at 2:09 PM, William Dunlap <wdunlap at tibco.com>
wrote:
> Read about the 'makepredictcall' generic function. There is a
method,
> makepredictcall.poly(), for poly() that attaches the polynomial
> coefficients
> used during the fitting procedure to the call to poly() that predict()
> makes.
> You ought to supply a similar method for your xpoly(), and xpoly() needs
> to return an object of a a new class that will cause that method to be
used.
>
> E.g.,
>
> xpoly <- function(x,degree=1,...){ ret <- poly(x,degree=degree,...);
> class(ret) <- "xpoly" ; ret }
> makepredictcall.xpoly <- function (var, call)
> {
> if (is.call(call)) {
> if (identical(call[[1]], quote(xpoly))) {
> if (!is.null(tmp <- attr(var, "coefs"))) {
> call$coefs <- tmp
> }
> }
> }
> call
> }
>
> g2 <- glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
> predict(g2,dc)
> # 1 2 3 4
> 5
> #-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
> #-0.01398928608
> # 6 7 8 9
> #-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
>
> You can see the effects of makepredictcall() in the 'terms'
component
> of glm's output. The 'variables' attribute of it gives the
original
> function
> calls and the 'predvars' attribute gives the calls to be used for
> prediction:
> > attr(g2$terms, "variables")
> list(lot1, log(u), xpoly(u, 1))
> > attr(g2$terms, "predvars")
> list(lot1, log(u), xpoly(u, 1, coefs = list(alpha = 40, norm2 = c(1,
> 9, 8850))))
>
>
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
> On Thu, Jul 16, 2015 at 12:35 PM, Kunshan Yin <yinkunshan at
gmail.com>
> wrote:
>
>> Hello, I have a question about the formula and the user defined
function:
>>
>> I can do following:
>> ###Case 1:
>> > clotting <- data.frame(
>> + u = c(5,10,15,20,30,40,60,80,100),
>> + lot1 = c(118,58,42,35,27,25,21,19,18),
>> + lot2 = c(69,35,26,21,18,16,13,12,12))
>> > g1=glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma)
>> > dc=clotting
>> > dc$u=1
>> > predict(g1,dc)
>> 1 2 3 4 5
>> 6 7 8 9
>> -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929
>> -0.01398929 -0.01398929 -0.01398929
>>
>> However, if I just simply wrap the poly as a user defined function ( in
>> reality I would have my own more complex function) then I will get
error:
>> ###Case 2:
>> > xpoly<-function(x,degree=1){poly(x,degree)}
>> > g2=glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family =
Gamma)
>> > predict(g2,dc)
>> Error in poly(x, degree) :
>> 'degree' must be less than number of unique points
>>
>> It seems that the predict always treat the user defined function in the
>> formula with I(). My question is how can I get the results for Case2
>> same
>> as case1?
>>
>> Anyone can have any idea about this?
>>
>> Thank you very much.
>>
>> Alex
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
[[alternative HTML version deleted]]