Example for reparametrisation of multivariate normal parameters in a mixed effects model

[Bob_Carpenter]

@mesudell

in my data example, M>300,000, N= roughly 30000, K=5,

This is a lot of data for Stan.

@mesudell

this is an exercise to allow me to understand what I need to do with a more complex model

Great! We always recommend people start as simply as possible and build up. It’s how we code our own models.

There are a bunch of things you can do to make the code a bit faster, but nothing that’s going to change the order by much. For example, you define Dl and Al but never use them. If you really really want them in the output but aren’t going to use them in the model, then move the definition to generated quantities and you won’t have all the autodiff overhead for unused variables. You can also inline definitions, so it’s

matrix[N, ql] b2_L = (diag_pre_multiply(tau_Dl, Omega_Dl) * z_b2_L)';

(I switched to the operator transposition.). Given that you’re transposing, which is expensive, you probably want to redeclare these things as

(diag_pre_multiply(A, B) * Z)'
    = Z' * diag_pre_multiply(A, B)'
    = Z' * diag_post_multiply(B', A')

So if you transpose all the variables Z, B, A in their declarations, this gets simpler. It may show up elsewhere—I didn’t try to track through other uses of these variables that may have other requirements.

@mesudell

Am I correct I don’t need one here (mathematical question apologies)

Yes—you don’t apply non-constant non-linear functions to variables and then put them on the left-hand side of sampling statements.

You only need a Jacobian adjustment when you apply a function f() to a variable and then sample that. Here’s the simplest case that needs a Jacobian adjustment, where I’m assuming a function f with derivative function df_dx.

  real y = f(x);
  target += log(abs(df_dx(x)));
  ...
  y ~ foo(theta);  // Jacobian needed

On the other hand, if you use it on the right you don’t need to adjust,

real theta = f(phi);
...
y ~ bar(theta);  // no Jacobian needed

@mesudell

matrix matrixlongver2(row_vector X, int m){
    matrix[m,cols(X)] Xout;
    for(i in 1:m){
      Xout[i,] = X;
    }
    return Xout;
  }

Things like this cause a lot of memory pressure with duplicated entires. Are you use you need a matrix here? Also, if you do need it, it can be replaced with the following (it will duplicate rows given row vector or duplicate cols given column vector).

matrix[m, cols(X)] Xout = rep_matrix(X, m);

Sorry that’s not helping with the reparameterization. I wasn’t clear on what you thought needed to be reparameterized here. With the amount of data you have, using the non-centered parameterizations may be worse than using centered ones. It comes down to trial and error, I’m afraid.