another rstanarm bug, this time with k-fold cross-validation, kfold()

[Andrew Gelman]

So. I wanted to crate a model where loo would fail. I did so using fake-data simulation, creating 5 random predictors that were correlated 0.8 with each other:

library("rstanarm")
options(mc.cores = parallel::detectCores())
library("arm")
library("MASS") # needed for mvrnorm()

n <- 200
k <- 5
rho <- 0.8
Sigma <- rho*array(1, c(k,k)) + (1-rho)*diag(k)
X <- mvrnorm(n, rep(0,k), Sigma)
b <- -2 + (1:k)
sigma <- 2
z <- rlogis(n, X %*% b, 1)
y <- ifelse(z>0, 1, 0)

## First I fit the full model:
fit_1 <- stan_glm(y ~ X, family=binomial(link="logit"), prior=normal(0, 10, autoscale=FALSE))
loo_1 <- loo(fit_1)
## I got this warning:

Warning message:
Found 23 observations with a pareto_k > 0.7. With this many problematic observations we recommend calling 'kfold' with argument 'K=10' to perform 10-fold cross-validation rather than LOO.

## OK, fine. So I do 10-fold cross-validation:
kfold_1 <- kfold(fit_1, K=10)
print(kfold_1)
## And this works just fine.

## Now I want to compare to a model removing predictor 1
Z <- X[,2:5]
fit_2 <- stan_glm(y ~ Z, family=binomial(link="logit"), prior=normal(0, 10, autoscale=FALSE))
kfold_2 <- kfold(fit_2, K=10)
print(kfold_2)
## This works fine too.

## Then I try putting X[,2:5] right into the model (instead of giving it the new name Z):
fit_2 <- stan_glm(y ~ X[,2:5], family=binomial(link="logit"), prior=normal(0, 10, autoscale=FALSE))
## This works fine.
## But then I try the next step:
kfold_2 <- kfold(fit_2, K=10)
## And I get an error message:

Fitting model 1 out of 10
Error in X[, 2:5] : incorrect number of dimensions

## Whassup?? The model fit fine. So why the problem with kfold()?

A

[Jonah Gabry]

So... kfold attempts to account for the fact that sometimes data is specified using the data argument to stan_glm and other times using the global environment. When no data argument is specified it needs to construct a data frame that it can use to refit the model K times (each time modifying that data frame to drop the right set of observation). So it constructs a data frame consisting of the 2nd through 5th columns of X because that's what you implicitly told it was your data. But now when we need to go refit the model how is the formula including X[, 2:5] supposed to be interpreted relative to this data frame? There are ways around this but so far they seem to introduce other complications that are even more annoying. 

[Jonah Gabry]

Also, it's really really really hard to write functions like kfold that handle all possible things a user can do when specifying a model. This is the curse of R's flexibility. It's nice to be able to write the same model in many ways but that also means it's a nightmare to make sure all of those specifications are compatible with every thing you might want to do after the fact. It would be much easier if we allowed fewer ways of specifying the same model. 

[Ben Goodrich]

On Wednesday, February 1, 2017 at 11:44:46 AM UTC-5, Jonah Gabry wrote:

There are ways around this but so far they seem to introduce other complications that are even more annoying. 

In general, it is better to follow the advice of Hadley and many others that the variables specified in a formula passed to a model fitting function in the rstanarm package should be contained in a data.frame that is passed to the data argument.

Ben
 

[Jonah Gabry]

If we required that then we could save ourselves a lot of problems. But if we want to mimic glm and similar functions we need to keep accepting formulas referring to objects in the global environment.  

[Jonah Gabry]

I would LOVE to require that the user specify the data argument. Objections?

[Bob Carpenter]

I really don't like this "read out of the global environment"
thing. I hope Andrew's not writing the new book to lean heavily
on that.

- Bob

> On Feb 1, 2017, at 11:44 AM, Jonah Gabry <jga...@gmail.com> wrote:
>
> So... kfold attempts to account for the fact that sometimes data is specified using the data argument to stan_glm and other times using the global environment. When no data argument is specified it needs to construct a data frame that it can use to refit the model K times (each time modifying that data frame to drop the right set of observation). So it constructs a data frame consisting of the 2nd through 5th columns of X because that's what you implicitly told it was your data. But now when we need to go refit the model how is the formula including X[, 2:5] supposed to be interpreted relative to this data frame? There are ways around this but so far they seem to introduce other complications that are even more annoying.
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Jonah Sol Gabry]

On Wed, Feb 1, 2017 at 1:17 PM, Bob Carpenter <ca...@alias-i.com> wrote:

