Exploratory Factor Analysis and the Cholesky Factor

[Rick Farouni]

Hi,

Hi how can I constrain  the matrix" L " (initialized as  cholesky_factor_cov[M,K]  L where M>K ) to have a diagonal of ones?  I am referring to the standard  "hierarchical" structural constraint for the loadings matrix used in Bayesian factor analysis (e..g. Aguilar and West, 2000). 

Aguilar, Omar, and Mike West. "Bayesian dynamic factor models and portfolio allocation." Journal of Business & Economic Statistics 18.3 (2000): 338-357.

[Christian Roy]

Hi Rick,

I don't think my solution will help you out but I managed to put constraints on a matrix by something like this. 

{

int counter;

for(i in 1:N_row){

  alpha_matrix[i,i] <- 1;

  for(j in (i+1):N_col){alpha_matrix[i,j] <- 0};

}

for(i in 2:N_row){

  for(j in 1:(i-1)){

  counter  <- counter+1;

  alpha_matrix[i,j] <- something[counter];

  }}

}

Obviously the problem is with this hack is the lower part of the matrix where you have to include an index to account for the "something" part. It will get uglier if there are constraints on the matrix as we would expect with Cholesky factors covariance matrices.

Hopefully someone will have a more intelligent solution

Chris

-------------- 

--
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.

[Michael Betancourt]

You should be able to do this by defining not using a cholesky_factor_cov as a parameter directly.

Instead define unconstrained variables for the off-diagonal elements and then construct a 

cholesky_factor_cov in the transformed parameters block, filled with ones on the diagonal and

the unconstrained variables off the diagonal.

[Bob Carpenter]

I'll try to add a working example the manual---we need some
more realistic, but complex constrained parameter examples.

You want to do what Michael Betancourt suggested. You can then use
it as the covariance parameter with multi_normal_cholesky().

If you transform it further to a covariance matrix and try to put
a prior on the covariance matrix, you'll need to include the Jacobian
adjustment as described in the changes of variables chapter in the
manual.

- Bob

[Rick Farouni]

Thanks everyone,

I think I managed to get the constraints working  for an KxD rectangular factor loading matrix, where D is the number of dimensions (i.e. factors).  Here is the code:

transformed parameters{

matrix[K,D]  alpha;

{

int index;

for (j in 1:D) {

    alpha[j,j] <- 1;

    for (i in (j+1):K) {

      index <- index + 1;

      alpha[i,j] <- offdiag[index];} // offdiag is a D*(K-D)+(D*(D-1)/2) vector specified a prior in the model block

}

for(i in 1:K){

     for(j in (i+1):D){

     alpha[i,j] <- 0;}

}

}

}

However, the constraints themselves don't seem to be sufficient to restrict the likelihood to a single mode, which somehow contradicts the literature. Also, why are the enclosing braces necessary to define" int index"? Without them, I get an error message "integer parameters or transformed parameters are not allowed.

Rick,

[terrance savitsky]

​

 befa-aoas-r5.pdf

​hi, for whatever it's worth, i attach a paper under review that performs bayesian efa estimation without any restriction on the structure of the loadings matrix to achieve a marginally exchangeable prior for the communality.  the benefit of such exchangeability is that it removes the dependence on variable ordering associated to the lower triangular restriction.  rotational identifiability of the draws is established in a simple post-processing step.   the whole approach was designed for implementation under black box solvers like STAN.

terrance​

 supplement-r5.pdf

​

Thank you, Terrance Savitsky

[Rick Farouni]

Thanks! It's worth a lot!!  I was also looking into post processing approaches, for example http://www.econstor.eu/bitstream/10419/67349/1/731510127.pdf and https://www.stat.washington.edu/research/reports/2011/tr589.pdf. Is your approach different? Anyways, I will try to implement the method described in your paper in STAN  and I will get back to you.

[Michael Betancourt]

We don’t allow integer parameters or transformed parameters because they do not 

worth with Hamiltonian Monte Carlo.  In your code the brace opens up a local

scope such that index is just a local variable used in the computation of alpha,

but not a transformed parameter in of itself.

[terrance savitsky]

