[R] Coefficients in a polynomial glm with family poisson/binomial

[sam_oi at yahoo.fr]

Dear R-users

When running a glm polynomial model with one explanatory variable (example Y~X+X^2), with a poisson or binomial error distribution, the predicted values obtained from using the predict() function and those obtained from using the coefficients from the summary table "as is" in an equation of the form Y=INTERCEPT+ XCoef x X + XCoef x X^2, differ considerably. The former are correct and the latter are wrong. 
This does not occur using lm() or in a glm with family as normal. I conclude that this is due to the link function, predict() having some way of back transforming the data. But if this is so, are the estimated coefficients wortheless in this case? 
I need to get accurate coefficients (for use in another model using offset), and have resorted to re-estimating them by running a second polynomial (lm() this time) on the predicted values from predict() of the glm. This is clearly not a nice way of doing things. 

Could anyone please inform me of why this is happening and of a better way around this? 


Code:

glm2<-glm(FEDSTATUS1~AGE+I(AGE^2), family=binomial(link="probit"))
summary(glm2) ### first set of "wrong coefficients"

nd1<-expand.grid(AGE=c(1:70))
Pred.Fed1<-predict(glm2,nd1,type="response")
points(predict(glm2,nd1,type="response")~nd1$AGE, col=2)


AGE11<-c(11:70)
Pred<-t(rbind(Pred.Fed1,AGE11))
Pred<-as.data.frame(Pred)
model<-lm(Pred$Pred.Fed1~Pred$AGE11+I(Pred$AGE11^2))
summary(model) ### "accurate coefficients"


Thanks

Samuel Riou
University of Leeds