return to lmer vs stan

[Douglas Bates]

I am dipping my toes into stan, coming from a background in maximum likelihood estimation of linear mixed-effects models (I am one of the authors of the lme4 package for R and also of the MixedModels package for Julia). I also have some familiarity with C++, Rcpp and Eigen having done the interface between R/Rcpp and Eigen for the RcppEigen package.

I have at least skimmed the Reference Manual for Stan 2.2.0 and read some of the sections pretty thoroughly.

I very much like the idea of the 

cholesky_factor_cov

and

cholesky_factor_cor

classes.  As the Stan authors may know, lmer uses the elements of the cholesky_factor_cov representation of the (relative) covariance matrices for the random effects as the parameters being directly optimized in the maximum likelihood estimation of parameters.  (The fixed-effects coefficients and the residual standard deviation are profiled out of the log-likelihood function.)

My first question is regarding a prior on the parameters in a cholesky_factor_cor object.  Sections 46 and 47 of the manual discuss distributions of covariance matrices and correlation matrices.  Is it possible to assign a prior distribution to the Cholesky factor directly? For example, in the lkj_corr_log method defined in lkj_corr.hpp the correlation matrix is factored before being passed to lkj_corr_choleksy, but lkj_corr_cholesky_log.  Have I missed a definition of lkj_corr_cholesky or lkj_cov_cholesky to be used in an expression like

L~lkj_corr_cholesky(eta)

or do I go the direct route and call lkj_corr_cholesky_log to increment the log probability?

[Bob Carpenter]

On Mar 17, 2014, at 6:32 PM, Douglas Bates <dmb...@gmail.com> wrote:

> I am dipping my toes into stan,

Welcome!

> coming from a background in maximum likelihood estimation of linear mixed-effects models (I am one of the authors of the lme4 package for R and also of the MixedModels package for Julia).

We've all heard a lot about lme4 from Andrew (who calls it "lmer").
We've been spending some time going the other way and thinking
about (penalized) MLE and (penalized) maximum marginal likelihood for
Stan.

> I also have some familiarity with C++, Rcpp and Eigen having done the interface between R/Rcpp and Eigen for the RcppEigen package.
>
> I have at least skimmed the Reference Manual for Stan 2.2.0 and read some of the sections pretty thoroughly.
>
> I very much like the idea of the
>
> cholesky_factor_cov
>
> and
>
> cholesky_factor_cor
>
> classes. As the Stan authors may know, lmer uses the elements of the cholesky_factor_cov representation of the (relative) covariance matrices for the random effects as the parameters being directly optimized in the maximum likelihood estimation of parameters. (The fixed-effects coefficients and the residual standard deviation are profiled out of the log-likelihood function.)

Only a few of us are versed in the lme4 code.

> My first question is regarding a prior on the parameters in a cholesky_factor_cor object. Sections 46 and 47 of the manual discuss distributions of covariance matrices and correlation matrices. Is it possible to assign a prior distribution to the Cholesky factor directly? For example, in the lkj_corr_log method defined in lkj_corr.hpp the correlation matrix is factored before being passed to lkj_corr_choleksy, but lkj_corr_cholesky_log. Have I missed a definition of lkj_corr_cholesky or lkj_cov_cholesky to be used in an expression like
>
> L~lkj_corr_cholesky(eta)
>
> or do I go the direct route and call lkj_corr_cholesky_log to increment the log probability?

That function no longer exists as of Stan 2.2 (though it's still there in C++ if
that's what you're after). We're going to bring it back when we've added the cholesky_corr data type,
which is nearly ready to go. Without the data type, the density's useless because you'd have to
do the decomposition and Jacobian adjustment yourself in the Stan code, which
is pretty awful.

For now (Stan 2.2), you have to go the roundabout way of defining
a correlation matrix parameter:

parameters {
...
corr_matrix[K] Omega;
...
}

and then give it an LKJ prior:

model {
matrix[K,K] L; // local variable

Omega ~ lkj_corr(nu); // prior directly on correlation matrix
L <- cholesky_decompose(Omega);

... use L ...
}