yes,  you'll just use an optimized orthogonal rotation criterion - like varimax - to rotate the draws and then resolve permutations in the manner we show in the paper.   the whole thing is much simpler than the alternatives i'd initially considered.  the real reason to dispense with the loadings restrictions is because it removes dependence on ordering of variables (unlike for MLE), rather than implementation in STAN, but you'll achieve better mixing - particularly if you use parameter expansion - and once you set up the two post-processing steps, you'll be able to easily do this in any future models you'd conceive.

[Ben Goodrich]

On Saturday, June 21, 2014 12:25:58 PM UTC-4, Rick Farouni wrote:

However, the constraints themselves don't seem to be sufficient to restrict the likelihood to a single mode, which somehow contradicts the literature.

The reason there is not a single mode is that you can permute the factors and permute the loadings accordingly without changing the likelihood. Also, post-processing the output of Stan tends not to work well when there are multiple equivalent modes. You might be better off restricting the diagonal of the loadings matrix to be a copy of a positive_ordered variable and restrict the factor variances to be one.

Ben

[terrance savitsky]

the two steps the paper i attached proposes to post-identify draws will work and the steps are easy to implement, so it's not hard to test.  the varimax step reduces general rotational invariance (which is uncountable) to a countable permutation invariance (with sign flips).  this is just a standard label switching problem and the second step in the paper resolves it.   this procedure will work with output from STAN, JAGS or any other solver.  the point is not the estimation algorithm (assuming the chains are converged), but the probability model definition.

mixing issues aside, the fundamental problem with using a lower triangular with positive diagonals restriction on the loadings matrix to fix a rotation is that this constraint induces a non-exchangeable prior on the communality - where founder factors receive a smaller prior variance than other loadings parameters.  the result (which is well known) is that the posterior is sensitive - varies with the ordering of variables, unlike EFA under MLE.  so it's really a good idea to dispense with the restriction. 

--
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.

--
Thank you, Terrance Savitsky

[Ben Goodrich]

But the initial problem is that NUTS sampler has a hard time adapting and moving through the space when different parameter proposals yield the same value of the log-posterior function, such as when you can rotate, permute, rescale, etc. without changing the model. So, post-processing schemes are post-processing a dubious set of initial draws.

Ben

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[terrance savitsky]

i see. thanks for clarifying.  i hadn't thought that NUTS wouldn't have problems exploring multiple equivalent modes, in part because HMC produces better mixing.  is that the gradient descent proposal step directions struggle under multiple equivalent modes?  

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.

[Ben Goodrich]

For multiple modes --- regardless of whether they are exactly equivalent --- HMC will have difficulty because the Hamiltonian proposal tends to be toward the closest mode, but that is true for a lot of MCMC algorithms. The adaptation of NUTS is further jeopardized for a variety of reasons including that it is assuming constant curvature, which wouldn't be the case under multimodality. Furthermore, if the posterior distribution has a long ridge, it might keep leapfrogging indefinitely along that ridge. Basically, it is not real good at all in situations like this. There is more somewhere in the manual in the subsection about difficult posteriors.

Ben

[Rick Farouni]

Hi Ben,

I have attached a STAN file for a simple Factor Analysis model  and an R code that generated the data. I have restricted the diagonal to positive _ordered, but it seems it made things worse. How can I modify the code so I can get a meaningful result in this case?  In the model I used the data (full information factor analysis), not the correlation matrix. Also, I read the  NxK data matrix using "vector[K]  y[N];" is that correct?

Thanks,

Rick

[Ben Goodrich]

I think the main problem here is that it is inappropriate to use an iid standard normal prior on a positive_ordered variable. That said, we don't have any good suggestions for such a prior. So, probably in this case, it would be better to make the diagonal elements positive but not necessarily ordered and put some appropriate prior on them.

