My model does not parse correctly

[Joe King]

Hello !

I am a newbie to stan.  I am trying to run 2 models. The first one runs in lme4/R as follows:

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

The data looks like this:

str(dt.all)

'data.frame':   30000 obs. of  4 variables:
$ pupil.ID : num  1 2 3 4 5 6 7 8 9 10 ...
$ school.ID: num  1 1 1 1 1 1 1 1 1 1 ...
$ outcome  : num  16.2 22.1 14.8 23.6 16.9 ...
$ year     : Factor w/ 3 levels "2004","2005",..: 1 1 1 1 1 1 1 1 1 1 ...

I tried to run it in rstan like this :

schools_code <- '

data {
	int<lower=1> n2 ;  
	int<lower=1> N ; 
	real y[N] ; 
	int school[N] ; 
}

parameters {
	real mu[N] ;
	real u2[n2] ;
	real<lower=0> tau ;
	real<lower=0> tau.u2 ;
	real beta ;
}

model
{
	# Level 1 definition
	for(i in 1:N) {
		y[i] ~ normal(mu[i],tau) ;
		mu[i] <- beta + u2[school[i]] ;
	}
	# level 2 definitions
	for (j in 1:n2) {
		u2[j] ~ normal(0,tau.u2) ;
	}
	# Priors for fixed effects
	beta[1] ~ uniform(0,1000) ;
	# Priors for random terms
	tau ~ dgamma(0.001000,0.001000) ;
	sigma2 <- 1/tau ;
	tau.u2 ~ dgamma(0.001000,0.001000) ;
	sigma2.u2 <- 1/tau.u2 ;
}

'

schools_dat <- list(N = 30000, n2 = 300,
                    y = dt.all$outcome,
                    school = dt.all$school.ID)

fit <- stan(model_code = schools_code, data = schools_dat, 
            iter = 1000, chains = 4)

but it does note parse. If anyone can tell me what is wrong I would be grateful.

The 2nd model is 

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

If anyone can provide some hints for how to specify this model I would be grateful.

Thanks
JK

[Ben Goodrich]

Hi Joe,

In the future, when asking about a parser error please quote the error message, which in your case is:

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

        mu[i] <- beta + u2[school[i]] ;

        ^-- here

DIAGNOSTIC(S) FROM PARSER:

attempt to assign variable in wrong block. left-hand-side variable origin=parameter
Parser expecting: "}"

This happens because mu is defined as a (free) parameter (which cannot be assigned to), rather than a transformed parameter or a local parameter in the model section. In general, your code is essentially BUGS code, and there are some differences in the Stan language that need to be kept in mind. These include: not using periods in variable names, bounding the support of a variable instead of using a uniform prior, using the standard deviation rather than the precision in the normal distribution, vectorizing instead of looping where possible, and defining a variable before using it. Most of these are mentioned in the Appendix entitled "Stan for Users of BUGS" in the Stan manual.

I'll rewrite your code so that it parses and add some comments with //.

data {
    int<lower=1> n2 ;  
    int<lower=1> N ;
    real y[N] ;
    int school[N] ;
}

parameters {

//    real mu[N] ; // moved to model{} section

    real u2[n2] ;
    real<lower=0> tau ;

//    real<lower=0> tau.u2 ;
    real<lower=0> tau_u2 ; // no periods in names
//    real beta ;
    real <lower=0,upper=1000> beta; // add bounds on the support of beta
}

model
{
    real mu; // can be scalar in this case
    real sigma;

    sigma <- sqrt(1/tau);

    # Level 1 definition
    for(i in 1:N) {

        mu <- beta + u2[school[i]] ;
        y[i] ~ normal(mu,sigma) ; // normal is parameterized in terms of the standard deviation
    }
    # level 2 definitions
//    for (j in 1:n2) {
//        u2[j] ~ normal(0,tau.u2) ;
//    }
    u2 ~ normal(0, sqrt(1/tau_u2)); // vectorize instead of looping where possible and use standard deviation

    # Priors for fixed effects

//    beta[1] ~ uniform(0,1000) ; // not necessary when the support of beta is bounded on both sides

    # Priors for random terms

    // use gamma, rather than dgamma
    tau ~ gamma(0.001000,0.001000) ;
    tau_u2 ~ gamma(0.001000,0.001000) ;
}

