Initialization failed: When trying to implement a LONGITUDINAL RANDOM EFFECTS MODELS WITH MEASUREMENT ERROR

hi all.

I would like to extend the LONGITUDINAL RANDOM EFFECTS MODELS WITH MEASUREMENT ERROR 1) to a Bayesian framework.
In the following master thesis 2), it is actually done and I want to reproduce it, but it does not converge well.
Could you please help me? please…

  1. Buonaccorsi, John, et al. “ESTIMATION IN LONGITUDINAL RANDOM EFFECTS MODELS WITH MEASUREMENT ERROR.” Statistica Sinica , vol. 10, no. 3, 2000, pp. 885–903. JSTOR , Accessed 8 July 2021.

  2. Yuan, Zheng . Linear Mixed Effects Models with Measurement Error: Bayesian Approach.;jsessionid=1F0C071F0C076887AA5EB52BFA599B58?sequence=6

  data {                                                                        
        int N; // num of obs (instance * year)                                                                    
        int J; // num of instance 
        int Y; // num of year                                                                    
        int year[N]; // year vector                                                                 
        vector[N] w; // error in variable 
        vector[N] y; // outcome                                                               
        matrix[5, 2] A; // desigh matrix 

     parameters {                                                                                                                           
       real  beta;                                                              
       real alpha;
       vector[1] epsiron;
       real<lower=0> s_epsiron;             
       vector[Y] x;
       vector[2] alpha_x;
       vector<lower=0>[2] phy_x;
       corr_matrix[2] corr_phy;
       vector[1] me;
       real <lower=0> s_me;
     transformed parameters {
        cov_matrix[2] cov_phy;
        matrix[Y, Y] x_cov;
        matrix[Y, Y] w_cov;
        matrix[Y, Y] y_cov;
        cov_matrix[Y] me_cov;
        cov_matrix[Y] epsiron_cov;
        vector[2] x_variables;
        cov_phy = quad_form_diag(corr_phy, phy_x);
        x_cov = A * cov_phy * A';
        w_cov = x_cov + me_cov;
        y_cov = beta^2 * x_cov + epsiron_cov;
        x_variables = alpha_x + phy_x;
        me_cov = s_me * diag_matrix(rep_vector(1,Y));
        epsiron_cov= s_epsiron * diag_matrix(rep_vector(1,Y));

     model {
       corr_phy ~ lkj_corr(2);
       phy_x ~ cauchy(0, 2.5);
       me~ normal(0, s_me);
       epsiron~ normal(0, s_epsiron);
       x ~ multi_normal(A * x_variables,  x_cov);  //exposure model  
       w ~ multi_normal(x[year], w_cov[year, year]);    //measurement error model                                                                                                           
       y ~ multi_normal(alpha + x[year] * beta, y_cov[year, year]); //disease model

df = read.csv('sample_1d_data.csv')

a_vec = c(1, 1, 1, 1, 1, 0, 1, 2, 3, 4)
A = matrix(a_vec, nrow=5)

stanfile <- 'stancode.stan'

data <-list(

error code

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is -6.66134e-16.  (in 'model354c455069e1_stancode_matrix' at line 49)

Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."  

Hi, sorry for not getting to you earlier, your question is relevant and well written.

The problem appears to be that the way you construct the covariance matrices does not ensure they are valid (I.e. positive definite). I sm on a phone so can’t run the code myself, but the indexing with year looks a bit suspicious as it will repeatwhole rows and columns from the smallercovariancematrixtocreate a big one and that is unlikely to result in a valid cov matrix . What are you trying to achieve there? Maybe you actually want to have multiple vectors office Y each month
MVN distributed? In any case a good way to debug would be to use print statements to see the intermediate results and check where exactly do the problematic values arise.

Also, it is generally substantially more efficient to working Cholesky decomposition of the matrices, so if you can rewrite your matrix algebra to work there from the start, it could help to make the model faster (but that’s a second-order concern)

Best of luck with your model!