Using the onion method in occupancy hierarchical model

Modeling
,
2

I don’t have time to understand your exact model but I’ll give you how to use the onion method. The linked thread just shows how to reconstruct the cov mat in gen quant but you can move that to model or tp such as (nb, haven’t tested this code so there may be errors)

data {
  int<lower = 0> K;    // dim of corr mat
  real<lower = 0> eta;  // lkj param
}
transformed data {
  real<lower = 0> alpha = eta + (K - 2) / 2.0;
  vector<lower = 0>[K-1] shape1;
  vector<lower = 0>[K-1] shape2;

  shape1[1] = alpha;
  shape2[1] = alpha;
  for(k in 2:(K-1)) {
    alpha -= 0.5;
    shape1[k] = k / 2.0;
    shape2[k] = alpha;
  }
}
parameters {
  row_vector[choose(K, 2) - 1]  l;         // do NOT init with 0 for all elements
  vector<lower = 0,upper = 1>[K-1] R2; // first element is not really a R^2 but is on (0,1)  
}
model {
  // priors that imply Sigma ~ lkj_corr(eta)
  l ~ std_normal();
  R2 ~ beta(shape1, shape2);

 matrix[K, K] L = rep_matrix(0, K, K); // cholesky_factor corr matrix
  {
    int start = 1;
    int end = 2;

    L[1,1] = 1.0;
    L[2,1] = 2.0 * R2[1] - 1.0;
    L[2,2] = sqrt(1.0 - square(L[2,1]));
    for(k in 2:(K-1)) {
      int kp1 = k + 1;
      row_vector[k] l_row = segment(l, start, k);
      real scale = sqrt(R2[k] / dot_self(l_row));
      L[kp1, 1:k] = l_row[1:k] * scale;
      L[kp1, kp1] = sqrt(1.0 - R2[k]);
      start = end + 1;
      end = start + k - 1;

    }
}
generated quantities {
  matrix[K,K] Sigma = multiply_lower_tri_self_transpose(L);
}
4 Likes