I really don't like this "read out of the global environment"
thing.  I hope Andrew's not writing the new book to lean heavily
on that.

- Bob

Me too!

[Andrew Gelman]

I'm only including code that runs!

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls.
A

[Bob Carpenter]

One man's crap is another man's fertilizer.

- Bob

[Jonah Gabry]

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs! 

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls. 
A 

Well, not "data=data" but maybe "data=wells" isn't so bad? Unlike the list you need to (or needed to) manually create for rstan, for rstanarm it would just be the name of the data frame where the variables in the formula can be found. I understand the appeal of not needing to do that, but I think that's outweighed by the mess it causes for software development. Here's a pros/cons list:

Cons of searching for data in global (or local) environment if not specified:

1) Takes developers a lot more time to add anything new to rstanarm because it's not easy to ensure that the data is handled appropriately. This is not an exaggeration: at most 1/4 of the time I spent implementing kfold was getting the algorithm right and the rest was spent on avoiding errors when people specify data without a data argument. And even after all that time there are apparently still cases when it doesn't work!) 

2) Inevitably more bugs because the internal code needs to be much more complicated and is harder to test

3) Easier for the user to make mistakes (need to be very vigilant about what's already in the global environment) 

4) I could probably come up with a bunch more!

Pros of searching for data in global (or local) environment if not specified:

1) A little bit less typing for the user

 

Am I missing any other advantages other than a bit less typing? I really don't want to sacrifice the quality of the software just to avoid typing something like data=wells, but if there are other advantages I'm missing then maybe those will tip the scale. 

Jonah

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs!

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls.
A

> On Feb 1, 2017, at 1:17 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> I really don't like this "read out of the global environment"
> thing.  I hope Andrew's not writing the new book to lean heavily
> on that.  
>
> - Bob
>
>> On Feb 1, 2017, at 11:44 AM, Jonah Gabry <jga...@gmail.com> wrote:
>>
>> So... kfold attempts to account for the fact that sometimes data is specified using the data argument to stan_glm and other times using the global environment. When no data argument is specified it needs to construct a data frame that it can use to refit the model K times (each time modifying that data frame to drop the right set of observation). So it constructs a data frame consisting of the 2nd through 5th columns of X because that's what you implicitly told it was your data. But now when we need to go refit the model how is the formula including X[, 2:5] supposed to be interpreted relative to this data frame? There are ways around this but so far they seem to introduce other complications that are even more annoying.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Andrew Gelman]

Hi, Jonah.  We can talk about it on Friday.  Short answer is that when working with data and fitting models it's a real pain in the ass to keep working with a data frame.  I actually don't like the whole "data frame" thing as it can act as a straitjacket, it's a Stata way of thinking.

Here's an example.  I start with this model (for simplicity I'm assuming we have a stan_logit() alias):

stan_logit(y ~ dist + arsenic)

or

stan_logit(y ~ dist + arsenic, data="wells")

Ok, fine so far.  Now I decide I want to code as dist/100:

dist100 <- dist/100

stan_logit(y ~ dist100 + arsenic)

That worked fine.  Now the alternative:

dist100 <- dist/100 INSIDE "wells"

stan_logit(y ~ dist100 + arsenic, data="wells")

So here's my question:  can I do this:

dist100 <- dist/100 INSIDE "wells"

in a reasonable way?

Right now I have to first pull dist out of wells, then compute dist100, then put it back into wells, something like this:

dist <- wells$dist

dist100 <- dist/100

APPEND wells with dist100

(I don't remember how to do that last command but I know it can be done in one line.)

Anyway, that's the point.  In the midst of running regressions I play with data, for example creating a new predictor:

x_avg <- (x1 + x2 + x3)/3

or whatever.

If I have to keep putting things into the data frame just to run my regressions, that's a pain in the ass.  It's _not_ just "a bit more typing," it's actual additional lines of code, it makes things harder to follow.

In Stata world, you have a big data matrix and then pull out columns of it.  I don't see this as being R world.

A

On Feb 1, 2017, at 4:47 PM, Jonah Gabry <jga...@gmail.com> wrote:

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs! 

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls. 
A 

Well, not "data=data" but maybe "data=wells" isn't so bad? Unlike the list you need to (or needed to) manually create for rstan, for rstanarm it would just be the name of the data frame where the variables in the formula can be found. I understand the appeal of not needing to do that, but I think that's outweighed by the mess it causes for software development. Here's a pros/cons list:

Cons of searching for data in global (or local) environment if not specified:

1) Takes developers a lot more time to add anything new to rstanarm because it's not easy to ensure that the data is handled appropriately. This is not an exaggeration: at most 1/4 of the time I spent implementing kfold was getting the algorithm right and the rest was spent on avoiding errors when people specify data without a data argument. And even after all that time there are apparently still cases when it doesn't work!) 

