Multivariate, hierarchical ordered probit mixture model - sharing nuisance parameters

I have developed some multivariate probit models - based on bgoodri’s parameterization of the MVP - to meta-analyse diagnostic test accuracy data without a gold standard, where studies report data at different thresholds.

For my datasets which do not have patient-level covariates, individuals with the same test response patterns contribute equally to the likelihood. This means I can assign the same latent vector across these individuals, and hence they can also share the same nuisance parameters and increase efficiency. These models are really slow in Stan, so this is important, and i don’t yet have the expertise to code my own sampler for a model this complex (such as calibrated data augmentation Gibbs sampling, which would be much faster)

Sharing the same nuisance parameters works fine with no correlation, and gives the same answers as the version which estimates a separate nuisance parameter for every subject. The problem is, sharing the same nuisance parameters doesn’t work when I model the correlation

I’ll show the code for models which assume correlation = 0 first. The difference is in the parameters “u_d” and “u_nd”. Basically, the model below (model #1 ) generates a nuisance parameter for every observation, whereas model #2 only generates num_patterns*nt nuisance parameters (where nt = numer of tests = dimension, and num_patterns = total number of possible test patterns).

Model 1:


data {
          int num_binary_tests;
          int<lower=1> nt; // number of tests
          int<lower=1> NS;	// total # of studies
          int<lower=1> ns[NS]; // number of individuals in each study 
          int<lower=0> y[max(ns), nt, NS]; // N individuals and nt tests, NS studies 
          int num_patterns;
          int<lower=0> pa[max(ns), NS]; 
          int numg;
          int num_thresh_params[NS];
          int max_cutpoint;
          int cutpoints[NS, max_cutpoint];
          int twoway_thr[NS];
          int<lower=0, upper=1> primary_care_ind[NS]; //indicator for whether the study comes from primary care or not
          int r[NS, 4, choose(nt,2)];
          int num_refs; // number of reference tests
          int ref[NS]; // ref. test used in each of the NS studies
          int re[NS];
}

transformed data {
              int r2[choose(nt,2), NS, 4];
  
       for (a in 1:4)
          for (s in 1:NS)
            r2[1,s,a] = r[s,a,1];
}

parameters {
              ordered[2] a1_m_raw[num_refs]; // Se > Fp for each reference test
              ordered[2] a2_m_raw;   // don't order this due to potential thr. interction
              ordered[max_cutpoint] C_dm;
              ordered[max_cutpoint] C_ndm;
              vector<lower=0>[2] sd1;  
              real<lower=0,upper=1> u_d[max(ns),nt, NS]; // nuisance that absorbs inequality constraints
              real<lower=0,upper=1> u_nd[max(ns),nt, NS]; // nuisance that absorbs inequality constraints
              cholesky_factor_corr[2] L_Omega_bs1;
              vector[2] z2[NS];
              vector[NS] z3;
              vector[NS] z4;
              real<lower=0> sd3;
              real<lower=0, upper=1> p[NS]; 
              vector[2] b_primary_care; // coefficients for primary care (1 for MMSE Se and 1 for MMSE Sp)
}

transformed parameters {
              vector[2] mu[num_refs + 1];  
              matrix[2, nt] nu[NS];
              matrix[NS, num_patterns] log_lik_i; 
              vector[max_cutpoint] C_d[NS];
              vector[max_cutpoint] C_nd[NS];
              vector[num_patterns] lp1;
              vector[num_patterns] lp0;
              vector[nt] z_d[num_patterns];
              vector[nt] z_nd[num_patterns];
              vector[nt] y1d[num_patterns];
              vector[nt] y1nd[num_patterns];

                          // index test (MMSE)
                          mu[num_refs + 1, 1] =  a2_m_raw[1];
                          mu[num_refs + 1, 2] =  a2_m_raw[2];
                          
              for (s in 1:NS) {
                          // ref test (Se > Fp)
                          mu[ref[s], 1] =  a1_m_raw[ref[s], 2];   
                          mu[ref[s], 2] =  a1_m_raw[ref[s], 1];

                 // between-study model
                     // fixed ref. tests
                  nu[s, 1:2 , 1] =  mu[ref[s], 1:2]  ;  


                      // index test random effect - w/ covariate for primary care
                   nu[s, 1:2 , 2] =  mu[num_refs + 1 , 1:2 ]   + primary_care_ind[s]*b_primary_care[1:2] 
                                    + diag_pre_multiply(sd1, L_Omega_bs1) * z2[s,1:2] ; 

                for (i in 1:max_cutpoint) { 
                   C_d[s,i] = C_dm[i]  ; 
                   C_nd[s,i] = C_ndm[i]; 
                }
              }

 // Parameters for likelihood function. Based on code upoaded by Ben Goodrich which uses the 
                      // GHK algorithm for generating TruncMVN. See:
                      // https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan#L11
                      {
                 for (s in 1:NS) {
                    for (n in 1:ns[s]) {
                      for (i in 1:nt) { // loop over each test
                           real u_d1 = u_d[n,i,s];
                           real u_nd1 = u_nd[n,i,s];
                           real bound_d_bin     = phi_logit_approx(    - (nu[s,1,i])   ); // se
                           real bound_nd_bin    = phi_logit_approx(    - (nu[s,2,i])   );  // fp
                           vector[ num_thresh_params[s] + 1 ] bound_d;
                           vector[ num_thresh_params[s] + 1 ] bound_nd;
                           bound_d[ num_thresh_params[s] + 1 ] = 1;
                           bound_nd[ num_thresh_params[s] + 1 ] = 1;
                           bound_d[1] =  phi_logit_approx(  -( nu[s,1,i] -  C_d[s,cutpoints[s,1]]));
                           bound_nd[1] = phi_logit_approx(  -( nu[s,2,i] -  C_nd[s,cutpoints[s,1]]));
                         
                           for (j in 2:(num_thresh_params[s])) {
                             bound_d[j] =  phi_logit_approx(  -( nu[s,1,i] -  C_d[s,cutpoints[s,j]] ) );
                             bound_nd[j] = phi_logit_approx(  -( nu[s,2,i] -  C_nd[s,cutpoints[s,j]]  )  );
                           }
                           for (j in 2:(num_thresh_params[s] + 1 )) {
                                  if (y[n,i,s] == j) {
                                    z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[j-1]  + (bound_d[j] - bound_d[j-1])* u_d1);      
                                    z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[j-1] + (bound_nd[j] - bound_nd[j-1])* u_nd1);    
                                    y1d[pa[n,s],i]   = log(bound_d[j] - bound_d[j-1]);  // Jacobian adjustment
                                    y1nd[pa[n,s],i]  = log(bound_nd[j] - bound_nd[j-1]); // Jacobian adjustment
                                  }
                                }
                                if (y[n,i,s] == 1 && i > num_binary_tests ) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[1]* u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[1]* u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d[1]); // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd[1]); // Jacobian adjustment
                                }
                                if (y[n,i,s] == 1 && i <= num_binary_tests ) {
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin + (1 - bound_d_bin)*u_d1);      
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin + (1 - bound_nd_bin)*u_nd1);    
                                  y1d[pa[n,s],i]   = log1m(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log1m(bound_nd_bin); // Jacobian adjustment
                                }
                               if (y[n,i,s] == 0) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin*u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin*u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd_bin); // Jacobian adjustment
                                }
                         }
                            lp1[pa[n,s]]  = sum(y1d[pa[n,s],])  +   bernoulli_lpmf(1 | p[s]); 
                            lp0[pa[n,s]]  = sum(y1nd[pa[n,s],]) +   bernoulli_lpmf(0 | p[s]);
                            log_lik_i[n,s] =  log_sum_exp(lp1[pa[n,s]] , lp0[pa[n,s]] ); 
                        }
                       }
                     }
}

