[Rd] model frames and update()

[Achim Zeileis]

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 at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>