2) Inevitably more bugs because the internal code needs to be much more complicated and is harder to test

3) Easier for the user to make mistakes (need to be very vigilant about what's already in the global environment) 

4) I could probably come up with a bunch more!

Pros of searching for data in global (or local) environment if not specified:

1) A little bit less typing for the user

 

Am I missing any other advantages other than a bit less typing? I really don't want to sacrifice the quality of the software just to avoid typing something like data=wells, but if there are other advantages I'm missing then maybe those well tip the scale. 

Jonah

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs!

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls.
A

> On Feb 1, 2017, at 1:17 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> I really don't like this "read out of the global environment"
> thing.  I hope Andrew's not writing the new book to lean heavily
> on that.  
>
> - Bob
>
>> On Feb 1, 2017, at 11:44 AM, Jonah Gabry <jga...@gmail.com> wrote:
>>
>> So... kfold attempts to account for the fact that sometimes data is specified using the data argument to stan_glm and other times using the global environment. When no data argument is specified it needs to construct a data frame that it can use to refit the model K times (each time modifying that data frame to drop the right set of observation). So it constructs a data frame consisting of the 2nd through 5th columns of X because that's what you implicitly told it was your data. But now when we need to go refit the model how is the formula including X[, 2:5] supposed to be interpreted relative to this data frame? There are ways around this but so far they seem to introduce other complications that are even more annoying.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[FHRB Toledo]

Dear All,

I am following this thread and in some sense agree with Jonah.

BTW, regarding the question from Professor Gelman:

> So here's my question:  can I do this:

> dist100 <- dist/100 INSIDE "wells"

> in a reasonable way?

I think you can do:

stan_logit(y ~ I(dist / 100) + arsenic, data="wells")

Where I() stands for: https://stat.ethz.ch/R-manual/R-devel/library/base/html/AsIs.html

Cheers,

FH

[Ben Goodrich]

On Wednesday, February 1, 2017 at 5:23:20 PM UTC-5, Andrew Gelman wrote:

So here's my question:  can I do this:

dist100 <- dist/100 INSIDE "wells"

in a reasonable way?

You can do that in several ways, such as

wells$dist100 <- wells$dist / 100

But it is more readable to do that inside the formula as FH suggests.

Ben

[Andrew Gelman]

Yes, I do know about that, but I like to just define intermediate quantities and go from there.  I find I() a bit awkward.

Also I think that things that can be done in lm() and glm() should be doable in stan_lm() and stan_glm().

That said, I recoginse that there are tradeoffs.

A

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Andrew Gelman]

Oh, I didn't know about this:  

wells$dist100 <- wells$dist / 100

Intersting.  Still a bit awkward, but, hmm, I'll have to thikn about this.

A

[Michael Betancourt]

+1.  People are going to copy and paste code from the book,

and if there’s no requirement of passing things through a data

frame then the pasted code will run even if the environment is

not correctly set up and result in bad fits.  Requiring a data=

argument, however, will fail immediately and guide the user

to set up their environment appropriately.  

The most annoying problems I’ve seen in workshops and

tutorials is exactly when the code tries to be smart and do

things automatically.  Easy debugging >>>> slightly more

compact code.

On Feb 1, 2017, at 4:47 PM, Jonah Gabry <jga...@gmail.com> wrote:

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs! 

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls. 
A 

Well, not "data=data" but maybe "data=wells" isn't so bad? Unlike the list you need to (or needed to) manually create for rstan, for rstanarm it would just be the name of the data frame where the variables in the formula can be found. I understand the appeal of not needing to do that, but I think that's outweighed by the mess it causes for software development. Here's a pros/cons list:

Cons of searching for data in global (or local) environment if not specified:

1) Takes developers a lot more time to add anything new to rstanarm because it's not easy to ensure that the data is handled appropriately. This is not an exaggeration: at most 1/4 of the time I spent implementing kfold was getting the algorithm right and the rest was spent on avoiding errors when people specify data without a data argument. And even after all that time there are apparently still cases when it doesn't work!) 

2) Inevitably more bugs because the internal code needs to be much more complicated and is harder to test

3) Easier for the user to make mistakes (need to be very vigilant about what's already in the global environment) 

4) I could probably come up with a bunch more!

Pros of searching for data in global (or local) environment if not specified:

1) A little bit less typing for the user

 