model {
      p ~ beta(1,5); 

      b_primary_care ~ std_normal();

      C_dm  ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);
      C_ndm ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);

      z3 ~ std_normal();
      z4 ~ std_normal();
      sd3 ~ std_normal();

            sd1  ~   std_normal();
            L_Omega_bs1 ~ lkj_corr_cholesky(2);


       for (j in 1:2) 
           a1_m_raw[,j]  ~ std_normal(); 

           a2_m_raw  ~ std_normal(); 

        for (s in 1:NS) 
            z2[s, ] ~ std_normal();  
          

          for (s in 1:NS) { // likelihood
            for (n in 1:ns[s]) 
              target +=  log_lik_i[n,s]; 
          }
   
}

Is equivalent to the this model below (model # 2) which generates u_d and u_nd vectors for every individual

Model 2:



data {
          int num_binary_tests;
          int<lower=1> nt; // number of tests
          int<lower=1> NS;	// total # of studies
          int<lower=1> ns[NS]; // number of individuals in each study 
          int<lower=0> y[max(ns), nt, NS]; // N individuals and nt tests, NS studies 
          int num_patterns;
          int<lower=0> pa[max(ns), NS]; 
          int numg;
          int num_thresh_params[NS];
          int max_cutpoint;
          int cutpoints[NS, max_cutpoint];
          int twoway_thr[NS];
          int<lower=0, upper=1> primary_care_ind[NS]; //indicator for whether the study comes from primary care or not
          int r[NS, 4, choose(nt,2)];
          int num_refs; // number of reference tests
          int ref[NS]; // ref. test used in each of the NS studies
          int re[NS];
}

transformed data {
              int r2[choose(nt,2), NS, 4];
  
       for (a in 1:4)
          for (s in 1:NS)
            r2[1,s,a] = r[s,a,1];
}

parameters {
              ordered[2] a1_m_raw[num_refs]; // Se > Fp for each reference test
              ordered[2] a2_m_raw;   // don't order this due to potential thr. interction
              ordered[max_cutpoint] C_dm;
              ordered[max_cutpoint] C_ndm;
              vector<lower=0>[2] sd1;  
              real<lower=0,upper=1> u_d[num_patterns,nt]; // nuisance that absorbs inequality constraints
              real<lower=0,upper=1> u_nd[num_patterns,nt]; // nuisance that absorbs inequality constraints
              cholesky_factor_corr[2] L_Omega_bs1;
              vector[2] z2[NS];
              vector[NS] z3;
              vector[NS] z4;
              real<lower=0> sd3;
              real<lower=0, upper=1> p[NS]; 
              vector[2] b_primary_care; // coefficients for primary care (1 for MMSE Se and 1 for MMSE Sp)
}

transformed parameters {
              vector[2] mu[num_refs + 1];  
              matrix[2, nt] nu[NS];
              matrix[NS, num_patterns] log_lik_i; 
              vector[max_cutpoint] C_d[NS];
              vector[max_cutpoint] C_nd[NS];
              vector[num_patterns] lp1;
              vector[num_patterns] lp0;
              vector[nt] z_d[num_patterns];
              vector[nt] z_nd[num_patterns];
              vector[nt] y1d[num_patterns];
              vector[nt] y1nd[num_patterns];

                          // index test (MMSE)
                          mu[num_refs + 1, 1] =  a2_m_raw[1];
                          mu[num_refs + 1, 2] =  a2_m_raw[2];
                          
              for (s in 1:NS) {
                          // ref test (Se > Fp)
                          mu[ref[s], 1] =  a1_m_raw[ref[s], 2];   
                          mu[ref[s], 2] =  a1_m_raw[ref[s], 1];

                 // between-study model
                     // fixed ref. tests
                  nu[s, 1:2 , 1] =  mu[ref[s], 1:2]  ;  


                      // index test random effect - w/ covariate for primary care
                   nu[s, 1:2 , 2] =  mu[num_refs + 1 , 1:2 ]   + primary_care_ind[s]*b_primary_care[1:2] 
                                    + diag_pre_multiply(sd1, L_Omega_bs1) * z2[s,1:2] ; 

                for (i in 1:max_cutpoint) { 
                   C_d[s,i] = C_dm[i]  ; 
                   C_nd[s,i] = C_ndm[i]; 
                }
              }

 // Parameters for likelihood function. Based on code upoaded by Ben Goodrich which uses the 
                      // GHK algorithm for generating TruncMVN. See:
                      // https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan#L11
                      {
                 for (s in 1:NS) {
                    for (n in 1:ns[s]) {
                      for (i in 1:nt) { // loop over each test
                           real u_d1 = u_d[pa[n,s],i];
                           real u_nd1 = u_nd[pa[n,s],i];
                           real bound_d_bin     = phi_logit_approx(    - (nu[s,1,i])   ); // se
                           real bound_nd_bin    = phi_logit_approx(    - (nu[s,2,i])   );  // fp
                           vector[ num_thresh_params[s] + 1 ] bound_d;
                           vector[ num_thresh_params[s] + 1 ] bound_nd;
                           bound_d[ num_thresh_params[s] + 1 ] = 1;
                           bound_nd[ num_thresh_params[s] + 1 ] = 1;
                           bound_d[1] =  phi_logit_approx(  -( nu[s,1,i] -  C_d[s,cutpoints[s,1]]));
                           bound_nd[1] = phi_logit_approx(  -( nu[s,2,i] -  C_nd[s,cutpoints[s,1]]));
                         
                           for (j in 2:(num_thresh_params[s])) {
                             bound_d[j] =  phi_logit_approx(  -( nu[s,1,i] -  C_d[s,cutpoints[s,j]] ) );
                             bound_nd[j] = phi_logit_approx(  -( nu[s,2,i] -  C_nd[s,cutpoints[s,j]]  )  );
                           }
                           for (j in 2:(num_thresh_params[s] + 1 )) {
                                  if (y[n,i,s] == j) {
                                    z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[j-1]  + (bound_d[j] - bound_d[j-1])* u_d1);      
                                    z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[j-1] + (bound_nd[j] - bound_nd[j-1])* u_nd1);    
                                    y1d[pa[n,s],i]   = log(bound_d[j] - bound_d[j-1]);  // Jacobian adjustment
                                    y1nd[pa[n,s],i]  = log(bound_nd[j] - bound_nd[j-1]); // Jacobian adjustment
                                  }
                                }
                                if (y[n,i,s] == 1 && i > num_binary_tests ) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[1]* u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[1]* u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d[1]); // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd[1]); // Jacobian adjustment
                                }
                                if (y[n,i,s] == 1 && i <= num_binary_tests ) {
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin + (1 - bound_d_bin)*u_d1);      
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin + (1 - bound_nd_bin)*u_nd1);    
                                  y1d[pa[n,s],i]   = log1m(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log1m(bound_nd_bin); // Jacobian adjustment
                                }
                               if (y[n,i,s] == 0) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin*u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin*u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd_bin); // Jacobian adjustment
                                }
                         }
                            lp1[pa[n,s]]  = sum(y1d[pa[n,s],])  +   bernoulli_lpmf(1 | p[s]); 
                            lp0[pa[n,s]]  = sum(y1nd[pa[n,s],]) +   bernoulli_lpmf(0 | p[s]);
                            log_lik_i[s,pa[n,s]] =  log_sum_exp(lp1[pa[n,s]] , lp0[pa[n,s]] ); 
                        }
                       }
                     }
}

