Difficulties with lkj_corr

[Markian Pahuta]

Hi,

I am trying to perform a meta-analysis of longitudinal studies using the approach described in http://www.ncbi.nlm.nih.gov/pubmed/19254930

My model data and code are available at https://drive.google.com/folderview?id=0B2_rKFnvrjMATG5JaVpuVjlPUVk&usp=sharing

Not all studies have observations at all time periods, so my model is clunky:

model{
  matrix[Ntm,Ntm] L_beta_0_prior;
  matrix[Ncovar,Ncovar] L_betas_prior;
 
  #impute missing standard deviations  
  s_obs~cauchy(0,scl);
  s_ms~cauchy(0,scl);


  #hyper priors for intercept and coefficients
  beta_0_prior~multi_normal_cholesky(hyper_beta_0,L_beta_0_hyper);
  betas_prior~multi_normal_cholesky(hyper_betas,L_betas_hyper);


  #heterogeneity parameter for random effect for intercept and coefficients
  tau_beta_0_sing ~ inv_gamma(gam_shape,gam_rate);
  tau_betas_sing ~ inv_gamma(gam_shape,gam_rate);


  #correlation matrices for prior and study
  Omega_prior ~ lkj_corr(1);
  Omega_study ~ lkj_corr(1);




  L_betas_prior<-cholesky_decompose(diag_matrix(rep_vector(sqrt(tau_betas_sing),Ncovar)));
  L_beta_0_prior <- cholesky_decompose(Omega_prior);


  for (i in 1:Ntm){
    for (j in 1:i){
      L_beta_0_prior[i,j] <- sqrt(tau_beta_0_sing) * L_beta_0_prior[i,j];
    }
  }


  #study level intercept and coefficients
  for(study in 1:Nstudy){
    beta_0[study]~multi_normal_cholesky(beta_0_prior,L_beta_0_prior);
    betas[study]~multi_normal_cholesky(betas_prior,L_betas_prior);
  }


  for(grp in 1:Ngrp){
    matrix [len_tms[grp],Ntm] des_mat;
    vector [len_tms[grp]] temp1;
    vector [len_tms[grp]] temp2;
    matrix [len_tms[grp],len_tms[grp]] temp3;
    matrix [len_tms[grp],len_tms[grp]] temp4;
   
    #Create a design matrix for each study, so that only observed time points are used
    des_mat<- rep_matrix(0,len_tms[grp],Ntm);
    for(i in 1:len_tms[grp]){
      des_mat[i,tms[grp,i]]<-1;
    }


    temp1<-des_mat*mns[grp];
    temp2<-des_mat*beta_0[study_ind[grp]];
    temp3<- des_mat*covar_grp[grp]*des_mat';
    temp4<-cholesky_decompose(temp3);
    temp1~multi_normal_cholesky(temp2,temp4);
  }
}

My model runs when the parameters Omega_prior and Omega_study are a diagonal matrix of 1 (uncomment variable definition in data block, and comment in parameter block). 

When I run the model using the lkj distribution I get the error:

Error in function stan::prob::lkj_corr_log(N4stan5agrad3varE): yy is not positive definite. y(0,0) is 1:0.Rejecting proposed initial value with zero density.

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

When I define the Omega_ parameters simply as matrices (not corr_matrix), I get the error:

Error in function stan::prob::lkj_corr_log(N4stan5agrad3varE): yy is not symmetric. y[0,1] is 1.2277783349104863:0, but y[1,0] element is -1.2975:0Rejecting proposed initial value with zero density.

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

I have gone over my code several times, and tried to replicate what is outlined in the STAN manual... but I can't get it to work.  I was hoping someone could point me in the right direction?

Thanks

Mark

[Ben Goodrich]

On Tuesday, October 14, 2014 9:33:32 PM UTC-4, Markian Pahuta wrote:

My model runs when the parameters Omega_prior and Omega_study are a diagonal matrix of 1 (uncomment variable definition in data block, and comment in parameter block). 

When I run the model using the lkj distribution I get the error:

Error in function stan::prob::lkj_corr_log(N4stan5agrad3varE): yy is not positive definite. y(0,0) is 1:0.Rejecting proposed initial value with zero density.

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

You can circumvent this error by commenting out

Omega_prior ~ lkj_corr(1);
Omega_study ~ lkj_corr(1);


if your shape parameter is 1, then this prior is uniform over positive definite correlation matrices. So it doesn't affect the posterior in principle, although the function has to check whether the matrix is in fact positive definite, which is numerically fallible even though a declared correlation matrix will be positive definite by construction in Stan. However, I would guess that this is a symptom of a problem with your .stan program, rather than the problem itself.

