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)?