Am I missing any other advantages other than a bit less typing? I really don't want to sacrifice the quality of the software just to avoid typing something like data=wells, but if there are other advantages I'm missing then maybe those well tip the scale. 

Jonah

On Wednesday, February 1, 2017 at 4:06:31 PM UTC-5, Andrew Gelman wrote:

I'm only including code that runs!

But, that's right, I don't want to include crap like "data=data" in all of my stan and stan_glm calls.
A

> On Feb 1, 2017, at 1:17 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> I really don't like this "read out of the global environment"
> thing.  I hope Andrew's not writing the new book to lean heavily
> on that.  
>
> - Bob
>
>> On Feb 1, 2017, at 11:44 AM, Jonah Gabry <jga...@gmail.com> wrote:
>>
>> So... kfold attempts to account for the fact that sometimes data is specified using the data argument to stan_glm and other times using the global environment. When no data argument is specified it needs to construct a data frame that it can use to refit the model K times (each time modifying that data frame to drop the right set of observation). So it constructs a data frame consisting of the 2nd through 5th columns of X because that's what you implicitly told it was your data. But now when we need to go refit the model how is the formula including X[, 2:5] supposed to be interpreted relative to this data frame? There are ways around this but so far they seem to introduce other complications that are even more annoying.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Jonah Sol Gabry]

Ok we'll talk more Friday about this but just a few quick things before then:

On Wed, Feb 1, 2017 at 5:23 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Also I think that things that can be done in lm() and glm() should be doable in stan_lm() and stan_glm().

I don't fully trust that this always works properly in lm and glm either. I believe that the lm and glm functions fit the models correctly but I have no doubt that there are methods written for those fitted model objects somewhere with subtle bugs because of this issue.

So here's my question:  can I do this:

dist100 <- dist/100 INSIDE "wells"

in a reasonable way?

Without having to load extra packages (e.g. Hadley's dplyr) you can use any of the following:

wells$dist100 <- wells$dist / 100

wells$dist100 <- with(wells, dist / 100)

wells <- within(wells, dist100 <- dist / 100)

Not super elegant, I admit, but there are options. 

If I have to keep putting things into the data frame just to run my regressions, that's a pain in the ass.  It's _not_ just "a bit more typing," it's actual additional lines of code, it makes things harder to follow.

Ok yeah that can definitely be a pain in the ass. There's probably some middle ground we can find when we discuss on Friday.

Jonah

P.S. I'm reminded of your blog post about how frustrating it was when you found there was a dataset "attitude" automatically loaded into your R session. That's similar to the feeling everyone who writes software packages has whenever we think about searching for variables in the global environment. You just never quite know what you're going to find!

[Jonah Gabry]

There's also the fact that there will be code included in the book and that requires some special considerations. In particular, I think Michael makes a good point when he says in his comment above:

[Andrew Gelman]

I really don't like this idea that I'm working "out of the global environment." I'd like to do everything within the local environment. That's what's standard in R, and I just want to do standard R. For example, if you're working in a function in R, you use the local environment.

I feel like you (Bob) have the impression that I want to use global variables (and maybe GOTO statements). But that's not it at all. In R there's always a well-defined local environment (or, at least, I think there is), and that's what I think should be used. NOT the global environment. I never wanted that.

A

> On Feb 1, 2017, at 1:17 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>

[Jonah Sol Gabry]

We've been (at least Ive been) a bit sloppy with terminology and I think that's leading to some misunderstanding. Some of it has to do with the global environment but some of it is still a problem with local environments too. 

The real issue is that R will look for variables in parent environments (and then keep looking if it doesn't find them).  Often the global environment is the only one but it can also include local environments that themselves are parents to other environments and this is often just as problematic. It's really useful for some things but anything that isn't taken directly from user specified input can lead to trouble when it needs to be used later. Any specification of data that isn't directly input by the user is messy even if the data exists in an intermediate environment between the modeling function and the global environment. So even behavior that is natural to R is actually problematic for development of modeling software  because the same problem can occur no matter if the sole parent environment is the global environment or if there are layers of environments in between.

Ultimately it's really just that it's much harder to make robust software when you have to deal with R's handling of environments (not just global) instead of relying on the user to just tell you exactly what to use. 

I definitely see why being stricter about this can be a pain from a user perspective, so hopefully we can find a good compromise. 

Jonah

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/gwaEpERn1Io/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Mike Lawrence]

On Wed, Feb 1, 2017 at 6:42 PM Andrew Gelman <gel...@stat.columbia.edu> wrote:

Oh, I didn't know about this:  

wells$dist100 <- wells$dist / 100

Oh my! Can someone please get Hadley and Andrew in a room? 

In the interim, Andrew, you might enjoy this: http://r4ds.had.co.nz/pipes.html

--

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Bob Carpenter]