Probably unrelated, but we do want to push people toward expressing their models in terms of a Cholesky factor of a correlation matrix, rather than declaring a correlation matrix in the parameters block. That avoids a lot of numerical problems and can be a lot faster. See the cholesky_factor_corr type and the lkj_cholesky_corr(1) density, because you do have to register the prior explicitly when dealing with the Cholesky factor.
 

When I define the Omega_ parameters simply as matrices (not corr_matrix), I get the error:

Error in function stan::prob::lkj_corr_log(N4stan5agrad3varE): yy is not symmetric. y[0,1] is 1.2277783349104863:0, but y[1,0] element is -1.2975:0Rejecting proposed initial value with zero density.

Initialization failed after 100 attempts.  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

That won't work for the reason you encountered: It isn't a symmetric matrix so it can't be a correlation matrix.

Ben

[Bob Carpenter]

The declarations need to be as corr_matrix, but even so,
Stan will run into trouble in high dimensions due to the
transform from unconstrained parameters and/or the
solver. If your correlation matrices are less than
10 dimensional, I don't see why it wouldn't at least
initialize as you wrote it (given that it works with
hard-coded unit correlation matrices).

But having said all that, what you really want to do for something
like this is follow Ben Goodrich's advice, which I've
written up for the 2.5 manual (attached). The final coding we recommend is
in the "Optimization through Cholesky Factorization" section.
I just realized in first introducing the model I should point
to this implementation, as it should be a lot faster and a lot
more stable arithmetically.

We're also not recommending very diffuse (inverse) gamma
priors like BUGS has in its example models, but I don't know
what your gam_shape or gam_rate values are. I'm just writing
up Andrew's advice on those in another section for the 2.5 manual.

- Bob

P.S. And if you haven't built up to this complicated model through
simpler models, I'd suggest doing that.

[Markian Pahuta]

Hi,

Thanks for your tips. I have built the model up from smaller chunks, and it works properly until I reach this point.

My correlation matrix is 16x16... could this be the problem?

> --
> 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.


--
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/eFMxEuV_ENA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

--- Mark

[Ben Goodrich]

On Wednesday, October 15, 2014 7:34:35 AM UTC-4, Markian Pahuta wrote:

Thanks for your tips. I have built the model up from smaller chunks, and it works properly until I reach this point.

My correlation matrix is 16x16... could this be the problem?

It could, although I would still wager on symptom of the problem. Either way, no one should ever use a corr_matrix that is declared in the parameters block; just use cholesky_factor_corr but be sure to put an lkj_cholesky_corr prior on it in the model block.

Ben

[Markian Pahuta]

Thank you for your tips. I got the model to work using the cholesky approach.

an additional question, if 

 matrix <lower=0> [Ngrp,Ntm] se_mat;

then is se_mat[grp] a column vector of Ntm items? Ie this will be all elements in the row indexed by grp?

I used print and it seems to confirm this, but I want to make sure.

thanks

mark

--

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/eFMxEuV_ENA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

--- Mark

[Bob Carpenter]

If a is a matrix, then a[n] is the n-th *row* vector.
That's why we can write a linear regression as:

y ~ normal(x * beta, sigma);

or as

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

where y is a vector, x the data matrix, and beta the
coefficient vector.

- 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.

[Markian Pahuta]

Hi, 

a follow-up question to the advice regarding the cauchy distribution for the standard deviation. I have been specifying my informative inverse gamma prior using the approach outlined in this paper: http://www.ncbi.nlm.nih.gov/pubmed/8619108 (bottom of page 5)

essentially, I find an IG distribution with mean as my plausible estimate of variance, and less than 5% probability of exceeding my high estimate of variance.

are there any guidelines on developing a similar informative prior with the cauchy distribution?

[Bob Carpenter]

[new title]

Andrew should really answer this one.

You could do the same thing with the Cauchy using the median
and inverse CDF (be careful to account for truncation which
doubles the tail area). I'm not sure how different that would be
from the inverse gamma in practice. Probably not much if the data
supports the value you're expecting.

- Bob

[Andrew Gelman]

Hi, I don’t think the inverse-gamma distribution works well, for reasons discussed in my 2006 paper. But, whatever prior you set, I think it probably makes sense to set it a bit more vague than from your plausible estimate. I have some general idea that the prior is there to regularize the inferences, to rule out “bad” values more than to bring the estimate toward “good” values. I do not yet have a precise formulation of this idea but it’s how I think about things.
A

[Markian Pahuta]

Hi,

so would a reasonable approach be to set the scale parameter on the cauchy to my high plausible estimate (ie 95%)?

[Andrew Gelman]

Yes, I think that would make sense.