In addition, while it is not wrong to define the tau and tau_u2 as parameters representing if you have an informative prior over a precision, it would probably be clearer to just define sigma and sigma_u2 as nonnegative parameters in the parameters {} block (and not define tau or tau_u2) and then use sigma and sigma_u2 in the model{} block (rather than defining them as local parameters as I have done).

Ben

[Joe King]

Hi Ben

Thank you for your reply, and I am sorry for not following the correct protocol in my first post to the group.

I have implemented your changes and the model now runs.

However it seems to be very slow - so far it has taken about 30 mins for 1000 iterations (chain 1) and so far on chain 2 it has taken 40 mins and only reached 500 iterations..... This compares to about 30 seconds for 5000 iterations with mlwin (MCMC estimation mode) and 40 seconds for 10,000 iterations with MCMCglmm - with the same data.

Is there anything I can do to improve the speed ? I already set  CXXFLAGS to 03 in the makeconf

Thanks again
JK

[Bob Carpenter]

On 10/8/12 1:18 PM, Joe King wrote:> Hi Ben

>
> I have implemented your changes and the model now runs.
>
> However it seems to be very slow - so far it has taken about 30 mins for 1000 iterations (chain 1) and so far on chain 2
> it has taken 40 mins and only reached 500 iterations..... This compares to about 30 seconds for 5000 iterations with
> mlwin (MCMC estimation mode) and 40 seconds for 10,000 iterations with MCMCglmm - with the same data.
>
> Is there anything I can do to improve the speed ? I already set CXXFLAGS to 03 in the makeconf

The main thing that would speed it up is to
vectorize the normal() calls. Even if you have to use
an intermediate array it'll be faster, but only
by a factor of 2--4 or so, not a factor
of 300.

Also, be sure to measure the effective sample size,
not just the number of draws. It's easy to generate
gazillions of Metropolis draws quickly by choosing
a small proposal scale, but the effective sample size
will be lousy.

Do you know what algorithms mlwin and MCMCglmm use
for sampling?

Stan is at a disadvantage compared to systems that can exploit
sufficient statistics or calculate posteriors analytically as
can be done for many GLMs.

We're going to focus more on this kind of GLM in the
future for speed, but there's nothing on the horizon
immediately to look forward to.

- Bob

[Ben Goodrich]

Hi Joe,

On Monday, October 8, 2012 1:18:31 PM UTC-4, Joe King wrote:

However it seems to be very slow - so far it has taken about 30 mins for 1000 iterations (chain 1) and so far on chain 2 it has taken 40 mins and only reached 500 iterations..... This compares to about 30 seconds for 5000 iterations with mlwin (MCMC estimation mode) and 40 seconds for 10,000 iterations with MCMCglmm - with the same data.

Is there anything I can do to improve the speed ? I already set  CXXFLAGS to 03 in the makeconf

You can also specify -march=native in CXXFLAGS, but that won't make too much of a difference.
 
For large values of N and n2, the loop in your model{} block is going to take a while. For a model that can be estimated by Gibbs sampling, a Gibbs sampler can take less time than a HMC sampler per iteration, although hopefully Stan will require fewer iterations than would a Gibbs sampler to get a reasonable effective sample size. Maybe someone else can jump in, but I don't see anything in your code that could be sped up substantially.

Ben

[Matt Hoffman]

I agree with Bob—vectorizing the normal() calls is the easiest thing
to do, and could produce a very significant speedup. That's
implemented below. Also, it's definitely not fair to just look at
iterations—not all MCMC samples are created equal.

One more way that the model could be sped up is by reparameterizing
the model so that tau_u2 and u2 are decoupled in the prior, which is
also implemented below. When you say something like

<code>
tau ~ gamma(0.001, 0.001);
u2[i] ~ normal(mu, 1/sqrt(tau));
</code>

the result is that
1. there is a very strong correlation in the prior (and often in the
posterior) between tau and the elements of u2, which is undesirable
because it demands small step sizes and
2. marginally u2[i] has a very, very heavy-tailed prior, which can
result in a heavy-tailed posterior, which is more difficult for
NUTS/HMC (and Gibbs, and Metropolis, and pretty much any MCMC
algorithm) to sample from.

You can instead say something like