> On Feb 1, 2017, at 5:23 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>
> Hi, Jonah. We can talk about it on Friday. Short answer is that when working with data and fitting models it's a real pain in the ass to keep working with a data frame. I actually don't like the whole "data frame" thing as it can act as a straitjacket, it's a Stata way of thinking.

We're not talking about a data frame here, but
a list, right? Don't data frames have to be square?

I have no idea what Stata does.

> Here's an example. I start with this model (for simplicity I'm assuming we have a stan_logit() alias):
>
> stan_logit(y ~ dist + arsenic)
> or
> stan_logit(y ~ dist + arsenic, data="wells")

Are you passing the name of a list or something here?

> Ok, fine so far. Now I decide I want to code as dist/100:
> dist100 <- dist/100
> stan_logit(y ~ dist100 + arsenic)

This is where I think it's super helpful for those *reading*
your work (which includes you in the future) to keep
everything together. At this point, you have defined y somewere
before the last fit and define a new dist100 here and everything's
just sitting in the global environment ready to be modified.

I'd prefer to keep the variables in a list. At least in
any code I wanted to show to anyone.

> That worked fine. Now the alternative:
> dist100 <- dist/100 INSIDE "wells"
> stan_logit(y ~ dist100 + arsenic, data="wells")
>
> So here's my question: can I do this:
> dist100 <- dist/100 INSIDE "wells"
> in a reasonable way?
>
> Right now I have to first pull dist out of wells, then compute dist100, then put it back into wells, something like this:
> dist <- wells$dist
> dist100 <- dist/100
> APPEND wells with dist100
> (I don't remember how to do that last command but I know it can be done in one line.)
>
> Anyway, that's the point. In the midst of running regressions I play with data,

Right. You define

dist100 <- dist/100

Which is playing with data. And now the definition of the data
you're passing into the regression is spread through your history.

What Michael was arguing for is keeping the definition
of the data for a model you're going to fit encapsulated
in one little block of code that produces a list as
a result.

> for example creating a new predictor:
> x_avg <- (x1 + x2 + x3)/3
> or whatever.
> If I have to keep putting things into the data frame just to run my regressions, that's a pain in the ass. It's _not_ just "a bit more typing,"

I agree, it's not just a bit more typing.

> it's actual additional lines of code, it makes things harder to follow.

I think it makes it *easier* to follow if you're coming
at it fresh. It's no longer a narrative you have to tease
apart to figure out where you're at, but an encapsulated
block of code.

Sure, you can save some keystrokes, but it means more
ground to search for your reader.

- Bob

[Bob Carpenter]

Part II is just the technical issue.

But it's more confusing than I thought. The core concept here is
that R passes by value and that the value of (y ~ x) is a formula.
But it looks like it behaves differently than ("y ~ x") in a function
context.

So this all looks OK:

```
> foo <- function(n) { x <- rnorm(n); y <- 100 * x + 200 + rnorm(n); return(lm(y ~ x)); }
> foo(10)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x
199.6 100.4
```

But now what happens when I quote that?

```
> bar <- function(n) { x <- rnorm(n); y <- 100 * x + 200 + rnorm(n); return(lm("y ~ x")); }
> bar(10)
Error in eval(expr, envir, enclos) : object 'y' not found
```

Hmm. I honestly don't know what's going on there. But a little
spelunking shows me it's looking further up the environment chain:

```
> bar <- function(n) { x <- rnorm(n); y <- 100 * x + 200 + rnorm(n); return(lm("y ~ x")); }
> x <- rnorm(10); y <- rnorm(10)
> bar(10)

Call:
lm(formula = "y ~ x")

Coefficients:
(Intercept) x
0.02121 -0.20474
```

If you're in the top-level environment, now that x and y are defined,
quotes and no quotes appear to work the same way:

```
> lm(y ~ x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x
0.02121 -0.20474

> lm("y ~ x")

Call:
lm(formula = "y ~ x")

Coefficients:
(Intercept) x
0.02121 -0.20474
```


I'm confused. And rather frightened by this behavior.

- Bob

> On Feb 1, 2017, at 5:23 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>

[Bob Carpenter]

Ack. My mailer's treating the R input > as an embedded message.
And this isn't Discourse, so those ticks didn't do anything.
Hopefully it's still readable.

- Bob

> On Feb 3, 2017, at 2:54 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> Part II is just the technical issue.
>
> But it's more confusing than I thought. The core concept here is
> that R passes by value and that the value of (y ~ x) is a formula.
> But it looks like it behaves differently than ("y ~ x") in a function
> context.
>
> So this all looks OK:
>
> ``

> $ foo <- function(n) { x <- rnorm(n); y <- 100 * x + 200 + rnorm(n); return(lm(y ~ x)); }
> $ foo(10)

[Jonah Gabry]

On Friday, February 3, 2017 at 2:48:25 AM UTC-5, Bob Carpenter wrote:

We're not talking about a data frame here, but
a list, right?  Don't data frames have to be square?

For rstan a list, for rstanarm preferably a data frame (although list can be made to work too, as data frames are just fancy lists). Data frames do not have to be square (if by square you mean same number of rows and columns). 

 

I have no idea what Stata does.

Somehow Stata is still in business despite only allowing one dataset in memory at a time!

 

> Here's an example.  I start with this model (for simplicity I'm assuming we have a stan_logit() alias):
>
> stan_logit(y ~ dist + arsenic)
> or
> stan_logit(y ~ dist + arsenic, data="wells")