This introduces a lot of redundant calculations because the underlying representation
for Omega in the unconstrained space is as its Cholesky factor, so we have to
do a redundant multiplication and Jacobian adjustment to get to the correlation
matrix itself in the constrained space, and then I believe lkj_corr does the Cholesky
factorization under the hood, and then cholesky_decompose does the factorization again.

All the redundant calculation will go away when we have the cholesky_corr data type, and
you'll be able to write something like this (we'll probably simplify the name of the
density):

parameters {
cholesky_corr L;
...
model {
L ~ lkj_cholesky(nu);
... use L ...

Ben Goodrich has a neat paper comparing Wishart and LKJ-style priors (he was
the one driving the design of these parts of the code):

http://www.stat.columbia.edu/~gelman/research/unpublished/Visualization.pdf

We're also adding a quad_form_diag(Omega,sigma) function in the next release
of Stan that will do

diag_matrix(sigma) * Omega * diag_matrix(sigma)

to make the kind of scaled covariance priors Andrew and Ben like easier to express.

- Bob

[Douglas Bates]

I should have mentioned that I am working on the master branch from github which is why I was writing about the cholesky_factor_corr class.

I had forgotten myself that I was working with a pre-release version of stan.  I tend to work off the github repository so that I can read the code if needed.

[Daniel Lee]

The master branch is actually our latest tagged version, so you're looking at v2.2.0. The develop branch is our current working version.

Daniel

On Mon, Mar 17, 2014 at 2:55 PM, Douglas Bates <dmb...@gmail.com> wrote:

I should have mentioned that I am working on the master branch from github which is why I was writing about the cholesky_factor_corr class.

I had forgotten myself that I was working with a pre-release version of stan.  I tend to work off the github repository so that I can read the code if needed.

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

The master branch is just the latest release. The
develop branch is where to go for the latest code
for the next release. There is also a ton of code
flying around in branches --- we don't merge into
develop until we have unit tests passing and end-to-end
integration tests passing. So develop should, theoretically,
always be in a consistent, "ready to release" state.

- Bob

On Mar 17, 2014, at 7:55 PM, Douglas Bates <dmb...@gmail.com> wrote:

> I should have mentioned that I am working on the master branch from github which is why I was writing about the cholesky_factor_corr class.
>
> I had forgotten myself that I was working with a pre-release version of stan. I tend to work off the github repository so that I can read the code if needed.
>

[Ben Goodrich]

Hi Doug,

Good to hear from you again. As Bob mentioned, the lkj_corr_cholesky() function is easier to use incorrectly than to use correctly, which is why it is slated to get disappeared for a while. In short, I intended it to be an internal-only C++ function but it accidentally got exposed to the Stan language in the initial release of Stan. The only reason to use it that is still valid (i.e. it corresponds to the posterior distribution you defined) is if a correlation matrix is a parameter and you need its Cholesky factor for some subsequent calculation. In that case,

model {
  L <- cholesky_decompose(correlation_matrix);
  L ~ lkj_choleksy_corr(eta); // this confusingly implies correlation_matrix ~ lkj_corr(eta)
  // use L in a subsequent calculation to avoid recalculating it
}

In all other cases in the Stan language,

foo ~ some_distribution();

means that you are putting a prior on foo itself.

What, you, I, and everybody else really needs is for the Cholesky factor of a correlation matrix to be a primitive type in Stan. This is possible with some convoluted hand-written Stan code that I am attaching. It constructs L via something equivalent to the onion method described in the LKJ (2008) paper such that

L * L' ~ lkj_corr(eta);

without actually forming that matrix product. To verify (a necessary condition that) L * L' has the LKJ distribution, you can do in R

library(rstan)
K <- 11
posterior <- stan(file = "onion.stan", data = list(K = K, eta = 1))
correlation_matrices <- extract(posterior, "Sigma")[[1]]
rho <- correlation_matrices[,1,2]
qqplot(rho, -1 + 2 * rbeta(length(rho), K / 2, K / 2), pch = 20)
abline(a = 0, b = 1, col = "red", lty = "dashed")

Hope that helps,
Ben