R help - Mar 2004 - using "unstack" inside my function: that old scope problem again

I've been reading the R mail archives and I've found  a lot of messages 
with this same kind of problem, but I can't understand the answers.  Can 
one of you try to explain this to me?

Here's my example. Given a regression model and a variable, I want to 
use unstack() on the vector of residuals and make some magic with the 
result. But unstack hates me.
 

PCSE <- function (tmodel,groupVar) {
  myres1 <- resid(tmodel)
  resUnstacked <- unstack(myres1, form = myres1 ~ groupVar));
  E <- as.matrix(resUnstacked)
  SIGMA <- (1/nrow(E))*(t(E) %*% E)
  OMEGA <- diag(x=1, nrow=nrow(E), ncol=nrow(E)) %x% SIGMA

  X <- model.matrix(tmodel)
  XPRIMEXINV <- solve(t(X) %*% X)
  PCSECOVB <- XPRIMEXINV %*%  (t(X) %*% OMEGA %*% X ) %*% XPRIMEXINV
}


The error is:
PCSE(eld.ols1,dat2$STATEID)
Error in eval(expr, envir, enclos) : Object "groupVar" not found

Here's what I don't understand the most.
 
If I hack this so that the "resUnstacked" is created by a
"matrix"
command, then the function works. Why would matrix() work when unstack 
does not.  And why does model.matrix() work if unstack does not.

Thanks in advance, as usual.

-- 
Paul E. Johnson                       email: pauljohn at ku.edu
Dept. of Political Science            http://lark.cc.ku.edu/~pauljohn
1541 Lilac Lane, Rm 504                              
University of Kansas                  Office: (785) 864-9086
Lawrence, Kansas 66044-3177           FAX: (785) 864-5700

Paul Johnson <pauljohn at ku.edu> writes:
> I've been reading the R mail archives and I've found  a lot of
> messages with this same kind of problem, but I can't understand the
> answers.  Can one of you try to explain this to me?
> 
> Here's my example. Given a regression model and a variable, I want to
> use unstack() on the vector of residuals and make some magic with the
> result. But unstack hates me.
>  PCSE <- function (tmodel,groupVar) {
>   myres1 <- resid(tmodel)
>   resUnstacked <- unstack(myres1, form = myres1 ~ groupVar));
>   E <- as.matrix(resUnstacked)
>   SIGMA <- (1/nrow(E))*(t(E) %*% E)
>   OMEGA <- diag(x=1, nrow=nrow(E), ncol=nrow(E)) %x% SIGMA
> 
>   X <- model.matrix(tmodel)
>   XPRIMEXINV <- solve(t(X) %*% X)
>   PCSECOVB <- XPRIMEXINV %*%  (t(X) %*% OMEGA %*% X ) %*% XPRIMEXINV
> }
> 
> 
> The error is:
> PCSE(eld.ols1,dat2$STATEID)
> Error in eval(expr, envir, enclos) : Object "groupVar" not found
> 
> Here's what I don't understand the most.
>  If I hack this so that the "resUnstacked" is created by a
"matrix"
> command, then the function works. Why would matrix() work when unstack
> does not.  And why does model.matrix() work if unstack does not.
> 
> Thanks in advance, as usual.
Look inside getAnywhere(unstack.default) and you'll find things like

x <- as.list(x)
...
res <- c(tapply(eval(form[[2]], x), eval(form[[3]], x), as.vector))


Now, "x" is myres1, form[[2]] is quote(myres1), and form[[3]] is
quote(groupVar). myres1 would appear to be a vector and I suspect it
doesn't have any elements of that name...

You might be able to unstack(data.frame(myres1,groupVar), ..etc..)

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907

The problem is the use of a formula, and hence non-standard evaluation 
(even non-standard for a formula).  Unstack contains

    res <- c(tapply(eval(form[[2]], x), eval(form[[3]], x), as.vector))

and those eval() calls are looking only in x and then along the search
path.  They should be something like

eval(form[[2]], x, parent.frame())

or

eval(form[[2]], x, environment(form))

depending on your preferences for scoping rules.

unstack.default is very simple: why not extract the part you want?


On Fri, 19 Mar 2004, Paul Johnson wrote:
> I've been reading the R mail archives and I've found  a lot of
messages
> with this same kind of problem, but I can't understand the answers. 
Can
> one of you try to explain this to me?
> 
> Here's my example. Given a regression model and a variable, I want to 
> use unstack() on the vector of residuals and make some magic with the 
> result. But unstack hates me.
>  
> 
> PCSE <- function (tmodel,groupVar) {
>   myres1 <- resid(tmodel)
>   resUnstacked <- unstack(myres1, form = myres1 ~ groupVar));
>   E <- as.matrix(resUnstacked)
>   SIGMA <- (1/nrow(E))*(t(E) %*% E)
>   OMEGA <- diag(x=1, nrow=nrow(E), ncol=nrow(E)) %x% SIGMA
> 
>   X <- model.matrix(tmodel)
>   XPRIMEXINV <- solve(t(X) %*% X)
>   PCSECOVB <- XPRIMEXINV %*%  (t(X) %*% OMEGA %*% X ) %*% XPRIMEXINV
> }
> 
> 
> The error is:
> PCSE(eld.ols1,dat2$STATEID)
> Error in eval(expr, envir, enclos) : Object "groupVar" not found
> 
> Here's what I don't understand the most.
>  
> If I hack this so that the "resUnstacked" is created by a
"matrix"
> command, then the function works. Why would matrix() work when unstack 
> does not.  And why does model.matrix() work if unstack does not.
Because the eval calls in model.matrix.default have been thought through.
> 
> Thanks in advance, as usual.
> 
> 
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595