<code>
tau ~ gamma(0.001, 0.001);
u2_z[i] ~ normal(mu, 1);
u2[i] <- u2_z * 1/sqrt(tau);
</code>

and sample on the space of the transformed variable u2_z instead of
u2. This prior is equivalent to the one above, but u2_z is
uncorrelated with tau (which is good) and has normal tails (which is
also good). I've found that this sort of reparameterization can
dramatically improve performance of hierarchical models.

The model as written below should give a significant speedup, I think.
(Assuming I didn't screw it up accidentally.) But again, please
remember that what counts is the quality of the samples, not the total
number of iterations.

Matt

<code>

data {
int<lower=1> n2 ;
int<lower=1> N ;
real y[N] ;
int school[N] ;
}

parameters {

real u2_z[n2] ; // drawn from a standard normal
real<lower=0> tau ;

real<lower=0> tau_u2 ; // no periods in names

real <lower=0,upper=1000> beta; // add bounds on the support of beta
}

transformed parameters {
real u2[n2]; // this is the u2 you care about

for (i in 1:n2) {
u2[i] <- u2_z[i] * sqrt(1/tau_u2);
}
}

model
{
real mu[N]; // should be vectorized

real sigma;

sigma <- sqrt(1/tau);

# Level 1 definition
for(i in 1:N) {

mu[i] <- beta + u2[school[i]] ;

// y[i] ~ normal(mu,sigma) ; // normal is parameterized in

terms of the standard deviation
}

y ~ normal(mu, sigma); // vectorize instead of looping

# level 2 definitions

u2 ~ normal(0, sqrt(1/tau_u2)); // vectorize instead of looping
where possible and use standard deviation

# Priors for random terms

// use gamma, rather than dgamma

tau ~ gamma(0.001,0.001) ;
tau_u2 ~ gamma(0.001,0.001) ;
}
</code>

> --
>
>

[Douglas Bates]

I don't think I would call this model a GLM.  I would say it is a hierarchical linear model (MLM) or a linear mixed model (LMM).  One could refer to it as a generalized linear mixed model but only in the sense that all linear mixed models are also generalized linear mixed models.

I'm not sure how easy (well it is more a question of whether or not it is exceptionally difficult) it would be to extend the sampler framework to take advantage of models with a linear predictor expression like this.  Conditional on the values of sigma2 and sigma2.u2 the distribution of beta and mu is a multivariate normal and can be vectorized except that you really want to use a sparse model matrix, not a dense model matrix, when you get to large numbers of observations.

[Joe King]

Hi Bob, and thanks again for your time !

As far as I know, the default in MCMCglmm is an adaptive Metropolis-Hastings algorithm. I'm not sure about mlwin (I don't use it very often and only dusted it off today to compare this model).

I am beginning to wonder if the models are equivalent. The stan model has resulted in a more output than I am used to seeing with MCMCglmm or other glmm packages (such as lme4).  Normally I expect a single variance estimate for the random intercept (school) and the residual plus an estimate for the overall mean (intercept). The intercept estimate from stan (beta) corresponds with the output I get with lme4/MCMCglmm, although the effective sample size for beta is only 3

.For the random intercept stan doesn't seem to produce the variance for it - instead I get 300 values u2[1] to u[300].  The effective samples for these seem quite good (1400-1500) but I don't understand what they are !

It also produces estimates of 0.0 for tau and tau_u2. 

FWIW this is the output I get from lme4, which I trust because these data are artificial: I simulated them (I could send you the code if you like ?) so I know what estimates are reasonable:

Linear mixed model fit by REML
Formula: outcome ~ 1 + (1 | school.ID)
........
Random effects:
 Groups    Name        Variance Std.Dev.
 school.ID (Intercept)  34.616   5.8835
 Residual              105.751  10.2835
Number of obs: 30000, groups: school.ID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)   20.669      0.591   34.97

Any further insight would be appreciated very much.

Thanks again !
JK

[Joe King]

Oh, I guess the u2[1] to u2[300] are the individual school effects and I just have to compute the variance of them ?

But it's strange that there are only 100 unique schools.

[Joe King]

Whoops. My bad, n2 should be 100......

[Ben Goodrich]

On Monday, October 8, 2012 2:13:58 PM UTC-4, Joe King wrote:

Oh, I guess the u2[1] to u2[300] are the individual school effects and I just have to compute the variance of them ?

But it's strange that there are only 100 unique schools.

The length of u2 is whatever n2 is. See Matt Hoffman's comments especially; I think that could make a big difference. If you want to store the standard deviation of u2, you can put sigma into the transformed parameters {} block instead of the model{} block like

transformed parameters {
  real<lower=0> sigma_u2;
  sigma_u2 <- sqrt(1/tau_u2);
}

or define sigma_u2 in the parameters {} block and drop tau_u2.

Ben

[Douglas Bates]

On Monday, October 8, 2012 1:13:58 PM UTC-5, Joe King wrote:

Oh, I guess the u2[1] to u2[300] are the individual school effects and I just have to compute the variance of them ?

Not really.  The observed variance of the conditional means of the random effects (what are sometimes called the BLUPs, best linear unbiased estimators) that are available in an lmer fitted model as

ranef(myLmerModel)

is not the same as the estimate of the variance of the random effects.

[Joe King]

Hi Ben

Thank you, although I am sorry that I don't understand. With the code you wrote below, u2 is then undefined and it doesn't parse.

Sorry if I'm being stupid......
JK

[Joe King]

Thank you Douglas - I think I understand that now. I am familiar with lmer. Actually I am trying to produce output similar to what I get with this:

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

Linear mixed model fit by REML
Formula: outcome ~ 1 + (1 | school.ID)

   Data: dt.all
    AIC    BIC  logLik deviance REMLdev
 225433 225458 -112713   225428  225427

Random effects:
 Groups    Name        Variance Std.Dev.
 school.ID (Intercept)  34.616   5.8835
 Residual              105.751  10.2835
Number of obs: 30000, groups: school.ID, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)   20.669      0.591   34.97

Do you know how I can achieve that with stan ?

Thanks again
JK

[Ben Goodrich]

On Monday, October 8, 2012 2:41:50 PM UTC-4, Joe King wrote:

Thank you, although I am sorry that I don't understand. With the code you wrote below, u2 is then undefined and it doesn't parse.

Can't quite tell what you mean. It is always easier if you just attach the file that does or doesn't parse. Here is what I would do, approximately incorporating Matt and Bob's advice (also attached):

data {
    int<lower=1> n2 ;  
    int<lower=1> N ;
    real y[N] ;
    int school[N] ;
}
parameters {

    real u2_z[n2] ;
    real<lower=0> sigma ;
    real<lower=0> 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) ;

    // or some other prior on sigma and sigma_u2
    sigma ~ gamma(0.001000,0.001000) ;
    sigma_u2 ~ gamma(0.001000,0.001000) ;
}

Ben



[Ben Goodrich]

On Monday, October 8, 2012 3:13:02 PM UTC-4, Ben Goodrich wrote:

On Monday, October 8, 2012 2:41:50 PM UTC-4, Joe King wrote:

Thank you, although I am sorry that I don't understand. With the code you wrote below, u2 is then undefined and it doesn't parse.

Can't quite tell what you mean. It is always easier if you just attach the file that does or doesn't parse. Here is what I would do, approximately incorporating Matt and Bob's advice (also attached):

The attachment didn't correspond to the code in the post. This attachment is the right one.

Ben
 

[Joe King]

Hi Ben

Thanks !  It has parsed now.  And it seems very much faster, However, it gives the following warning several times during iterations:

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!

Also from the final output from print(fit) sigma and sigma_u2 have an n_eff of only 2 and beta of only 15

Please advise

Thanks again
JK

[Joe King]

Also some of these:

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Location parameter[0]  is -1.#INF:0, but must be finite!

[Ben Goodrich]

On Monday, October 8, 2012 3:49:38 PM UTC-4, Joe King wrote:

Thanks !  It has parsed now.  And it seems very much faster, However, it gives the following warning several times during iterations:

Warning (non-fatal rejection): Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is 0:0, but must be > 0!

A few of those, particularly early in the run, are not so bad. But ...
 

Also from the final output from print(fit) sigma and sigma_u2 have an n_eff of only 2 and beta of only 15

I'm guessing that sigma and sigma_u2 are close to zero? That may mean that there is not much point in having a multilevel model, or it may mean Stan is having a hard time sampling, or it may mean something else. I would try bounding sigma and sigma_u2 away from zero with something like

