[R] scoping problem

Duncan Murdoch murdoch at stats.uwo.ca
Wed Oct 24 22:09:44 CEST 2007

On 10/24/2007 3:43 PM, Sandy Weisberg wrote:
> I would like to write a function that computes Tukey's 1 df for 
> nonadditivity.  Here is a simplified version of the function I'd like to 
> write:  (m is an object created by lm):
> tukey.test <- function(m) {
>   m1 <- update(m, ~.+I(predict(m)^2))
>   summary(m1)$coef
>   }
> The t-test for the added variable is Tukey's test.  This won't work:
> data(BOD)
> m1 <- lm(demand~Time,BOD)
> tukey.test(m1)
> Error in predict(m) : object "m" not found
> This function doesn't work for two reasons:
> 	1. The statement m1 <- update(m, ~.+I(predict(m)^2)) can't see 'm' in 
> the call to predict.

It is looking in the environment of m$terms, and failing to find m 
there.  Here's an ugly way to do what you want; there might be a cleaner 
way, but I don't know it:

tukey.test <- function(m) {
   # Create a new environment whose parent is where update is currently
   # looking
   newenv <- new.env(parent=environment(m$terms))

   # Put the new predictor there.  You could instead put m there, but
   # that will create a loop whose evaluation will destroy the universe
   newenv$predm2 <- predict(m)^2

   # Attach the new environment in place of the old one
   environment(m$terms) <- newenv

   # Proceed as before
   m1 <- update(m, ~.+predm2)

Another variation is just to assign predm2 into the existing 
environment, i.e.

   environment(m$terms)$predm2 <- predict(m)^2

but that's potentially destructive of an existing variable named predm2 
so I didn't do it.  Environments are reference objects, so you would 
kill the real one, and have side effects outside your scope:  bad idea.

> 	2. If in creating m missing values had been present, then predict(m), 
> even if it could be computed, could be of the wrong length.

Not sure about the best approach here, but the prediction vector has 
names indicating which observations were used, so you can probably do 
something with those.

Duncan Murdoch


More information about the R-help mailing list