Random slopes model

[Joe King]

Hello !

I have been using the code below to fit a random intercepts (variance components) model.

data {
    int<lower=1> n2 ; 
    int<lower=1> N ;
    real y[N] ;
    int school[N] ;
}
parameters {
    real u2_z[n2] ;
    real<lower=.1,upper=100> sigma ;
    real<lower=.1,upper=100> sigma_u2;
    real <lower=0,upper=1000> beta;
}
model {
    real mu[N];
    for (i in 1:N)
      mu[i] <- beta + u2_z[school[i]] * sigma_u2;

    y ~ normal(mu,sigma);
    u2_z ~ normal(0, 1) ;
}

This would be fitted in lme4 as

lmer(outcome ~1 + (1|school.ID),data=dt.all)

I now want to include a term for random slope, which would be fitted in lme4 as

lmer(outcome ~1 + (year|school.ID),data=dt.all)
or
lmer(outcome ~1 + (year-1|school.ID),data=dt.all)

The variable year is a factor with 3 levels., Please can anyone help me modify the stan model to implement these ?

Thanks  !
JK

[Bob Carpenter]

You just need to add another term to mu[i]
for the slope * predictor.

Assuming the predictor is x[i], the loop's just:

for (i in 1:N) {
real intercept;
real slope;
intercept <- beta + u2_z[school[i]] * sigma_u2;
slope <- gamma + v2_z[school[i]] * sigma_v2;
mu[i] <- intercept + slope * x[i];
}

And you can use the same implicit and explicit priors
for gamma, v2_z and sigma_v2 as you used for beta, u2_z
and sigma_u2 (assuming you also want to restrict gamma to
be positive).

You've already applied what we've taken to calling the "Matt trick"
(after Matt Hoffman, who also invented NUTS) to your code, which perhaps
makes it a little harder to see what's going on (but much more efficient
to sample from). If you instead wrote your model as:

intercept ~ normal(beta,sigma_u2);
y ~ normal(intercept[school[i]],sigma);

it should be clearer that the slope + intercept model is:

intercept ~ normal(beta,sigma_u2);
slope ~ normal(gamma,sigma_v2);
y ~ normal(intercept[school[i]] + slope[school[i]] * x[i], sigma);

Ideally, Stan would figure out the transforms on its own,
but not yet.

- Bob

[Joe King]

Hi Bob

Thanks for your reply, though I am a bit confused. There are no predictors in this model (yet). I am trying to replicate the output I get from lme4 thus:

> fit.lme4.all <- lmer(outcome ~1 + (year|school.ID),data=dt.all)
> summary(fit.lme4.all)
Linear mixed model fit by REML
Formula: outcome ~ 1 + (year | school.ID)
   Data: dt.all
    AIC    BIC logLik deviance REMLdev
 183343 183410 -91664   183325  183327
Random effects:
 Groups    Name        Variance  Std.Dev. Corr         
 school.ID (Intercept)   0.94242  0.97079              
           year2005     20.83882  4.56496 -0.197       
           year2006    334.04831 18.27699 -0.283 -0.073
 Residual               25.26760  5.02669              
Number of obs: 30000, groups: school.ID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  19.9286     0.1021   195.2

Thanks again
JK

[Daniel Lee]

On Oct 21, 2012, at 6:31 AM, Joe King <joeki...@yahoo.com> wrote:

Hi Bob

Thanks for your reply, though I am a bit confused. There are no predictors in this model (yet).

I think Bob meant to type "year[i]" instead of "school[i]" or "x[i]" in his last post. So the loop would look something like:

for (i in 1:N) {
   real intercept;
   real slope;
   intercept <- beta + u2_z[school[i]] * sigma_u2;
   slope <- gamma + v2_z[year[i]] * sigma_v2;
   mu[i] <- intercept + slope;
 }

--
 
 

[Joe King]

Thanks Daniel

But now my model does not parse:


TRANSLATING MODEL 'schools_code' FROM Stan CODE TO C++ CODE NOW.
Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
  failed to parse Stan model 'schools_code' with error message:
LOCATION:  file=input; line=27, column=14

   intercept <- beta + u2_z[school[i]] * sigma_u2;

             ^-- here

DIAGNOSTIC(S) FROM PARSER:

Parser expecting: ";"

Here's the code:


schools_code <- '

data {
    int<lower=1> n2 ; 
    int<lower=1> N ;

    int<lower=1> N_year;

    real y[N] ;
    int school[N] ;

    int year[N] ;
}
parameters {
    real u2_z[n2] ;
    real v2_z[N_year] ;

    real<lower=.1,upper=100> sigma ;
    real<lower=.1,upper=100> sigma_u2;
    real <lower=0,upper=1000> beta;

    real<lower=.1,upper=100> gamma ;
    real<lower=.1,upper=100> sigma_v2;

}
model {

      real mu[N];

    for (i in 1:N) {
       real intercept;
       real slope;
       intercept <- beta + u2_z[school[i]] * sigma_u2;
       slope <- gamma + v2_z[year[i]] * sigma_v2;
       mu[i] <- intercept + slope;
     }

    y ~ normal(mu,sigma);
    2_z ~ normal(0, 1) ;
    2_z ~ normal(0, 1) ;
}

'
schools_dat <- list(N = 30000, n2 = 100, N_year = 3,
                    y = dt.all$outcome,
              year=dt.all$year,
                    school = dt.all$school.ID)

fit <- stan(model_code = schools_code, data = schools_dat,
            iter = 800, chains = 1)

Thanks in advance
JK

[Joe King]

Whoops, the last 2 lines of  model {} should be

    u2_z ~ normal(0, 1) ;
    v2_z ~ normal(0, 1) ;

On Sunday, October 21, 2012 8:49:56 AM UTC+1, Joe King wrote:

[Joe King]

Sorry, I remember now the guidance about posting attachments with code that doesn't work

So please find the code and data attached

---

From: Joe King <joeki...@yahoo.com>
To: stan-...@googlegroups.com
Sent: Sunday, October 21, 2012 12:57 PM
Subject: [stan-users] Re: Random slopes model

--
 
 

[Ben Goodrich]

This is weird. It looks like a parser bug to me. @Bob, the attached doesn't parse. But if you change the loop in the model block to

for (i in 1:N) {

  mu[i] <- beta + u2_z[school[i]] * sigma_u2 + gamma + v2_z[year[i]] * sigma_v2 ;

}

it parses and compiles. You will then run into a problem where your year variable is a factor rather than a numeric variable. Then it seems to run slowly, but maybe that is to be expected with an N of 30000.

Ben

[Joe King]

Hello Ben

After making the change you suggested, I got the following compilation error:

Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) :
  variable year not found.

Is this what you meant by "You will then run into a problem where your year variable is a factor rather than a numeric variable" ? I didn't know how to specify a factor in the data{} section : can you tell me how to do that ?

In the meantime, I changed the year variable to numeric, and the prior on v2_z:
    v2_z ~ gamma(0.001000,0.001000) ;

Then it compiled and seemed to be running, until I got the popup:

Microsoft Visual C++ Runtime Library

This application has requested the Rutime to terminate in an unusual way.
Please contact the application's support team for more information.

And the message that R had stopped working, to which I answered "close", though it didn't close down and the process seemed to be still in a loop (thrashing the CPU but unresponsive) so I had to kill the process.

Please advise

Thanks
JK

[Daniel Lee]

Yes, this is a weird parser bug. I'll send a short model to stan-dev soon.

> --
>
>

[Daniel Lee]

Actually, the last model Ben sent will do the trick. If you change these lines

    intercept <- beta + u2_z[school[i]] * sigma_u2 ;

    slope <- gamma + v2_z[year[i]] * sigma_v2 ;

with

    slope <- gamma + v2_z[year[i]] * sigma_v2 ;
    intercept <- beta + u2_z[school[i]] * sigma_u2 ;

it compiles.

I haven't found the minimal model to show this behavior, but once we narrow it down, we should add it to our parser test.

[Ben Goodrich]

On Sunday, October 21, 2012 9:51:00 AM UTC-4, Joe King wrote:

After making the change you suggested, I got the following compilation error:

Error in new_CppObject_xp(fields$.module, fields$.pointer, ...) :
  variable year not found.

Is this what you meant by "You will then run into a problem where your year variable is a factor rather than a numeric variable" ? I didn't know how to specify a factor in the data{} section : can you tell me how to do that ?

That is what I meant. You can't specify a factor in the data{} section of a .stan file. You would need to convert the factor inside R to integer codes. So, if for example, you had a factor in R like

data(InsectSprays)
summary(InsectSprays$spray)

You would need to coerce it to an integer vector like

int_spray <- as.integer(InsectSprays$spray)
summary(int_spray)

  

Then it compiled and seemed to be running, until I got the popup:

Microsoft Visual C++ Runtime Library

This application has requested the Rutime to terminate in an unusual way.
Please contact the application's support team for more information.

And the message that R had stopped working, to which I answered "close", though it didn't close down and the process seemed to be still in a loop (thrashing the CPU but unresponsive) so I had to kill the process.

Please advise

According to the Stan documentation, the MSVC compiler can't compile Stan models, which your experience seems to be inconsistent with. In R, can you tell us the output of

system("R CMD config CXX")

Hopefully, someone who uses Windows will know more than I do about this.

Ben

 

[Joe King]

Thanks Ben

The result of system ("R CMD config CXX") is

g++ -m64

The problem is different when I use as.integer instead of as.numeric. In that case I get this error

SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).
Error in sampler$call_sampler(args_list[[i]]) : std::bad_alloc

Note: it runs OK with
v2_z ~ normal(0, 1) ;

but not
v2_z ~ gamma(0.001000,0.001000) ;

Thanks again
JK

---

From: Ben Goodrich <goodri...@gmail.com>
To: stan-...@googlegroups.com
Cc: Joe King <joeki...@yahoo.com>
Sent: Sunday, October 21, 2012 3:20 PM
Subject: Re: [stan-users] Re: Random slopes model

On Sunday, October 21, 2012 9:51:00 AM UTC-4, Joe King wrote:

After making the change you suggested, I got the following compilation error:

Error in new_CppObject_xp(fields$. module, fields$.pointer, ...) :
  variable year not found.

Is this what you meant by "You will then run into a problem where your year variable is a factor rather than a numeric variable" ? I didn't know how to specify a factor in the data{} section : can you tell me how to do that ?

[Bob Carpenter]

Parsing and terminology issues inline below.

On 10/21/12 8:59 AM, Ben Goodrich wrote:

...
> model {
> real intercept;
> real slope;
> real mu[N];
>
> for (i in 1:N) {
> intercept <- beta + u2_z[school[i]] * sigma_u2 ;
> slope <- gamma + v2_z[year[i]] * sigma_v2 ;

> mu[i] <- intercept + slope ;
> }

...

Ben's right. This should parse. And it does if
the variable "intercept" is renamed to something
that doesn't start with "int". This will be a problem
for any variable whose name has a prefix that matches
a ype (like "real").

The problem is that the "int" in "intercept" is
getting parsed as the beginning of a declaration
of variable "ercept"! I need to only match "int"
as starting a declaration if the next char is not
alphanumeric or "_".

I'll commit a patch for the 1.0.3 release next week.

On 10/21/12 7:02 AM, Daniel Lee wrote:> On Oct 21, 2012, at 6:31 AM, Joe King <joeki...@yahoo.com

<mailto:joeki...@yahoo.com>> wrote:
>
>> Thanks for your reply, though I am a bit confused. There are no predictors in this model (yet).
>
> I think Bob meant to type "year[i]" instead of "school[i]" or "x[i]" in his last post.

I may have the terminology wrong here, but I thought
a "slope" was something you multiplied by a predictor.
The model Daniel suggested (see above) looks to me
like a multilevel model with two varying intercepts.

- Bob

[Joe King]

Dear Bob

Please accept my apologies for wrong use of terminology. I only began learning statistics in the last few months and I find much of the terminology very confusing. I learned about multilevel (or is it "hierarchical" ?)  models from my maths teacher in June and the method I learned was to first fit a "null model" with no covariates, so I am sorry for calling this a "random slope" model - but it will be after I add the covariates in the next stage I think.

Thank you for helping me ! I appreciate it very much !
JK

---

From: Bob Carpenter <ca...@alias-i.com>
To: stan-...@googlegroups.com
Sent: Sunday, October 21, 2012 9:02 PM
Subject: Re: [stan-users] Re: Random slopes model

--

[Bob Carpenter]

On 10/21/12 4:38 PM, Joe King wrote:
> Dear Bob
>
> Please accept my apologies for wrong use of terminology. I only began learning statistics in the last few months and I
> find much of the terminology very confusing. I learned about multilevel (or is it "hierarchical" ?) models from my
> maths teacher in June and the method I learned was to first fit a "null model" with no covariates, so I am sorry for
> calling this a "random slope" model - but it will be after I add the covariates in the next stage I think.

They're called "hierarchical" when the effects are
all nested inside each other, such as students
within classrooms within school districts within
states.

They're called "multilevel" when the effects are not
nested. Years and schools are cross-cutting categories,
so a model with both is multilevel, not hierarchical.

As to "slope" and "intercept", a typical univariate linear
regression of y on x would be:

for (n in 1:N)
y[n] ~ normal(alpha * x[n] + beta, sigma);

where alpha is the "slope" and beta the "intercept" of the
line defining the location of the normal distribution generating
y[n] based on x[n].

> Thank you for helping me ! I appreciate it very much !

Sorry the program ran into a bug with our syntax right away!
Thanks for letting us know, because I was able to fix it
for the next release.

I love the idea of starting learning stats with
multilevel models in Stan! There's a section in the
user's guide, but for more info, you should really
read Gelman and Hill's 2009 regression book. We
more or less built Stan to deal with extending
the gnarly models in their book to even more complex
structures.

- Bob

[Ben Goodrich]

On Sunday, October 21, 2012 10:45:59 AM UTC-4, Joe King wrote:

Thanks Ben

The result of system ("R CMD config CXX") is

g++ -m64

The problem is different when I use as.integer instead of as.numeric. In that case I get this error

SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).
Error in sampler$call_sampler(args_list[[i]]) : std::bad_alloc

I am surprised that it makes a difference. Maybe Stan is expecting everything to be numeric even if things are declared as int?

 

Note: it runs OK with
v2_z ~ normal(0, 1) ;

but not
v2_z ~ gamma(0.001000,0.001000) ;

That makes sense because the point of the Matt trick is to use standard normal priors wherever possible, because they lead to the most efficient sampling. Putting a gamma(0,0) prior on v2_z is just wrong as far as I can tell.

I managed to get it to run 32 chains each for 2000 iterations (on Linux) without getting the runtime error you got before. It didn't converge, but this example seems to have the same problem as the example from Shige: The default initial values are very far from the truth. Standardizing y or providing reasonable initial values should be sufficient.

Ben

[Bob Carpenter]

On 10/21/12 7:17 PM, Ben Goodrich wrote:
> On Sunday, October 21, 2012 10:45:59 AM UTC-4, Joe King wrote:

> The result of system ("R CMD config CXX") is
>
> g++ -m64
>
> The problem is different when I use as.integer instead of as.numeric.

Stan just has a base reader class from which a model
reads data.

