On Thu, 23 Apr 2015, Therneau, Terry M., Ph.D. wrote:
This issue has arisen within my anova.coxph routine, but is as easily
illustrated with glm.
testdata <- data.frame(y= 1:5,
n= c(8,10,6,20,14),
sex = c(0,1,0,1,1),
age = c(30,20,35,25,40))
fit <- glm(cbind(y,n) ~ age + sex, binomial, data=testdata, model=TRUE)
saveit <- fit$model
update(fit, .~. - age, data=saveit)
Error in cbind(y, n) : object 'y' not found
One would hope that a saved model frame is precisely the thing that
would work best. The issue of course is that "cbind(y, n)" is the name
of the first variable in saveit, and it is not being properly quoted
somewhere down the line. The same issue can occur on the right hand
side. "Save the model frame in case you need to refit something next
month" is does not appear to be a safe approach to reproducable
research.
If you want to re-run the same call (which is what the default update
method does), then you need to save the original data not the
pre-processed model frame. However, the model frame still has all the
information you need - but has already evaluated all transformations
(cbind, log, etc.).
fit2 <- glm(y ~ sex + log(age), poisson, testdata)
save2 <- fit2$model
update(fit2, . ~ . - sex, data=save2) # fails
glm(y ~ log(age), poisson, save2) # fails
I can work around this in my anova, but I wanted to not rebuild the
frame if it is already present. It looks like model.matrix plus attr(x,
'assign') time -- a bit harder to read, but that looks like what
anova.glm is doing. Is there a way to make update work?
If you have the model.frame() plus an update formula you could internally
extract the response, model matrix and weights, e.g.,
## original fit
fit <- glm(cbind(y,n) ~ log(age) + sex, family = binomial,
data = testdata)
## original model frame
mf <- fit$model
## update formula
up <- . ~ . - sex
## new terms
mt <- update(terms(mf), up)
## response y, matrix x, weights w (NULL here)
y <- model.response(mf)
x <- model.matrix(mt, mf)
w <- model.weights(mf)
## offset could be added similarly
And then you can call glm.fit() or coxph.fit() if you can get the family
from the original fit.
One should still check that the update formula does not change the
response (to something which may or may not exist in the model frame).
Possibly, one could try to get certain control arguments from the original
glm fit (but probably not 'start').
Maybe this helps... (But I can understand that the issues of data frame
vs. model frame can be quite a nuisance when programming model utilities
:))
The current code, BTW, starts by building its own frame using results of
terms.inner, which solves the above issue nicely and update() works as
expected. But it isn't robust to scoping issues. (As pointed out
yesterday by a user: lapply of a function that contained coxph followed
by anova gives a variable not found error.) It also ignores saved model
frames; thus the rewrite.
Terry T
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
