[R] using "unstack" inside my function: that old scope problem again

Peter Dalgaard p.dalgaard at biostat.ku.dk
Fri Mar 19 10:09:39 CET 2004


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




More information about the R-help mailing list