New User's Guide Intro, SEM, and user-defined constrained parameters chapters

General
1

After feedback from @WardBrian and @mitzimorris, I just merged two new chapters for the User’s Guide right up front in a new introductory part.

  1. The Stan Ecosystem
  2. How Stan Works

I’m not sure what the schedule is for releasing doc after we merge, but I think it should be live before our next release.

Chapter (1) tries to explain where everything is, so it’s a rehash of a lot of what’s on our web pages, but spelled out linearly with a bit more explanation.

Chapter (2) is a short overview of how the language works, including variables, blocks, and transforms.

I’m nearly done writing two more chapters:

(n) Structure equation models
(n + 1) Custom constrained parameters

Hopefully (n) is self explanatory. I’m just going over notation and showing how to code in Stan explicitly and with soft zero constraints. I met Ken Bollen (author of the main textbook in the area) at the Modern Modeling Methods conference and have been getting clarification as to whether I correctly understand the notation (that shrinking non-linked correlations to zero still feels wrong but I don’t plan to editorialize in the chapter). Of course, I mention blavaan as the thing to use if your model fits it.

For (n + 1), I’m writing about a bunch of different constraints that we don’t talk about elsehwere.

A. Orthonormal, orthogonal, and rotation matrices
B. Matrices matching specified margins with sum-to-zero matrix (thanks, @spinkney)
C. Fun shapes (from BUGS) like discs, hollow squares, rings, paralellograms, hexagons, etc.

I plan to add Sarah Heaps’s stationarity constraints for vector autoregressive (VAR) models to (n + 1). Then once these are in the User’s Guide, we can move some of them into the language itself. The User’s Guide is like our Boost :-).

I’m open to suggestions on what else to include in both new chapters. I’m particularly looking for constraints that might be widely used or illustrate a general technique.

2

You can parameterize exact 0s in the correlation matrix. I showed this as the 2024 stancon. I can help write that if you need.

On tables with given row and column sums, an additional constraint is that the interior cells be non-negative. A sum to zero matrix doesn’t enforce that last constraint. If the marginal sums are sufficiently large then the largest negative deviation from the sum to zero matrix is typically small enough to keep the cell contained to be non-negative. Otherwise, there’s a couple ways to parameterize this additional constraint.

How much more do you want to add? There’s linear equalities, linear inequalities (you can solve Ax < b or do it telescopically if you have parameters that depend on others), and ordered matrix.

3

I mean that, sometimes one wants an additional positivity constraint.

4

Thanks for the reminder, @spinkney.

I found the video:

and the matching paper:

https://arxiv.org/pdf/2405.07286

After equation (2), it allows range based bounds, but explicitly excludes the case where the lower bound (a) is equal to the upper bound (b), -1 \leq a_{i,j} < C_{i, j} < b_{i,j} \leq 1. The paper says it handles “known zeroes”, but I don’t see where that’s done. Wouldn’t it have to reduce the dimensionality of the inputs below n \choose 2, for example, if all the off-diagonal values were constrained to zero.

5

Yes, that’s correct. I have R and Stan code for this in open-talks/structured-corr-stancon24/R/code.R at main · spinkney/open-talks · GitHub

6

Related to (n), I recently wrote a review of popular Bayesian approaches to SEM estimation. Some parts will not be helpful to you, but Section 4.1 has a bit of discussion about notation and choice of likelihood.

7

Thanks for the link. It’s been fun coming to grips with this notation (I am first and foremost, a formal systems semantics geek and only secondarily a statistician and linear algebraist).

The author of the reply and paper, @edm, is also the author of the blavaan package for Bayesian SEM, so this should be a definitive treament. It covers Gibbs (for historical reasons?) and HMC in a reasonable way.

The good news is that I’ve sorted @spinkney’s structured correlation matrix construction so I can correctly code SEM models in Stan for the User’s Guide chapter and give them LKJ priors with coupled scale priors. I’m first adding a new chapter on user-defined constrained types which will go through @spinkney’s construction, which I can then reference when coding up SEM.

P.S. For anyone considering inverse Wishart priors for covariances, I’d recommend the paper “Visualizing distributions of covariance matrices” from @bgoodri and others.

8

This will do more for getting people to use the method than I could ever do. Thanks @Bob_Carpenter for tirelessly supporting Stan documentation and resources!

9

Here’s a draft of the explicit zeros code, which is made up of two functions. I made it more general and less efficient by treating all the rows with the general case. I’m still not sure about representing the stick on the log scale—it’s probably overkill.

But I must have introduced a bug, because it just locks up under sampling :-(.

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }
  
  matrix cholesky_corr_zeros_jacobian(int D,
                                      vector raw,
                                      array[,] int zeros) {
    int raw_idx = 1;
    int zero_idx = 1;
    matrix[D, D] L = rep_matrix(0, D, D);
    for (i in 1:D) {
      real log_sq_stick = 0;
      for (j in 1:(i - 1)) {
        if (is_zero(i, j, zeros, zero_idx)) {
          real b = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
          L[i, j] = -b / L[j, j];  // implies Omega[i, j] == 0
          zero_idx += 1;
        } else {
          real stick = exp(0.5 * log_sq_stick);
          // print("stick = ", stick);
          L[i, j] = lower_upper_bound_jacobian(raw[raw_idx], -stick, stick);  // inside stick
          raw_idx += 1;
        }
        log_sq_stick = log_diff_exp(log_sq_stick, 2 * log(abs(L[i, j])));
      }
      L[i, i] = exp(0.5 * log_sq_stick);
    }
    return L;
  }

  real lkj_cholesky_corr_zeros_lpdf(matrix L,
                                    real nu,
                                    array[,] int zeros) {
    real lp = lkj_corr_cholesky_lpdf(L | nu);  // over-adjusts
    int N_zero = size(zeros);
    for (n in 1:N_zero) {                      
      int col_idx = zeros[n, 2];
      lp -= log(L[col_idx, col_idx]);          // correct over-adjustment
    }
    return lp;
  }
}
data {
  int<lower=2> D;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(D, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(D, 2) - N_zero] raw;           // raw parameters
}
transformed parameters {
  matrix[D, D] L_Omega = cholesky_corr_zeros_jacobian(D, raw, zeros);
}
model {
  L_Omega ~ lkj_cholesky_corr_zeros(eta, zeros);
}
generated quantities {
  matrix[D, D] Omega = multiply_lower_tri_self_transpose(L_Omega);
}
10

Man, I think I once implemented linear equalities on my own, but I can’t recall ever figuring out how to do linear inequalities in any kind of general way. The stuff Bob mentioned and the stuff you mentioned are some of the least straightforward to do in Stan. And even though I don’t write as much Stan code as I used to, it will be much appreciated to see documentation and/or recommendations on best practices on them.

11

I don’t find the Stan part of it hard, but the actual transforms are very involved, especially if you aren’t already familiar with the component pieces like building Cholesky factors for correlation matrices. I remember my chagrin when @bgoodri gave me the Lewandowski, Kurowicka, and Joe paper (that led to the LKJ prior Ben developed)—I couldn’t understand any of it then.

Anyway, I’m going to finish this up and release this week now that the Walnutpie release is out.

12

I agree with and relate to every word of this.