You can speed up your code a little bit by making F a matrix and y a concatenated vector, in which case the guts of the model could be written as

  offdiag ~ normal(0,sigma_L);
  diag ~ something;    
  sigma_L ~ cauchy(0,5);
  to_vector(F) ~ normal(0,1);
  y ~ normal(to_vector(F * L'), 1);

That still converges slowly with all the factor scores. So, for a pure EFA model the way psychologists estimate it, I would specify it in terms of the manifest covariance matrix having a Wishart likelihood.

Ben

[Bob Carpenter]

On Jun 23, 2014, at 6:13 PM, Ben Goodrich <goodri...@gmail.com> wrote:

> I think the main problem here is that it is inappropriate to use an iid standard normal prior on a positive_ordered variable.

Inappropriate how?

I just made a note to drop in some advice in the next
version of the manual, but I'm not sure what to tell people.

If you give an ordered vector variable components an i.i.d. normal,
then I believe the result is the same as if you'd provided
the same prior for an unconstrained vector and then sorted.

Another alternative would be to put a prior on the first
element and then on the differences (probably two different
priors, and the prior is probably different for positive_ordered
and ordered).

I remember asking Andrew what he'd prefer here without getting
an answer.

- Bob

[Ben Goodrich]

On Monday, June 23, 2014 1:05:40 PM UTC-4, Bob Carpenter wrote:

On Jun 23, 2014, at 6:13 PM, Ben Goodrich wrote:

> I think the main problem here is that it is inappropriate to use an iid standard normal prior on a positive_ordered variable.

Inappropriate how?

The elements are not independent due to the ordering, so iid anything is not great. Furthermore, the normal distribution isn't appropriate because each element is bounded on one or both sides, and the normalization factor is a function of the adjacent random variables. And in this case, having a prior mean zero for a positive random variable probably puts too much mass at zero when those diagonal elements need to be bounded away from zero for likelihood-identification.

Ben

[Bob Carpenter]

It makes sense if you take the generative process
to be i.i.d. followed by sorting.

For example, 5 i.i.d. normal variates sorted:


parameters {
vector[5] mu_raw;
}
transformed parameters {
ordered[5] mu;
mu <- sort_asc(mu_raw);
}
model {
mu_raw ~ normal(0,1);
}

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu_raw[1] 0.0 0 1.0 -2.0 -0.7 0.0 0.7 2.0 3240 1
mu_raw[2] 0.0 0 1.0 -1.9 -0.6 0.0 0.7 2.0 3012 1
mu_raw[3] 0.0 0 1.0 -1.9 -0.6 0.0 0.7 1.9 3602 1
mu_raw[4] 0.0 0 1.0 -1.9 -0.6 0.0 0.7 1.9 3478 1
mu_raw[5] 0.0 0 1.0 -1.9 -0.6 0.0 0.7 2.0 3191 1
mu[1] -1.1 0 0.6 -2.5 -1.5 -1.1 -0.7 0.0 2162 1
mu[2] -0.5 0 0.5 -1.6 -0.8 -0.5 -0.1 0.6 2687 1
mu[3] 0.0 0 0.5 -1.0 -0.3 0.0 0.4 1.0 3407 1
mu[4] 0.5 0 0.5 -0.5 0.1 0.5 0.9 1.6 2653 1
mu[5] 1.2 0 0.7 0.0 0.7 1.1 1.6 2.6 2159 1
lp__ -2.4 0 1.6 -6.3 -3.2 -2.1 -1.3 -0.4 1359 1



This is the same result as


parameters {
ordered[5] mu;
}
model {
mu ~ normal(0,1);
}

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] -1.2 0.0 0.7 -2.6 -1.6 -1.1 -0.7 0.0 1060 1
mu[2] -0.5 0.0 0.5 -1.6 -0.8 -0.5 -0.1 0.5 2311 1
mu[3] 0.0 0.0 0.5 -1.0 -0.3 0.0 0.4 1.1 2627 1
mu[4] 0.5 0.0 0.5 -0.5 0.1 0.5 0.8 1.6 2953 1
mu[5] 1.1 0.0 0.7 0.0 0.7 1.1 1.6 2.5 2995 1
lp__ -6.6 0.1 1.7 -10.9 -7.6 -6.3 -5.3 -4.3 883 1



I ran both of these in R using

> fit <- stan("modelname.stan"));

and maximum split R-hat was 1.002 in both cases.

- Bob