parameters {
  real<lower = .1, upper = 100> sigma;
  real<lower = .1, upper = 100> sigma_u2;
  ...
}

and comment out the priors on sigma and sigma_u2 in the model{} block. Then run it again. If the posterior mass piles up near .1, then maybe the error variances are just small.

Ben



[Joe King]

No, they are not close to zero. I simulated these data so I know what the values "should" be - for this dataset they should be around 10 and 5 respectively

With your latest amendments, the warnings are no more, but the 95% credible intervals from stan are (5.0, 5.9) and (0.1, 1.1) with n_eff of 21 and 2 respectively.

The 95% CI for beta is , (19.7 , 21.1) with n_eff of 246 which compares well to the fixed value I used to simulate the data of 20....

[Andrew Gelman]

You've heard this before from me, but I'll say it again: I think gamma(.001,.001) is bad news.
A

[Douglas Bates]

The REML estimates of the standard deviations in the lmer fit are 5.9 and 10.3 for between and within groups, respectively.  My guess is that the between groups standard deviation (sigma_u2) is very poorly determined, which may be causing the chain to wander.  When the between groups standard deviation is dominated by the residual variability like this, it is very hard to estimate it.

Joe, would you be willing to make the data available somewhere so that we could check on some of these issues?

 

[Ben Goodrich]

OK, that suggests something weird is going on with the sampling. Can you attach the R code that generates the data?

Ben

[Andrew Gelman]

When the between-groups sd is low (as in the 8-schools example), Matt's reparameterization (for an example, see here: http://code.google.com/p/stan/wiki/RStanGettingStarted#Example_1:_Eight_Schools  )

should help a lot.

On Oct 8, 2012, at 4:32 PM, Douglas Bates wrote:

The REML estimates of the standard deviations in the lmer fit are 5.9 and 10.3 for between and within groups, respectively.  My guess is that the between groups standard deviation (sigma_u2) is very poorly determined, which may be causing the chain to wander.  When the between groups standard deviation is dominated by the residual variability like this, it is very hard to estimate it.

Joe, would you be willing to make the data available somewhere so that we could check on some of these issues?

 

--
 
 

[Joe King]

Hi Douglas

Yes ! I can make the data available, because I simulated it myself. I shall post it here, below, but first I would like to explain a little about it so you can understand what I did.

In my "real-life" data (which I have not analysed with stan yet) I have observations on pupils, nested in schools. There is 1 observation per pupil in the data. Moreover, the data are over several years. Interest is centred on how the school-level random effects have changed over time.

So, I simulated 3 years worth of data, where each year had a different variance for the school-level effects and I wanted to be able to "get back" those different variances in my model. I was able to do that very well, using lme4 with this:

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

So my ultimate goal is to run a model with stan which does the same. But first I am just trying out the simple variance components model, for which I gave the output for earlier.

So, this is the code that generates the data:

set.seed(10)
n.schools <- 100
n.pupils <- n.schools^2

pupil.ID <- seq(1,n.pupils,1)     # pupils with unique IDs
school.ID <- 1+trunc((pupil.ID-1)/n.schools) # schools with unique IDs
school.effect <- rnorm(n.schools,0,1)  # school-level effects, SD=1
outcome <- numeric(n.pupils)       
for (i in 1:n.pupils) {   
    outcome[i] <- 20+school.effect[school.ID[i]]+rnorm(1,0,5)
}

dt.y1 <- as.data.frame(cbind(pupil.ID, school.ID, outcome))

# Year 2
dt.y1$year <- rep(2004,n.pupils)  # year ID for first period

pupil.ID <- seq(n.pupils+1,n.pupils*2,1)     # add more pupils with unique IDs
school.effect <- rnorm(n.schools,0,4)  # school-level effects - SD=4
for (i in 1:n.pupils) { 
    outcome[i] <- 20+school.effect[school.ID[i]]+rnorm(1,0,5)
}
dt.y2 <- as.data.frame(cbind(pupil.ID, school.ID, outcome))
dt.y2$year <- rep(2005,n.pupils)  # year ID for first period


# Year 3

