[Rd] request for patch in "drop1" (add.R)

[Ben Bolker]

On 11-02-23 06:12 PM, Prof Brian Ripley wrote:
> residuals() and $residuals are often very different:  residuals() is
> generic, but even the default method is *not* simple extraction. Their
> values can be of different lengths: think about an lm fit with na.action
> = na.exclude.  That is precisely the sort of thing the tests in add.R
> were designed to detect, hence the use of $residuals.
> 
> None of this is used in drop1()!  One of the places $residuals was used
> in that file was the default method for drop1(), others being step() and
> the default method for add1().  As default methods these have to
> continue to work for any class of object that has previously been thrown
> at them over the last 10+ years, and even all CRAN packages will not be
> a good enough test suite.
> 
> In any case, the current code in R-devel makes use of the new generic
> nobs(), which is intended to help sort out the many versions of
> functions called 'BIC" in packages, but is also useful here.  (It is
> still under development.)
> 
> terms() is also generic so there is also some danger that its
> substitution could give an inappropriate result.  But as it used in
> several other places in add.R the breakage will probably occur elsewhere
> already.

  Thanks Prof. Ripley.
  (I will say that, while I understand why residuals(x) and x$residuals
could be different, I am happy that a more transparent form of coding is
being introduced ...)
> 
> On Wed, 23 Feb 2011, Martin Maechler wrote:
> 
>>>>>>> Ben Bolker <bbolker at gmail.com>
>>>>>>>     on Wed, 23 Feb 2011 09:14:37 -0500 writes:
>>
>>    >   By changing three lines in drop1 from access based on $
>>    > to access based on standard accessor methods (terms() and
>>    > residuals()), it becomes *much* easier to extend drop1 to
>>    > work with other model types.  The use of $ rather than
>>    > accessors in this context seems to be an oversight rather
>>    > than a design decision, but maybe someone knows better ...
>>
>>    >   In particular, if one makes these changes (which I am
>>    > pretty certain will not break anything, as the original
>>    > code basically mimicked the default methods anyway), it
>>    > becomes possible to make drop1() work with mer objects
>>    > (Doug Bates's new mixed model code) merely by defining:
>>
>>    > terms.mer <- function(x, ...) {
>>    > attr(x at frame,"terms")
>>    > }
>>
>>    > extractAIC.default <- function(fit, scale=0, k=2, ...) {
>>    > L <- logLik(fit)
>>    > edf <- attr(L,"df")
>>    > c(edf,-2*L+2*edf)
>>    > }
>>
>>    > Adding this definition of extractAIC.default also makes drop1() work
>>    > with lme fits ...
>>
>>    > Comments?  Should I submit to the bug database as "enhancement
>>    > request"?  Are there any hidden downsides to this?
>>
>> drawback: a possible very small performance cut for the cases
>> where the "$" access is ok.  But that should not count.
>>
>> I like the idea.... {it's a pity that only S3 methods work that way,
>> because residuals(), terms(), etc... are unfortunately not
>> general (implicit) S4 generics but just S3 ones..
>>
>> I'm currently testing your change, plus some more in step().
>> However, for step() to work automagically there is more needed.
>> It currently relies in more places on 'object' being a list to
>> which you can append new components, basically by
>>   fit <- object
>>   fit$new1 <- ...
>>   fit$new2 <- ...
>>
>> That would have to be changed to something like
>>   fit <- list(obj = object)
>>   fit$new1 <- ...
>>   fit$new2 <- ...
>> and more changes where 'fit' has to be replaced by 'fit$obj'.
>> Would that not be of interest as well?
>>
>> Martin
>>
>>
>>    > Ben Bolker
>>
>>    >
>> ----------------------------------------------------------------------
>>    > Index: add.R
>>    > ===================================================================
>>    > --- add.R    (revision 54562)
>>    > +++ add.R    (working copy)
>>    > @@ -330,7 +330,7 @@
>>    > drop1.default <- function(object, scope, scale = 0,
>> test=c("none", "Chisq"),
>>    > k = 2, trace = FALSE, ...)
>>    > {
>>    > -    tl <- attr(object$terms, "term.labels")
>>    > +    tl <- attr(terms(object), "term.labels")
>>    > if(missing(scope)) scope <- drop.scope(object)
>>    > else {
>>    > if(!is.character(scope))
>>    > @@ -344,7 +344,7 @@
>>    > ans <- matrix(nrow = ns + 1L, ncol = 2L,
>>    > dimnames =  list(c("<none>", scope), c("df", "AIC")))
>>    > ans[1, ] <- extractAIC(object, scale, k = k, ...)
>>    > -    n0 <- length(object$residuals)
>>    > +    n0 <- length(residuals(object))
>>    > env <- environment(formula(object))
>>    > for(i in seq(ns)) {
>>    > tt <- scope[i]
>>    > @@ -356,7 +356,7 @@
>>    > evaluate = FALSE)
>>    > nfit <- eval(nfit, envir=env) # was  eval.parent(nfit)
>>    > ans[i+1, ] <- extractAIC(nfit, scale, k = k, ...)
>>    > -        if(length(nfit$residuals) != n0)
>>    > +        if(length(residuals(nfit)) != n0)
>>    > stop("number of rows in use has changed: remove missing values?")
>>    > }
>>    > dfs <- ans[1L , 1L] - ans[, 1L]
>>
>>    >
>> ----------------------------------------------------------------------
>>    > ______________________________________________
>>    > R-devel at r-project.org mailing list
>>    > https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>