Are you passing the name of a list or something here?

I think that's just a typo. "wells" wouldn't work here. It should just be the object wells, which is a data frame. 

 

> Ok, fine so far.  Now I decide I want to code as dist/100:
> dist100 <- dist/100
> stan_logit(y ~ dist100 + arsenic)

This is where I think it's super helpful for those *reading*
your work (which includes you in the future) to keep
everything together.  At this point, you have defined y somewere
before the last fit and define a new dist100 here and everything's
just sitting in the global environment ready to be modified.

I completely agree, but I guess you were addressing to Andrew not me!

 

Jonah

[Jonah Gabry]

On Friday, February 3, 2017 at 2:54:41 AM UTC-5, Bob Carpenter wrote:

I'm confused.  And rather frightened by this behavior.

- Bob

Indeed, it's very confusing! Imagine trying to write the internals for rstanarm or lme4 or any modeling package when users are allowed to specify unquoted formulas or quoted formulas and either with or without explicitly passing data. Of the time I spend on rstanarm, probably 75% is related to problems stemming from this behavior and it's really frustrating. We could have a package that is both more stable and has a lot more features if we put more restrictions than R typically does on how the user can specify models. 

Jonah

[Bob Carpenter]

> On Feb 3, 2017, at 10:41 AM, Jonah Gabry <jga...@gmail.com> wrote:
>
> On Friday, February 3, 2017 at 2:48:25 AM UTC-5, Bob Carpenter wrote:
>
> We're not talking about a data frame here, but
> a list, right? Don't data frames have to be square?
>
> For rstan a list, for rstanarm preferably a data frame (although list can be made to work too, as data frames are just fancy lists). Data frames do not have to be square (if by square you mean same number of rows and columns).

I learn something new and surprising about R every day.

>
> I have no idea what Stata does.
>
> Somehow Stata is still in business despite only allowing one dataset in memory at a time!

So it's Andrew who wants the Stata-like behavior after all.

Here are three ways I see of proceeding:

I: Define data objects in the environment

y1 <- list(a = a1, b = b1);
y2 <- list(a = a2, b = b2);

fit1 <- stan("foo.stan", data=y1);
fit2 <- stan("foo.stan", data=y2);


II: Define data objects in the arguments

fit1 <- stan("foo.stan", data = list(a = a1, b = b1));
fit2 <- stan("foo.stan", data = list(a = a2, b = b2));

If you insist on calling only stan("foo.stan"), then you have
to do something destructive between calling the fits.


III: Define data objects globally

# load one global data set
a <- a1;
b <- b1;
fit1 <- stan("foo.stan")

# load a second global data set
a <- a2;
b <- b2;
fit2 <- stan("foo.stan");


I would never in a gazillion years want to recommend (III).
Personally, I prefer I or II. Michael's going further and
saying that at the point you have the y1 and y2 in the
first example, you should save them using rstan_dump or whatever
it's called so that you (a) have a record of the data you used
that's reproducible rather than the result of fiddling in
a REPL environment (https://en.wikipedia.org/wiki/Read–eval–print_loop),
and (b) have a way of sharing data and Stan file in a way that's
portable across Stan interfaces.

> > Here's an example. I start with this model (for simplicity I'm assuming we have a stan_logit() alias):
> >
> > stan_logit(y ~ dist + arsenic)
> > or
> > stan_logit(y ~ dist + arsenic, data="wells")
>
> Are you passing the name of a list or something here?
>
> I think that's just a typo. "wells" wouldn't work here. It should just be the object wells, which is a data frame.