pupil.ID <- seq((n.pupils*2)+1,n.pupils*3,1)     # add more pupils with unique IDs
school.effect <- rnorm(n.schools,0,16)  # school-level effects, SD=16
for (i in 1:n.pupils) { 
    outcome[i] <- 20+school.effect[school.ID[i]]+rnorm(1,0,5)
}
dt.y3 <- as.data.frame(cbind(pupil.ID, school.ID, outcome))
dt.y3$year <- rep(2006,n.pupils)  # year ID for first period

fit.lme4.y1 <- lmer(outcome ~1 + (1|school.ID),data=dt.y1)
summary(fit.lme4.y1)

fit.lme4.y2 <- lmer(outcome ~1 + (1|school.ID),data=dt.y2)
summary(fit.lme4.y2)

fit.lme4.y3 <- lmer(outcome ~1 + (1|school.ID),data=dt.y3)
summary(fit.lme4.y3)

dt.all <- rbind(dt.y1,dt.y2,dt.y3)
fit.lme4.all.1 <- lmer(outcome ~1 + (1|school.ID),data=dt.all)
summary(fit.lme4.all.1)

dt.all$year <- as.factor(dt.all$year)
fit.lme4.all.2 <- lmer(outcome ~1 + (year|school.ID),data=dt.all)
summary(fit.lme4.all.2)

Thanks !
JK

[Joe King]

In another project I had problems with gamma(.001,.001) using MCMCglmm and instead I used a parameter-expanded representation of a half-cauchy priory with better results. Do you think that would be appropriate here ? Unfortunately I didn't fully understand how the the scaling worked.

[Ben Goodrich]

For command-line Stan users, a joe.stan and a joe.Rdata file are attached, so that

> make /path/to/joe

should work from the STAN_HOME directory. It seems that the stepsize__ is small and not adapting for some reason:

> summary(samples[,!grepl("u2_z", colnames(samples))])
      lp__         treedepth__      stepsize__           sigma      
 Min.   :-21258   Min.   :4.000   Min.   :0.006126   Min.   :4.941  
 1st Qu.:-21230   1st Qu.:4.000   1st Qu.:0.006126   1st Qu.:5.019  
 Median :-21224   Median :4.000   Median :0.006126   Median :5.046  
 Mean   :-21224   Mean   :4.226   Mean   :0.006126   Mean   :5.046  
 3rd Qu.:-21217   3rd Qu.:4.000   3rd Qu.:0.006126   3rd Qu.:5.074  
 Max.   :-21200   Max.   :7.000   Max.   :0.006126   Max.   :5.172  
    sigma_u2          beta      
 Min.   :14.85   Min.   :18.60  
 1st Qu.:17.10   1st Qu.:21.13  
 Median :17.89   Median :22.15  
 Mean   :18.11   Mean   :21.81  
 3rd Qu.:19.13   3rd Qu.:22.60  
 Max.   :21.08   Max.   :23.90  

Ben

[Bob Carpenter]

> On Monday, October 8, 2012 9:29:05 PM UTC+1, Andrew Gelman wrote:
>
> You've heard this before from me, but I'll say it again: I think gamma(.001,.001) is bad news.

On 10/8/12 5:09 PM, Joe King wrote:
> In another project I had problems with gamma(.001,.001) using MCMCglmm and instead I used a parameter-expanded
> representation of a half-cauchy priory with better results. Do you think that would be appropriate here ? Unfortunately
> I didn't fully understand how the the scaling worked.

You can read up on Andrew's full set of recommendations here:

Bayesian Analysis (2006) 1, Number 3, pp. 515�533
Prior distributions for variance parameters in
hierarchical models
Andrew Gelman

http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf

We've been using half-Cauchy distributions for deviations and
full Cauchy for regression coefficients. In
Stan, you don't even need to truncate because we only need
the probability up to a proportion. So you could replace this:

> tau ~ gamma(0.001, 0.001);
> u2_z[i] ~ normal(mu, 1);
> u2[i] <- u2_z * 1/sqrt(tau);

with a proper z-score and a prior directly on deviation:

sigma ~ cauchy(0,5);
u2_z[i] ~ normal(0,1);
u2[i] <- mu + u2_z[i] * sigma;

I also pulled mu out for reasons Matt outlined in previous
messages (and which is going into the next manual).

Because Stan doesn't exploit conjugacy, the Cauchy will
be roughly the same efficiency per iteration as the gamma
prior.