In the command-line code, it's implemented
with a dump reader (which doesn't yet handle all the
creative variation in R's dump() output).

In RStan, it's implemented in RStan to read from
an R object. So that'd be the place to look for
what's going on with the types.

I'm guessing the following intended behavior may have
something to do with it:

The model will issue calls to read integer arrays
(for integer arrays in a Stan model) or double arrays
(for everything else). If the input is integers, it
should still be readable as a double by casting.

That's what the dump() reader implementation does and
is what the RStan implementation should do.

- Bob

[Jiqiang Guo]

Could someone let me know if there is still a difference between as.integer and as.numeric used on the data? 

However, there should not be a problem as in RStan, there is data preprocess, which will turn double to int if no precision is lost.  For example, 

> f <- gl(3,4)

> f

 [1] 1 1 1 1 2 2 2 2 3 3 3 3

Levels: 1 2 3

> lst <- list(a=as.integer(f), b=as.numeric(f))

> lst

$a

 [1] 1 1 1 1 2 2 2 2 3 3 3 3

$b

 [1] 1 1 1 1 2 2 2 2 3 3 3 3

> is.integer(lst$a)

[1] TRUE

> is.integer(lst$b)

[1] FALSE

> lst2 <- rstan:::data_preprocess(lst)

> identical(lst2$a, lst2$b)

[1] TRUE

> is.integer(lst2$a)

[1] TRUE

> is.integer(lst2$b)

[1] TRUE

--

Jiqiang

- Bob

--

[Joe King]

Hi Ben

Why is it wrong to put a gamma prior on v2_z ? It's a variance so it can't be negative.....

Sorry to keep going on about it, but this is what I am trying to replicate (I highlighted them in red)

> fit.lme4.all <- lmer(outcome ~1 + (year-1|school.ID),data=dt.all)

> summary(fit.lme4.all)
Linear mixed model fit by REML

Formula: outcome ~ 1 + (year - 1 | school.ID)

   Data: dt.all
    AIC    BIC logLik deviance REMLdev
 183343 183410 -91664   183325  183327
Random effects:
 Groups    Name     Variance  Std.Dev. Corr         

 school.ID year2004   0.94242  0.97078              
           year2005  20.03281  4.47580  0.016       
           year2006 324.96381 18.02675 -0.233 -0.137

 Residual            25.26760  5.02669              
Number of obs: 30000, groups: school.ID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  19.9286     0.1021   195.2

Thanks again !

JK

---

From: Ben Goodrich <goodri...@gmail.com>
To: stan-...@googlegroups.com
Cc: Joe King <joeki...@yahoo.com>

Sent: Monday, October 22, 2012 12:17 AM

Subject: Re: [stan-users] Re: Random slopes model

On Sunday, October 21, 2012 10:45:59 AM UTC-4, Joe King wrote:

Thanks Ben

The result of system ("R CMD config CXX") is

g++ -m64

The problem is different when I use as.integer instead of as.numeric. In that case I get this error

SAMPLING FOR MODEL 'schools_code' NOW (CHAIN 1).

Error in sampler$call_sampler(args_ list[[i]]) : std::bad_alloc

I am surprised that it makes a difference. Maybe Stan is expecting everything to be numeric even if things are declared as int?

 

Note: it runs OK with
v2_z ~ normal(0, 1) ;

but not
v2_z ~ gamma(0.001000,0.001000) ;

That makes sense because the point of the Matt trick is to use standard normal priors wherever possible, because they lead to the most efficient sampling. Putting a gamma(0,0) prior on v2_z is just wrong as far as I can tell.

I managed to get it to run 32 chains each for 2000 iterations (on Linux) without getting the runtime error you got before. It didn't converge, but this example seems to have the same problem as the example from Shige: The default initial values are very far from the truth. Standardizing y or providing reasonable initial values should be sufficient.

Ben

--
 
 

[Bob Carpenter]

Joe said:

> v2_z ~ gamma(0.001000,0.001000) ;

Then Ben said:

> .. Putting a gamma(0,0) prior on v2_z is just wrong as far as I can tell.


And Joe then asked:

> Why is it wrong to put a gamma prior on v2_z ? It's a variance so it can't be negative.....

It's not "wrong" in the sense that it's illegal,
but wrong in the sense that it's too concentrated
around zero.

Check out the attached plot of log gamma(x|0.001,0.001)
vs. x.

Ben also said:

> I managed to get it to run 32 chains each for 2000 iterations (on Linux) without getting the runtime error you got
> before. It didn't converge, but this example seems to have the same problem as the example from Shige: The default
> initial values are very far from the truth. Standardizing y or providing reasonable initial
> values should be sufficient.

This is a separate issue. You see the same
kind of thing in some of R's regression if you
give it bad starting points.

- Bob

[Ben Goodrich]

On Monday, October 22, 2012 2:51:35 AM UTC-4, Joe King wrote:

Hi Ben

Why is it wrong to put a gamma prior on v2_z ? It's a variance so it can't be negative.....

Are we still still talking about the same model here? Originally you had

parameters {
    real u2_z[n2] ;
    real v2_z[N_year] ;
...

If you want v2_z to be positive, you could impose that constraint here. But it wouldn't make sense to do so otherwise.

Ben

 

[Joe King]

Hello Bob

Thank you for your reply and your explanation about multilevel/hierarchical. Do you mean this book:

http://www.amazon.co.uk/Analysis-Regression-Multilevel-Hierarchical-Models/dp/052168689X/sr=8-2/qid=1157677172/ref=pd_bbs_2/102-8078425-4555328?ie=UTF8&s=books

?

I asked my brother and he said he can borrow it for me from his Uni library, though he won't come home again for a few weeks so I'll have to wait a bit.

Thanks again

Joe

---

From: Bob Carpenter <ca...@alias-i.com>
To: stan-...@googlegroups.com

Sent: Sunday, October 21, 2012 10:32 PM

Subject: Re: [stan-users] Re: Random slopes model

--

[Joe King]

Hi again Bob

What is the recommended way to make a "uninformative" prior for a variance ? My maths teacher told me about gamma (from what you said it seems she is wrong, so I will print the attached and ask her about it today), but my brother said I should use a Jeffrey's prior. I also read about a half-cauchy distribution !

OK off to school now. Will check back here this afternoon.

Thanks again

JK

---

From: Bob Carpenter <ca...@alias-i.com>
To: stan-...@googlegroups.com

Sent: Monday, October 22, 2012 8:05 AM

Subject: Re: [stan-users] Re: Random slopes model

--

[Ben Goodrich]

On Monday, October 22, 2012 3:05:16 AM UTC-4, Bob Carpenter wrote:

Joe said:

 >     v2_z ~ gamma(0.001000,0.001000) ;

Then Ben said:

> .. Putting a gamma(0,0) prior on v2_z is just wrong as far as I can tell.


And Joe then asked:

 > Why is it wrong to put a gamma prior on v2_z ? It's a variance so it can't be negative.....

It's not "wrong" in the sense that it's illegal,
but wrong in the sense that it's too concentrated
around zero.

Something is wrong here. In the original .stan file, real v2_z[N_year] isn't constrained to be positive. So either the declaration is wrong or the gamma prior isn't appropriate for v2_z. If that prior makes sense to you Joe, then Stan will sample a lot better if you declare it as real<lower=0> v2_z[N_year].

Ben

[Bob Carpenter]

On 10/22/12 3:11 AM, Joe King wrote:

> ... Do you mean this book:
>
> http://www.amazon.co.uk/Analysis-Regression-Multilevel-Hierarchical-Models/dp/052168689X/sr=8-2/qid=1157677172/ref=pd_bbs_2/102-8078425-4555328?ie=UTF8&s=books

Yes.

- Bob

[Bob Carpenter]

On 10/22/12 3:20 AM, Joe King wrote:

> What is the recommended way to make a "uninformative" prior for a variance ?

It's a red herring as far as modeling goes. Any choice
of prior is making a modeling decision. We tend to prefer
weakly informative priors.

> My maths teacher told me about gamma (from
> what you said it seems she is wrong, so I will print the attached and ask her about it today), but my brother said I

> should use a Jeffrey's prior.I also read about a half-cauchy distribution !

Half Cauchy is the easiest in Stan, because you don't even
need the half part because we only need to compute the probability
up to a proportion. For instance, you can use something like:

real<lower=0> sigma;

...

sigma ~ cauchy(0,5);

- Bob

[Joe King]

Thanks again Bob

Can you tell me how I provide initial values in RStan ? I looked in the manual but didn't see anything except in 4.2 about command line options...

JK

---

From: Bob Carpenter <ca...@alias-i.com>
To: stan-...@googlegroups.com

Sent: Monday, October 22, 2012 6:04 PM

Subject: Re: [stan-users] Re: Random slopes model

--

[Bob Carpenter]

On 10/22/12 3:57 PM, Joe King wrote:> Thanks again Bob

>
> Can you tell me how I provide initial values in RStan ? I looked in the manual but didn't see anything except in 4.2
> about command line options...

The manual doesn't cover RStan arguments. Use
help(rstan) and follow the index for a list of
functions and help(stan) for the arguments to stan().

There are multiple ways to do it, but the "standard"
way to do it is with a function.

Here's an example from a model I just ran:

Model params (Stan)

...
parameters {
simplex[K] pi; // prevalence of categories
simplex[K] theta[J,K]; // response of anotator j to category k
}
...

R inits (R)

piVote <- rep(0,K);
...
thetaVote <- array(NA,c(J,K,K));
...

Define the function (R):

init_fun = function() {
return(list(pi=piVote,
theta=thetaVote));
}

Call the fit function (R):

fit <- stan(file="dawid-skene.stan",
init=init_fun,
data=data_list,
iter=500, chains=4);

- Bob

[Ben Goodrich]

On Monday, October 22, 2012 2:51:35 AM UTC-4, Joe King wrote:

Sorry to keep going on about it, but this is what I am trying to replicate (I highlighted them in red)

> fit.lme4.all <- lmer(outcome ~1 + (year-1|school.ID),data=dt.all)
> summary(fit.lme4.all)
Linear mixed model fit by REML
Formula: outcome ~ 1 + (year - 1 | school.ID)
   Data: dt.all
    AIC    BIC logLik deviance REMLdev
 183343 183410 -91664   183325  183327
Random effects:
 Groups    Name     Variance  Std.Dev. Corr         
 school.ID year2004   0.94242  0.97078              
           year2005  20.03281  4.47580  0.016       
           year2006 324.96381 18.02675 -0.233 -0.137
 Residual            25.26760  5.02669              
Number of obs: 30000, groups: school.ID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  19.9286     0.1021   195.2

 
I have finally been able to replicate this result reasonably well. The fundamental problem was that the above lmer() model is fairly different from the original .stan code. This model has a single, fixed intercept and the coefficients on the 3 year variables vary by school (of which there are 100). I tend to think of a 100x3 matrix of coefficients whose rows are assumed to be multivariate normal with an unknown mean vector, an unknown diagonal matrix of standard deviations, and an unknown correlation matrix. So, in order to get something that resembles the above, we have to parameterize the .stan file accordingly. I'll attach the entire code but discuss it in blocks because it is a useful example.

The data are the same, except that I centered the outcome variable in R.

data {
    int<lower=1> n2 ;    // number of schools
    int<lower=1> N ;     // number of students
    int<lower=1> N_year; // number of years
    real y[N] ;          // outcome
    int school[N] ;      // school id
    int year[N] ;        // year {1,2 or 3}
}

The parameters are a bit different than what you originally specified. The year variable has three levels, and I'm treating the first as the baseline. Thus, I'm constraining the expectation of the first coefficient to be zero and letting the other two contrast with it.

parameters {
    real alpha;            // intercept
    vector[N_year-1] beta; // 2nd and 3rd expected slopes (first is constrained to be 0.0)
    real<lower=0> sigma ;  // standard error of regression
    matrix[n2,N_year] slopes_z; // z-score of the slopes
    corr_matrix[N_year] Tau; // correlation matrix of the slopes
    real<lower=0> tau[N_year]; // standard deviations of the slopes
}

I now use the "Multivariate Matt Trick" (also mentioned by Terrance) where quantities that have a multivariate normal distribution under the theoretical model are in practice, expressed as a function of standard normals and a Cholesky factor. See

http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Drawing_values_from_the_distribution

model {
  real mu[N];
  int k;
  matrix[n2,N_year] slopes_error;
  slopes_error <- slopes_z * transpose(cholesky_decompose(Tau)) ;

  for (i in 1:N) {

    k <- year[i];
    mu[i] <- alpha + (if_else(k - 1, beta[k - 1], 0.0) +
                       tau[k] * slopes_error[school[i],k]);
  }

What I call slopes_error (almost) corresponds to transpose(Az) in the Wikipedia notation. The difference is that that Tau is a correlation matrix rather than a covariance matrix, so we have to scale by the standard deviations (tau). The last line says that if k > 1, add beta[k-1] plus tau[k] times the appropriate element of schools_error, and if k = 1, then just use 0.0 for the expected baseline. The rest of the model {} block is straightforward, just a likelihood and a bunch of weakly informative priors priors (except on Tau, which is implicitly jointly uniform):

  y ~ normal(mu,sigma);
  for (j in 1:N_year) {
     tau[j] ~ cauchy(0.0,5.0) ;
     col(slopes_z,j) ~ normal(0, 1) ;
  }
  alpha ~ cauchy(0.0,5.0);
  beta ~ cauchy(0.0,5.0);
  sigma ~ cauchy(0.0,5.0);
  // implicit: Tau ~ lkj_corr(1.0);
}

Even with all that, it is still a bit difficult to sample from this posterior efficiently. I get

Inference for Stan model: schools_code.
8 chains: each with iter=2000; warmup=1000; thin=1; 2000 iterations saved.

           mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff  Rhat
alpha    -0.786   0.002 0.112 -1.001 -0.862 -0.786 -0.712 -0.561  2106 1.003
beta[1]   0.405   0.025 0.466 -0.518  0.090  0.416  0.722  1.311   336 1.037
beta[2]   1.629   0.074 1.783 -1.895  0.412  1.592  2.816  5.213   587 1.018
sigma     5.027   0.000 0.021  4.986  5.012  5.027  5.041  5.068  8000 1.000
tau[1]    0.993   0.002 0.090  0.836  0.929  0.987  1.051  1.185  3077 1.002
tau[2]    4.579   0.012 0.338  3.979  4.341  4.554  4.791  5.289   740 1.012
tau[3]   18.080   0.057 1.237 15.751 17.248 18.049 18.902 20.589   466 1.019
Tau[1,1]  1.000   0.000 0.000  1.000  1.000  1.000  1.000  1.000  8000   NaN
Tau[1,2]  0.013   0.006 0.110 -0.200 -0.067  0.014  0.093  0.218   315 1.038
Tau[1,3] -0.217   0.005 0.105 -0.415 -0.290 -0.220 -0.147 -0.008   388 1.015
Tau[2,1]  0.013   0.006 0.110 -0.200 -0.067  0.014  0.093  0.218   315 1.038
Tau[2,2]  1.000   0.000 0.000  1.000  1.000  1.000  1.000  1.000  7254 0.999
Tau[2,3] -0.136   0.004 0.093 -0.318 -0.201 -0.135 -0.072  0.040   633 1.011
Tau[3,1] -0.217   0.005 0.105 -0.415 -0.290 -0.220 -0.147 -0.008   388 1.015
Tau[3,2] -0.136   0.004 0.093 -0.318 -0.201 -0.135 -0.072  0.040   633 1.011
Tau[3,3]  1.000   0.000 0.000  1.000  1.000  1.000  1.000  1.000  6610 0.999

I marked the relevant comparisons to the lmer() results in red. Although the result for alpha here looks very different from the 19.93 produced by lmer(), keep in mind that I centered the outcome, which unfortunately often seems to be necessary at this point in Stan's development.

Ben

[Joe King]

Hi Ben

This is great - thank you so much !

Do you mean I should centre the outcome in R first with something like
outcome <- outcome - mean(outcome)

?

Also, would it make much difference if I use a implicit uniform prior on the parameters instead of cauchy ? I know how to do that with <lower=, upper> in the parameter section but how do I do that for the vector ?

Thanks !!
JK

---

From: Ben Goodrich <goodri...@gmail.com>
To: stan-...@googlegroups.com
Cc: Joe King <joeki...@yahoo.com>

Sent: Tuesday, October 23, 2012 10:38 AM

Subject: Re: [stan-users] Re: Random slopes model

[Ben Goodrich]

On Tuesday, October 23, 2012 3:12:52 PM UTC-4, Joe King wrote:

Do you mean I should centre the outcome in R first with something like
outcome <- outcome - mean(outcome)

Yes, for a fully observed outcome (i.e. the sample mean can be estimated), that is often a good idea with Stan and is what I did. In principle, it shouldn't be necessary but in practice it makes the default initial values of the parameters more plausible. There is not much cost to it, as long as you remember that it shifts the estimated intercept of the model.

Also, would it make much difference if I use a implicit uniform prior on the parameters instead of cauchy ? I know how to do that with <lower=, upper> in the parameter section but how do I do that for the vector ?

You can do it the same way, you just indicate the length of the vector at the end. In fact, I did so in the previous post

real<lower=0> tau[N_year]; // standard deviations of the slopes

An implicit uniform prior on the taus over some reasonable bounded interval would probably yield about the same results. An improper or barely proper prior probably wouldn't yield the same results.

Ben

[Joe King]

Hi Ben

I tried (in the model{})

//  alpha ~ cauchy(0.0,5.0);  //now implicitly uniform

//  sigma ~ cauchy(0.0,5.0); //now implicitly uniform

with (in the parameters {} )

real<lower=-10,upper=10> alpha;            // intercept
real<lower=0,upper=20> sigma ;  // standard error of regression

That seemed to be working, but then I put

//  beta ~ cauchy(0.0,5.0);  //now implicitly uniform

and

vector<lower=0,upper=20>[N_year-1] beta; // expected slopes (first is constrained to be 0.0)

which caused it to fail to parse. That's what I meant when I was talking about constraining the vector.....

Please advise

Thanks

JK

---

From: Ben Goodrich <goodri...@gmail.com>
To: stan-...@googlegroups.com
Cc: Joe King <joeki...@yahoo.com>

Sent: Tuesday, October 23, 2012 8:26 PM

Subject: Re: [stan-users] Re: Random slopes model

[Andrew Gelman]

OK, now I have a question (forgive me if the answer's already in the manual).  Suppose I define a parameter theta to have a normal (mu, sigma) distribution but then theta is declared as real<lower=0, upper=0.5>.  Does Stan then work with the truncated normal density (i.e., getting the correct normalizing factor which is now a function of mu and sigma) or does it just use the normal density without adjusting for the truncation?

A

On Oct 23, 2012, at 3:26 PM, Ben Goodrich wrote:

An implicit uniform prior on the taus over some reasonable bounded interval would probably yield about the same results. An improper or barely proper prior probably wouldn't yield the same results.

Ben

--
 
 

[Bob Carpenter]

On 10/23/12 4:32 PM, Andrew Gelman wrote:
> OK, now I have a question (forgive me if the answer's already in the manual). Suppose I define a parameter theta to
> have a normal (mu, sigma) distribution but then theta is declared as real<lower=0, upper=0.5>. Does Stan then work with
> the truncated normal density (i.e., getting the correct normalizing factor which is now a function of mu and sigma) or
> does it just use the normal density without adjusting for the truncation?

First, you never get the full normalization term in Stan.
All constant terms are dropped. In the normal, we never
include the constant normalizing term -log(sqrt(2*pi)).

The difference between

x ~ normal(mu,sigma); // no space before '('!!!

and

x ~ normal(mu,sigma) T[0,0.5];

is that the latter applies the truncated density and
the former does not.

If the variable x is declared as real<lower=0,upper=0.5>,
that is, with fixed bounds, then the first is OK as
we only need the density up to a proportion. If the
bounds were variable, you'd need the second form.

Assuming the bounds are properly declared for x,
we should be able to implicitly drop the truncation
in the future, but the code doesn't do that yet (it also
doesn't check ranges for x and return -inf if x is not
in the specified bounds and it probably should).

- Bob

[Ben Goodrich]

On Tuesday, October 23, 2012 4:18:18 PM UTC-4, Joe King wrote:

Hi Ben

I tried (in the model{})

//  alpha ~ cauchy(0.0,5.0);  //now implicitly uniform

//  sigma ~ cauchy(0.0,5.0); //now implicitly uniform

with (in the parameters {} )

real<lower=-10,upper=10> alpha;            // intercept
real<lower=0,upper=20> sigma ;  // standard error of regression

That seemed to be working, but then I put

//  beta ~ cauchy(0.0,5.0);  //now implicitly uniform

and

vector<lower=0,upper=20>[N_year-1] beta; // expected slopes (first is constrained to be 0.0)

which caused it to fail to parse. That's what I meant when I was talking about constraining the vector.....

That will work in the next version of Stan. It doesn't work in Stan v1.0.2. But I don't think beta needs to be a vector in the Stan sense. You could declare it now as

real<lower=0,upper=20> beta[N_year-1];

Ben

[Ben Goodrich]

On Tuesday, October 23, 2012 4:42:45 PM UTC-4, Bob Carpenter wrote:

On 10/23/12 4:32 PM, Andrew Gelman wrote:
> OK, now I have a question (forgive me if the answer's already in the manual).  Suppose I define a parameter theta to
> have a normal (mu, sigma) distribution but then theta is declared as real<lower=0, upper=0.5>.  Does Stan then work with
> the truncated normal density (i.e., getting the correct normalizing factor which is now a function of mu and sigma) or
> does it just use the normal density without adjusting for the truncation?

First, you never get the full normalization term in Stan.
All constant terms are dropped.  In the normal, we never
include the constant normalizing term -log(sqrt(2*pi)).

The difference between

    x ~ normal(mu,sigma);   // no space before '('!!!

and

    x ~ normal(mu,sigma) T[0,0.5];

is that the latter applies the truncated density and
the former does not.

If the variable x is declared as real<lower=0,upper=0.5>,
that is, with fixed bounds, then the first is OK as
we only need the density up to a proportion.  If the
bounds were variable, you'd need the second form.

Are you sure? I messed this once was in the FAQ and then Mav corrected me, but I may have messed it up again. If mu or sigma were unknown parameters, then I would think as Andrew said the "normalizing factor is now a function of mu and sigma" rather than a constant. My current attempt at the FAQ says omitting the explicit truncation is okay if the amount of truncated mass is a constant, which would be the case if the bounds were known and both mu and sigma were known.

Ben

[Andrew Gelman]

Hmm, now I'm confused. Here's what I think you're saying. Consider 3 possibilities:

1.
real theta
theta ~ normal (mu,sigma)

2.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma)

3.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma) T[0,0.5]

4.
real theta
theta ~ normal (mu, sigma) T[0,0.5]



As I understand it, option 1 does the log of the normal density (ignoring the -log(sqrt(2*pi)), which I'll ignore in the rest of this discussion).
Option 2 also uses the log of the normal density, but theta is constrained. So it's not doing the truncated normal distribution, it's doing something else (I suppose it's equivalent to the truncated distribution with a weird prior on (mu,sigma)).
Option 3 is doing the log of the truncated normal density (i.e., including the normalizing factor that depends on (mu,sigma)).
Option 4 . . . I don't know what that's doing. Perhaps it is also doing the truncated density, just with an inefficient implementation in that theta is not being subjected to a transformation for computational purposes?

> --
>

[Bob Carpenter]

On 10/23/12 5:02 PM, Ben Goodrich wrote:
> On Tuesday, October 23, 2012 4:42:45 PM UTC-4, Bob Carpenter wrote:

...

> If the variable x is declared as real<lower=0,upper=0.5>,
> that is, with fixed bounds, then the first is OK as
> we only need the density up to a proportion. If the
> bounds were variable, you'd need the second form.
>
>
> Are you sure?

No, but I just fact checked.

1. There ARE tests in place for being outside truncation bounds
that returns -inf.

I can't even remember if I added that or someone else
did -- it's only a couple lines of code.

> I messed this once was in the FAQ and then Mav corrected me, but I may have messed it up again. If mu or
> sigma were unknown parameters, then I would think as Andrew said the "normalizing factor is now a function of mu and
> sigma" rather than a constant. My current attempt at the FAQ says omitting the explicit truncation is okay if the amount
> of truncated mass is a constant, which would be the case if the bounds were known and both mu and sigma were known.

2. D'oh! It's only OK to remove the bounds
if mu and sigma are also known, because otherwise they
show up in the truncation terms!

- Bob

[Bob Carpenter]

On 10/23/12 5:08 PM, Andrew Gelman wrote:
> Hmm, now I'm confused. Here's what I think you're saying. Consider 3 possibilities:
>
> 1.
> real theta
> theta ~ normal (mu,sigma)
>
> 2.
> real<lower=0,upper=0.5> theta
> theta ~ normal (mu, sigma)
>
> 3.
> real<lower=0,upper=0.5> theta
> theta ~ normal (mu, sigma) T[0,0.5]
>
> 4.
> real theta
> theta ~ normal (mu, sigma) T[0,0.5]
>
>
>
> As I understand it, option 1 does the log of the normal density (ignoring the -log(sqrt(2*pi)), which I'll ignore in the rest of this discussion).
> Option 2 also uses the log of the normal density, but theta is constrained. So it's not doing the truncated normal distribution, it's doing something else (I suppose it's equivalent to the truncated distribution with a weird prior on (mu,sigma)).
> Option 3 is doing the log of the truncated normal density (i.e., including the normalizing factor that depends on (mu,sigma)).
> Option 4 . . . I don't know what that's doing. Perhaps it is also doing the truncated density, just with an inefficient implementation in that theta is not being subjected to a transformation for computational purposes?

That's right on all options, including 4 (the issue with
which will be illegal inits most likely).

Ben just corrected me, too.

Despite what I said earlier today, this is a great
list on which to get coaching for basic stats :-)

- Bob

[Ben Goodrich]

On Tuesday, October 23, 2012 5:08:11 PM UTC-4, Andrew Gelman wrote:

Hmm, now I'm confused.  Here's what I think you're saying.  Consider 3 possibilities:

or 4 :). My original FAQ email was a little off on this part. Read the edited version at

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/VonPUkdcZuUJ

and tell me whether it makes sense. But my understanding of your 4 possibilities are ...

1.
real theta
theta ~ normal (mu,sigma)

fine, as expected
 

2.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma)

will yield the wrong results if mu or sigma were unknown
 

3.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma) T[0,0.5]

will yield the right results, but truncation is not implemented for all distributions yet
 

4.
real theta
theta ~ normal (mu, sigma) T[0,0.5]

should be avoided because it generates proposals for theta with the log-posterior is not defined or is negative infinity.

Ben

[Jiqiang Guo]

On Tue, Oct 23, 2012 at 5:21 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Tuesday, October 23, 2012 5:08:11 PM UTC-4, Andrew Gelman wrote:

Hmm, now I'm confused.  Here's what I think you're saying.  Consider 3 possibilities:

or 4 :). My original FAQ email was a little off on this part. Read the edited version at

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/VonPUkdcZuUJ