This is the kind of thing where I never know what R's going to
do because of its flexibility with quotes in many situations.

- Bob

[Andrew Gelman]

Bob:

You keep talking about "global" things. R has local environments. For example, if you call Stan within a function, Stan will take the data from the local environment within the function call, not the global environment.

Here's a simple example:

foo <- function(){
a <- 3
print(a)
}

a <- 5
foo()

If you do the above, it will print 3, not 5.

P.S. stan() takes data in the form of a list. stan_glm() takes data in the form of a "data frame" (a matrix).

[Bob Carpenter]

> On Feb 3, 2017, at 4:14 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> Bob:
>
> You keep talking about "global" things. R has local environments.

Yes it does. And there's a top-level one called the global environment
about which R's doc says:

The global environment .GlobalEnv, more often known as the user's workspace,
is the first item on the search path. It can also be accessed by globalenv().
On the search path, each item's enclosure is the next item.

If you give R a variable (*not* a variable name), then
it searches up the environment chain from the current local
environment up to the global environment until it finds a definition
or fails.

Previously, RStan only read out of the global environment.
I didn't realize it got fixed so it'd read locally. I just
verified that you're right that it looks in the global environment
with RStan 2.14.

But lm() does *not* read out of the local environment when the
formula's given as a string.

$ x <- rnorm(100)

$ y <- 10 * x + 100 + rnorm(100)

$ foo <- function(n) { x <- rnorm(n); y <- 50 * x - 200 + rnorm(100); lm("y ~ x"); }

$ foo(200)

Call:
lm(formula = "y ~ x")

Coefficients:
(Intercept) x

99.89 10.03

So you can see it's reading out of the global environment, not
the function's local environment. If I define the formula without
quotes, it reads out of the local environment:

> bar <- function(n) { x <- rnorm(n); y <- 50 * x - 200 + rnorm(100); lm(y ~ x); }
> bar(200)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x

-200.0 49.9


No clue as to what's going on but this is why I'm talking about
local vs. global behavior and why Jonah says it's been such a pain
in the ass to mirror the behavior of lm() etc.

A language like C or C++ would let you read a global variable in a function,
but not if it's declared after the function. Everything has to
exist when a function's defined---that's static lexical scoping.
R uses dynamic lexical scoping, which is why it's so damn confusing.
It waits until run time, then sees if it can find a definition floating
around in the environment somewhere. It's also one of the reasons
it's so slow, but that's another issue.

The main reason static lexical scope is so desirable is that it
allows you to reason about what your program's going to do from
just looking at the program. With dynamic lexical scoping, you have
to wait until run time. At that point, someone could've accidentally
dropped in another definition of a variable between your global and
that function call, and that one would get picked up. This is why
we all think it's so dangerous as a general practice.

> For example, if you call Stan within a function, Stan will take the data from the local environment within the function call, not the global environment.
>
> Here's a simple example:
>
> foo <- function(){
> a <- 3
> print(a)
> }
>
> a <- 5
> foo()
>
> If you do the above, it will print 3, not 5.

Yes, I understand that. See above. You're passing
the a ***by value*** here, not by name. This is the
distinction we're trying to get through to you.

> P.S. stan() takes data in the form of a list. stan_glm() takes data in the form of a "data frame" (a matrix).

And apparently more general than a matrix. It has
types for columns and apparently doesn't have to be
rectangular, which was news to me.

- Bob

[Andrew Gelman]

Yes, R has local and global environments. It was just that I had the impression that you thought I wanted everything to be global.

[Andrew Gelman]

Mike:

This pipe thing looks cool but I don't know if I can easily integrate it into my workflow.  I maybe should be writing more functions.  But that creates problems too.  I've been having the naming problems that Hadley's talking about, so I'd like some solution, as my code is getting ugly.  My workspace gets filled up with objects like fit_1, stan_fit_1, sim_1, fit_2, stan_fit_2, etc.  It's just a mess.  Maybe I'll talk with Hadley and see if he can solve my problems.

A

[Jonah Sol Gabry]

It definitely takes some time to get used to the pipe, but it does make some (definitely not all) things in R much simpler, both practically and conceptually. It wouldn't do away with the problem at issue here though because it doesn't actually change how R works, it just provides a different interface. 

Jonah

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/gwaEpERn1Io/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I think the point is that with pipes you don't need
to write into intermediate results:

foo() %>% bar()

vs.

a = foo()
bar(a)

where you have to create a variable a to hold the
intermediate result of foo().

- Bob

[Bob Carpenter]

> On Feb 3, 2017, at 7:36 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>
> Yes, R has local and global environments. It was just that I had the impression that you thought I wanted everything to be global.