Our standard advice to do some sensitivity analysis also
stands -- if you don't want the prior to have a big effect,
make sure it doesn't by testing alternatives.

- Bob

[Bob Carpenter]

The default in Stan is to adapt step size during warmup,
but then to use a fixed step size in sampling. So you
see the same step size in every sample.

If you use the varying step command-line flag, epsilon_pm,
you get different step sizes per iteration. The "_pm" is
missing from the doc in the manual, so I just fixed that.
Here's how it should read:

Model execution option (page 23 or so of the manual):

--epsilon_pm=<[0,1]>
Sample epsilon uniformly from the interval
[epsilon * (1 - epsilon_pm), epsilon * (1 + epsilon_pm)]
(default = 0.0)


- Bob

[Ben Goodrich]

Ah, I had forgotten that I was only looking at the post-warmup output.

But we still have the question of why sigma_u2 is far from the true value with n_eff in the single digits, which seems to persist whether the priors on the standard deviations are half-Cauchy, bounded uniform, or gamma.

Ben

[Matt Hoffman]

Is sigma_u2 really so far from the true value? It looks to me like its
true value (16) is within the range of the simulations.

As for why the ESSs are on the low side, I have a hypothesis that I'll
test and report on shortly.

Matt

> --
>
>

[Ben Goodrich]

On Monday, October 8, 2012 6:01:33 PM UTC-4, Matt Hoffman wrote:

Is sigma_u2 really so far from the true value? It looks to me like its
true value (16) is within the range of the simulations.

You're right. I just looked at the top of the file where the sigma_u2 was 1 and didn't realize that it had changed to 16 when I dumped the data.

Ben

[Matt Hoffman]

Never mind about that hypothesis I mentioned.

The reason for the slow mixing we're seeing is pretty basic—in this
model, the intercept beta (just out of curiosity—why call it beta
instead of alpha?) is very highly anti-correlated with the
school-level effects u2. On average, the posterior correlation is
around -0.9. So that's the source of the slow mixing.

Matt

[Douglas Bates]

On Monday, October 8, 2012 5:19:55 PM UTC-5, Matt Hoffman wrote:

Never mind about that hypothesis I mentioned.

The reason for the slow mixing we're seeing is pretty basic—in this
model, the intercept beta (just out of curiosity—why call it beta
instead of alpha?) is very highly anti-correlated with the
school-level effects u2. On average, the posterior correlation is
around -0.9. So that's the source of the slow mixing.

One way around this is to sample from the coefficient vector, which in this case can be orthogonalized.  In general the conditional distribution of coefficients, even in the case of mixed-effects models can be orthogonalized.

IIRC JAGS has a specialized sampler for this.

I'm not sure how difficult it would be but the ability to use sparse matrices, which are supported in Eigen and RcppEigen, would be very helpful here.

[Douglas Bates]

For the record, I would use something like

library(lme4)

library(rstan)

set.seed(1234321)

samp <- within(data.frame(schoolID=gl(300,100)),

               outcome <- 20 + rnorm(30000, sd=10) + rnorm(300,sd=5)[schoolID])

(lm1 <- lmer(outcome ~ 1|schoolID, samp))

## The next line can only be run with the development lme4 installed as

## install.packages("lme4", repos="http://R-forge.R-project.org")

densityplot(pr1 <- profile(lm1), ylab=NULL, layout=c(3,1))

schools_code <- '

  data {

    int<lower=1> N;                   // number of observations

    int<lower=1> J;                   // number of schools

    int<lower=1,upper=J> school[N];   // school for each obs

    real         y[N];          // observed outcome

  }

  parameters {

    real mu; 

    real<lower=0> sigma;

    real<lower=0> sigma1;

    real eta[J];

  }

  model {

    real mm[N];

    for (i in 1:N) 

      mm[i] <- mu + eta[school[i]];

    y   ~ normal(mm, sigma);

    eta ~ normal(0, sigma1);

    sigma ~ gamma(0.001000,0.001000) ;

    sigma1 ~ gamma(0.001000,0.001000) ;

  }

'

fit <- stan(model_name="schools", model_code=schools_code,

            data=with(samp, list(N=nrow(samp), J=length(levels(schoolID)),

            y=outcome, school=as.integer(schoolID))),

            iter = 1000, chains = 4)