model {
      p ~ beta(1,5); 

      b_primary_care ~ std_normal();

      C_dm  ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);
      C_ndm ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);

      z3 ~ std_normal();
      z4 ~ std_normal();
      sd3 ~ std_normal();

            sd1  ~   std_normal();
            L_Omega_bs1 ~ lkj_corr_cholesky(2);


       for (j in 1:2) 
           a1_m_raw[,j]  ~ std_normal(); 

           a2_m_raw  ~ std_normal(); 

        for (s in 1:NS) 
            z2[s, ] ~ std_normal();  
          

          for (s in 1:NS) { // likelihood
            for (n in 1:ns[s]) 
              target +=  log_lik_i[s,pa[n,s]]; 
          }
}

Both of these work and give the same answers.

The problem is, when I allow for correlation, only the first version (which is less efficient) works. If i try to share nuisance parameters it doesn’t work - it runs quickly with 100% divergent transitions and the trace plots show that the parameters only explore a very narrow space, and some of them almost look like straight lines. So the model below (model # 3) does not work:

functions {
          // Function to do partial pooling on correlation matrices
          // See https://discourse.mc-stan.org/t/hierarchical-prior-for-partial-pooling-on-correlation-matrices/4852/27
          matrix convex_combine_cholesky(matrix global_chol_cor, matrix local_chol_cor, real alpha){
            int dim = rows(local_chol_cor);
            int global_dim = rows(global_chol_cor);
            matrix[global_dim,global_dim] global_cor = multiply_lower_tri_self_transpose(global_chol_cor);
            matrix[dim,dim] local_cor = multiply_lower_tri_self_transpose(local_chol_cor);
            matrix[dim,dim] global_cor_use;
            matrix[dim,dim] combined_chol_cor;
            
            if(dim < rows(global_cor)){
              global_cor_use = global_cor[1:dim,1:dim];
            } else {
              global_cor_use = global_cor;
            }
            combined_chol_cor = cholesky_decompose((1 - alpha)*global_cor_use + (alpha)*local_cor);
            return(combined_chol_cor);
  }
          real corr(int num, int[] vec1, int[] vec2) {
            real pearson_corr =      (num*sum(to_vector(vec1) .* to_vector(vec2)) -  sum(to_vector(vec1))*sum(to_vector(vec2)) )  /  
                              sqrt( (num*sum(to_vector(vec1) .* to_vector(vec1))   -  sum(to_vector(vec1))^2) * (num*sum(to_vector(vec2) .*  to_vector(vec2))   - sum(to_vector(vec2))^2 ) ); 
            return(pearson_corr);
          }
  real phi_logit_approx(real b) {
          real a;
          a = inv_logit( 1.7*b );
            return(a);
          }
  real inv_phi_logit_approx(real b) {
          real a;
          a = (1/1.7)*logit(b);    
            return(a);
          }

 real induced_dirichlet_lpdf(vector c, vector alpha, real phi) {
    int K = num_elements(c) + 1;
    vector[K - 1] sigma = inv_logit(phi - c);
    vector[K] p;
    matrix[K, K] J = rep_matrix(0, K, K);
    
    // Induced ordinal probabilities
    p[1] = 1 - sigma[1];
    for (k in 2:(K - 1))
      p[k] = sigma[k - 1] - sigma[k];
    p[K] = sigma[K - 1];
    
    // Baseline column of Jacobian
    for (k in 1:K) J[k, 1] = 1;
    
    // Diagonal entries of Jacobian
    for (k in 2:K) {
      real rho = sigma[k - 1] * (1 - sigma[k - 1]);
      J[k, k] = - rho;
      J[k - 1, k] = rho;
    }
return   dirichlet_lpdf(p | alpha)
           + log_determinant(J);
  }
}

data {
          int num_binary_tests;
          int<lower=1> nt; // number of tests
          int<lower=1> NS;	// total # of studies
          int<lower=1> ns[NS]; // number of individuals in each study 
          int<lower=0> y[max(ns), nt, NS]; // N individuals and nt tests, NS studies 
          int num_patterns;
          int<lower=0> pa[max(ns), NS]; 
          int numg;
          int num_thresh_params[NS];
          int max_cutpoint;
          int cutpoints[NS, max_cutpoint];
          int twoway_thr[NS];
          int<lower=0, upper=1> primary_care_ind[NS]; //indicator for whether the study comes from primary care or not
          int r[NS, 4, choose(nt,2)];
          int num_refs; // number of reference tests
          int ref[NS]; // ref. test used in each of the NS studies
          int re[NS];
          int total_n;
          int ind[NS];
          int ns_cumsum[NS];
}

transformed data {
              int r2[choose(nt,2), NS, 4];
  
       for (a in 1:4)
          for (s in 1:NS)
            r2[1,s,a] = r[s,a,1];
}

parameters {
              ordered[2] a1_m_raw[num_refs]; // Se > Fp for each reference test
              ordered[2] a2_m_raw;   
              ordered[max_cutpoint] C_dm;
              ordered[max_cutpoint] C_ndm;
              vector<lower=0>[2] sd1;  
              real<lower=0,upper=1> u_d[num_patterns,nt]; // nuisance that absorbs inequality constraints
              real<lower=0,upper=1> u_nd[num_patterns,nt]; // nuisance that absorbs inequality constraints
              cholesky_factor_corr[2] L_Omega_bs1;
              vector[2] z2[NS];
              real<lower=0, upper=1> p[NS]; 
              cholesky_factor_corr[nt] L_Omega_d[NS];
              cholesky_factor_corr[nt] L_Omega_nd[NS];
              vector[2] b_primary_care; // coefficients for primary care (1 for MMSE Se and 1 for MMSE Sp)    
}

transformed parameters {
              vector[2] mu[num_refs + 1];  
              matrix[2, nt] nu[NS]; 
              matrix[NS, num_patterns] log_lik_i; 
              vector[max_cutpoint] C_d[NS];
              vector[max_cutpoint] C_nd[NS];
              vector[num_patterns] lp1;
              vector[num_patterns] lp0;
              vector[nt] z_d[num_patterns];
              vector[nt] z_nd[num_patterns];
              vector[nt] y1d[num_patterns];
              vector[nt] y1nd[num_patterns];

                          // index test (MMSE)
                          mu[num_refs + 1, 1] =  a2_m_raw[1];
                          mu[num_refs + 1, 2] =  a2_m_raw[2];
                          
              for (s in 1:NS) {
                          // ref test (Se > Fp)
                          mu[ref[s], 1] =  a1_m_raw[ref[s], 2];   
                          mu[ref[s], 2] =  a1_m_raw[ref[s], 1];

                 // ref. tests
                   nu[s, 1:2 , 1] =  mu[ref[s], 1:2] ; 

                      // index test random effect - w/ covariate for primary care
                   nu[s, 1:2 , 2] =  mu[num_refs + 1 , 1:2 ]   + primary_care_ind[s]*b_primary_care[1:2] 
                                     + diag_pre_multiply(sd1*y1, L_Omega_bs1) * z2[s,1:2] ; 

                for (i in 1:max_cutpoint) { 
                   C_d[s,i] = C_dm[i];
                   C_nd[s,i] = C_ndm[i] ;
                }
              }

 // Parameters for likelihood function. Based on code upoaded by Ben Goodrich which uses the 
                      // GHK algorithm for generating TruncMVN. See:
                      // https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan#L11
                      {
                 for (s in 1:NS) {
                    for (n in 1:ns[s]) {
                      for (i in 1:nt) { // loop over each test
                           real u_d1 = u_d[pa[n,s],i];
                           real u_nd1 = u_nd[pa[n,s],i];
                           real prev_d = 0;
                           real prev_nd = 0;
                           real bound_d_bin     = phi_logit_approx(    - (nu[s,1,i] + prev_d)   / L_Omega_d[s,i,i] ); // se
                           real bound_nd_bin    = phi_logit_approx(    - (nu[s,2,i] + prev_nd ) / L_Omega_nd[s,i,i] );  // fp
                           vector[ num_thresh_params[s] + 1 ] bound_d;
                           vector[ num_thresh_params[s] + 1 ] bound_nd;
                           bound_d[ num_thresh_params[s] + 1 ] = 1;
                           bound_nd[ num_thresh_params[s] + 1 ] = 1;

                           for (j in 1:(num_thresh_params[s])) {
                             bound_d[j] =  phi_logit_approx(  -( nu[s,1,i] -  C_d[s,cutpoints[s,j]] + prev_d ) / L_Omega_d[s,i,i] );
                             bound_nd[j] = phi_logit_approx(  -( nu[s,2,i] -  C_nd[s,cutpoints[s,j]] + prev_nd) /L_Omega_nd[s,i,i] );
                           }

                           for (j in 2:(num_thresh_params[s] + 1 )) {
                                  if (y[n,i,s] == j) {
                                    z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[j-1]  + (bound_d[j] - bound_d[j-1])* u_d1);      
                                    z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[j-1] + (bound_nd[j] - bound_nd[j-1])* u_nd1);    
                                    y1d[pa[n,s],i]   = log(bound_d[j] - bound_d[j-1]);  // Jacobian adjustment
                                    y1nd[pa[n,s],i]  = log(bound_nd[j] - bound_nd[j-1]); // Jacobian adjustment
                                  }
                                }
                                if (y[n,i,s] == 1 && i > num_binary_tests ) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d[1]* u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd[1]* u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d[1]); // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd[1]); // Jacobian adjustment
                                }
                                if (y[n,i,s] == 1 && i <= num_binary_tests ) {
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin + (1 - bound_d_bin)*u_d1);      
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin + (1 - bound_nd_bin)*u_nd1);    
                                  y1d[pa[n,s],i]   = log1m(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log1m(bound_nd_bin); // Jacobian adjustment
                                }
                               if (y[n,i,s] == 0) { 
                                  z_d[pa[n,s],i]   = inv_phi_logit_approx(bound_d_bin*u_d1);
                                  z_nd[pa[n,s],i]  = inv_phi_logit_approx(bound_nd_bin*u_nd1);
                                  y1d[pa[n,s],i]   = log(bound_d_bin);  // Jacobian adjustment
                                  y1nd[pa[n,s],i]  = log(bound_nd_bin); // Jacobian adjustment
                                }
                              if (i < nt) prev_d  = L_Omega_d[s,i +1,1:i ] * head(z_d[pa[n,s],] ,i);
                              if (i < nt) prev_nd = L_Omega_nd[s,i +1,1:i ] * head(z_nd[pa[n,s],] ,i);
                         }
                            lp1[pa[n,s]]  = sum(y1d[pa[n,s],])  +   bernoulli_lpmf(1 | p[s]); 
                            lp0[pa[n,s]]  = sum(y1nd[pa[n,s],]) +   bernoulli_lpmf(0 | p[s]);
                            log_lik_i[s,pa[n,s]] =  log_sum_exp(lp1[pa[n,s]] , lp0[pa[n,s]] ); 
                        }
                        }
                       }
                     }
}

