[R] accessing variables inside a function, inside a loop

Joshua Wiley jwiley.psych at gmail.com
Wed Mar 2 06:43:33 CET 2011


Hi Alan,

Other more knowledgeable people may have better opinions on this than
I do.  Manipulating language and call objects is seriously stretching
my skills.  \

In any case, two ways come to mind, both of them sufficiently
cumbersome I would seriously question the value (btw, this is a
completely different question, right?).  To borrow from Barry
Rowlingson, I'd like to prefix all these 'solutions' with "Here's how
to do it, but don't actually do it."  The first option would be to
manually construct a call to glm() and then evaluate it.  That is,
rather than pass "frm" to the formula argument, construct a text
string of the entire glm call.  Something like:

#####################################
foo <- function(y) {
  paste("glm(formula = ", y, " ~ hp * wt, data = mtcars)", sep = '')
}
foo("mpg")
eval(parse(text = foo("mpg")))
#####################################

The other thought would be to just update the part of the model object
containing the call (I actually like this better than my first option,
though I'm still not fond of it).  Assuming you do not actually need
the entire call, you could easily use the formula.  Here's an example:

#####################################
xReg <- function(y) {
  frm <- eval(substitute(p ~ hp * wt, list(p = as.name(y))))
  mod <- glm(frm, data = mtcars)
  mod$call <- frm
  return(mod)
}

xReg("mpg")
#####################################


If you want to know what the formula used for a model is, my
suggestion would be to simply have your function xReg() return both
the model object AND the formula you used (i.e., "frm").  Here is an
example:

#####################################
## *my* preference

xReg <- function(y) {
  frm <- eval(substitute(p ~ hp * wt, list(p = as.name(y))))
  mod <- glm(frm, data = mtcars)
  output <- list(formula = frm, model = mod)
  attributes(output$formula) <- NULL
  return(output)
}

xReg("mpg")
#####################################

Side note, Dr. Bates (author of lme4, genius, and a nice, helpful
person to boot) taught me how to use substitute() for something I
tried once on the ggplot2 list.

Cheers,

Josh

On Tue, Mar 1, 2011 at 7:04 PM, zListserv <zlistserv at gmail.com> wrote:
> Joshua
>
> Great solution.  Taking off on your code, the following works but does not display the names of the variables in the formula.  Any suggestions about how to modify the function so that it displays the correct formula (e.g., "glm(formula = y1 ~ x1 * x2, data = dat)" instead of "glm(formula = frm, data = dat)")?
>
> R> x = runif(2000)
> R> y = runif(2000)
> R> z = runif(2000)
> R>
> R> y1 = y * x
> R> y2 = y * sqrt(x)
> R>
> R> x1 = y1 / y2 + z
> R> x2 = y2 / y1 * z + z
> R>
> R> dat = data.frame(y1,y2,x1,x2)
> R>
> R> xReg = function(y) {
> +
> +       frm = eval(substitute(p ~ x1 * x2, list(p = as.name(y))))
> +       mod = glm(frm, data=dat)
> +       }
> R>
> R> lapply(names(dat[,1:2]), xReg)
> [[1]]
>
> Call:  glm(formula = frm, data = dat)
>
> Coefficients:
> (Intercept)           x1           x2        x1:x2
> -0.1882452    0.4932059    0.0667401   -0.1310084
>
> Degrees of Freedom: 1999 Total (i.e. Null);  1996 Residual
> Null Deviance:      99.15032
> Residual Deviance: 67.71775     AIC: -1085.354
>
> [[2]]
>
> Call:  glm(formula = frm, data = dat)
>
> Coefficients:
> (Intercept)            x1            x2         x1:x2
> -0.005464627   0.386937367   0.037363416  -0.094136334
>
> Degrees of Freedom: 1999 Total (i.e. Null);  1996 Residual
> Null Deviance:      112.7078
> Residual Deviance: 90.24796     AIC: -510.9287
>
> ---
>
> Thanks,
>
> Alan
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/



More information about the R-help mailing list