I am not sure what the message to take from the following.  It gives me the impression that you are saying the prior distribution for real<lower=0> sigma can only be a distribution that is symmetric around zero.  But in this case, it could be any distribution whose support includes the positive real line. 

cauchy(0.0,5.0) T[0.0,] \propto cauchy(0.0,5.0)
and the same is true for other distributions that are symmetric around zero.

 

and tell me whether it makes sense. But my understanding of your 4 possibilities are ...

1.
real theta
theta ~ normal (mu,sigma)

fine, as expected
 

2.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma)

will yield the wrong results if mu or sigma were unknown
 

3.
real<lower=0,upper=0.5> theta
theta ~ normal (mu, sigma) T[0,0.5]

will yield the right results, but truncation is not implemented for all distributions yet
 

4.
real theta
theta ~ normal (mu, sigma) T[0,0.5]

should be avoided because it generates proposals for theta with the log-posterior is not defined or is negative infinity.

Ben

--
 
 

[Bob Carpenter]

On 10/23/12 7:49 PM, Jiqiang Guo wrote:

>
>
> On Tue, Oct 23, 2012 at 5:21 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> On Tuesday, October 23, 2012 5:08:11 PM UTC-4, Andrew Gelman wrote:
>
> Hmm, now I'm confused. Here's what I think you're saying. Consider 3 possibilities:
>
>
> or 4 :). My original FAQ email was a little off on this part. Read the edited version at
>
> https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/VonPUkdcZuUJ
>
>
> I am not sure what the message to take from the following. It gives me the impression that you are saying the prior
> distribution for real<lower=0> sigma can only be a distribution that is symmetric around zero. But in this case, it
> could be any distribution whose support includes the positive real line.

It can be anything. You don't need the centered Cauchy here.

you can say write "y ~ cauchy(5,2) T[0,];"

- Bob

[Sachintha Karunaratne]

Hi Bob,

Can you kindly advice whether the following have the same effect in stan (if truncation is allowed for the cauchy distribution)?

option 1:

real<lower=0> sigma; 

 sigma ~ cauchy(0,5); 

option 2:

real<lower=0> sigma; 

sigma ~ cauchy(0,5) T[0,]; 

Regards,

Sachintha.

[Bob Carpenter]

They have the same effect on parameters, but will have
different log densities by a constant because of the truncation.
That constant isn't necessary in expressing a model in Stan
because we only need the log density up to an additive constant.

They won't be the same if the truncation point on sigma is
not a constant or if the parameters to the cauchy are not constant.

- Bob

> --
> 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.
> For more options, visit https://groups.google.com/d/optout.