model {

      p ~ beta(a1,b1); 

      b_primary_care ~ std_normal();

      C_dm  ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);
      C_ndm ~ induced_dirichlet(rep_vector(1, max_cutpoint + 1), 0);

            sd1  ~   std_normal();
            L_Omega_bs1 ~ lkj_corr_cholesky(2);


       for (j in 1:2) 
           a1_m_raw[,j]  ~ std_normal(); 

           a2_m_raw  ~ std_normal(); 


         for (s in 1:NS) {
          L_Omega_d[s,] ~ lkj_corr_cholesky(x1);
          L_Omega_nd[s,] ~ lkj_corr_cholesky(x2);
         }

        for (s in 1:NS) 
            z2[s, ] ~ std_normal();  
          
           for (s in 1:NS) { // likelihood
            for (n in 1:ns[s]) 
              target +=    log_lik_i[s,pa[n,s]] 
          } 
}

I was wondering whether anybody would have any ideas as to why would the model with no correlation and shared nuisance parameters (model #2) works, but not the same model with correlation (model # 3)?

1 Like

Hi, your question is relevant, but it is honestly very hard to answer: the models you posted are huge and I am unable to understand what’s going on in them without significant time investment. I think you are more likely to get an answer if you find a way to strip the models down to the smallest case that shows the issue (often this process in preparing a question helps me to find the answer myself :-)

Best of luck with the model!

1 Like

Thanks, I’ll work on a simpler version and update