I want it to be passed by value from a local environment
as a function argument.

y = ...
stan(..., data=y)

You want it to be picked up non-locally from within
the function reaching back to the nested calling environment.

y = ...
stan(...)

Now you may think this isn't very different, but in the first
version, the function stan() receives the data y as an argument.
That is, it gets the value.

In the second version, what happens is that the function stan()
has to go back into the environment that called it and look
for the variable y. This is *completely* different.

The thing that can cause subtle bugs in a language that
dynamically scopes like R is that the path for searching
for a value *will* go all the way out to the global
environment. So if you're in some function body and forget
to define y, then Stan will keep looking further out.
But you don't know what that sequence of nested
environments are when you write the function---that's what
it means to be dynamic. So you might get things by mistake.

So the long story is that yes, you want to look in
the environment and that implicitly looks all the way up
to the global environment. I only want to look at an
argument. Only local scoping and no nasty run-time surprises.

- Bob

[Jonah Sol Gabry]

That's true, but for rstanarm it wouldn't really avoid the need to collect things into a data frame somewhere between creating the variables and calling the modeling function. For example, it would be cumbersome if not impossible to pipe y, x1, x2, etc., all into a modeling function unless they're combined into an object first. Or am I mistaken?

Jonah 

[Andrew Gelman]

Yes, the problem is that I do use the intermediate quantities, thus I'm stuck in naming hell, creating objects such as fit_1, sim_1, fit_2, sim_2, etc.

A

[Bob Carpenter]

I'd love to hear any suggestions anyone has. I have
this problem when fitting models, too. I try to be
careful to give all of my models meaningful names then
stick to that, but they wind up being too similar to
one another or too long.

One way to tame some of this complexity is to create
intermediate objects (even lists), as in:

run1 <- list(fit = fit1, sim = sim1)

The big win is that you then pass run1 around instead
of (fit1, sim1) and there are fewer names to memorize.
You may not like that you have to use run1$fit, but you
could always run it in an environment defined by run1.
Jonah cleaned up a lot of my case study code that way.

I actually find that clearer if there's a declared
structure somewhere than having all the bits and pieces
flying around. The use is much more consistent in that
it enforces the right naming. You can't make the mistake
of having sim1, fit1, sims2, fit2 (see the problem?).

- Bob

[Andrew Gelman]

Yes. maybe I should try this. Currently I'm saying:

fit_1 <- stan_glm(...)
sim_1 <- as.matrix(fit_1)
beta <- fit1_1[1:k,]

But I could just do
fit_1 <- stan_glm(...)
beta <- as_matrix(sim_1)[1:k,]

It's not clear what's better. The first formulation introduces this extra, unnecessary variable, but on the plus side it makes the steps more transparent. The second formulation has a lot happening in one line, which I think can be confusing to students, as they tend to read the code as a mysterious sequence of characters that they have to type in exactly right.

[Ben Goodrich]

On Saturday, February 4, 2017 at 9:37:46 PM UTC-5, Andrew Gelman wrote:

Yes. maybe I should try this.  Currently I'm saying:

fit_1 <- stan_glm(...)
sim_1 <- as.matrix(fit_1)
beta <- fit1_1[1:k,]

But I could just do
fit_1 <- stan_glm(...)
beta <- as_matrix(sim_1)[1:k,]

What is it that you are doing with beta that doesn't already have a convenience function to do it? We expose the as.matrix function because we can't think of everything that a user might want to do with the posterior draws but we can think of everything that an author wants to put in a book.

[Andrew Gelman]

I'm doing some predictions "by hand" (y = X*beta), also plotting some fitted regression lines on top of data.  I will send you book draft and I hope you can suggest some ways for me to make the code clearer!

A

[Ben Goodrich]

On Saturday, February 4, 2017 at 10:25:21 PM UTC-5, Andrew Gelman wrote:

I'm doing some predictions "by hand" (y = X*beta), also plotting some fitted regression lines on top of data.  I will send you book draft and I hope you can suggest some ways for me to make the code clearer!

Presumably the book can do that just once and then utilize posterior_predict() thereafter, right?

[Andrew Gelman]

Probably.  I'll just have to go thru the examples one at a time and see!

[Bob Carpenter]

> fit_1 <- stan_glm(...)
> beta <- as_matrix(sim_1)[1:k,]

If readers have problem with this, then they've got bigger
problems than learning Stan.

I was suggesting grouping variables together into objects
so they can be passed together and have internally consistent
naming access. It's easier to think about chunks than all
these individual arguments, especially if the work's being
repeated all over the place.

- Bob