Tree-structured covariance in Stan

Modeling
1

Hi all,

I’d like to impose a tree structure on a covariance matrix I’m estimating from some data. The quantity of interest here is the rooted tree (topology and vector of branch lengths). I’m working off of Bravo et al 2009 (link to PDF here), which specifies two recipes for imposing a tree structure on a covariance matrix, but I don’t know how to implement these recipes in Stan, as I think they’d involve dynamically declaring parameter bounds. Any tips on whether I can and/or how to implement either of these recipes in Stan appreciated.

Terminology:
B is the tree-structured covariance, of dimension p x p

Recipe 1:
B = VDV^T, where:

  • D is a diagonal matrix of size p x p that has nonnegative entries
  • V is a binary basis matrix of size p x 2p-1 with unique columns that satisfies the following properties:
    • (1) V contains the vector of all 1s as a column
    • (2) for every column w in V with more than one non-zero entry, it contains exactly 2 columns u and v such that u+v=w

Recipe 2:
The entries of B satisfy the following constraints:
B_{ij} \geq 0 \: \forall{i,j}
B_{ii} \geq B_{ij} \: \forall i \neq j
B_{ij} \geq B_{ik} - (1 - \rho_{ijk1}M)
B_{ik} \geq B_{jk} - (1 - \rho_{ijk1}M)
B_{jk} \geq B_{ik} - (1 - \rho_{ijk1}M)

B_{ik} \geq B_{ij} - (1 - \rho_{ijk2}M)
B_{ij} \geq B_{jk} - (1 - \rho_{ijk2}M)
B_{jk} \geq B_{ij} - (1 - \rho_{ijk2}M)

B_{jk} \geq B_{ij} - (\rho_{ijk1} - \rho_{ijk2}M)
B_{ij} \geq B_{ik} - (\rho_{ijk1} - \rho_{ijk2}M)
B_{ik} \geq B_{ij} - (\rho_{ijk1} - \rho_{ijk2}M)

\rho_{ijk1} + \rho_{ijk2} \leq 1
\rho_{ijk1},\rho_{ijk2} \in \{0,1\} \: \forall i > j > k

where M is a very large positive constant and \rho_{ijk1} and \rho_{ijk2} are integer variables

2

This is what prevents you from being able to declare V in the parameters block of a Stan program because the posterior is not differentiable with respect to V. In principle, you could enumerate every possible V in the transformed data block and marginalize over them, but that is a huge number of matrices if p is not very small.

3

Sorry @bgoodri fingers slipped and posted the original post prematurely. Edited to include Recipe 2. But also thanks for the fast response!!!

4

What is M?

5

Sorry - M is a very large positive constant, and \rho_{ijk1} and \rho_{ijk2} are integer variables. [edited original post to include]

6

If those are unknown integer variables, then you have the same problem with lack of differentiability as in recipe 1, although the marginalization route may be more feasible than before.