# How to test alternative dimension reduction strategy to PCA?

I have a system of six co-varying variables which, through principle component analysis, can be reduced to two principle components while retaining >95% of variance explained.

The issue with PCA is that the components themselves do not make physical sense within the context of the data they're intended to fit. Some of the entries in the components (but not all, so no NMF) don't "make sense" to have negative weights. To solve this issue, previous work has derived a method for deriving alternative components through a mixture of expert knowledge and simulated derivation data. While they are presented as alternatives to PCA components, these new components do NOT fit the constraints of PCA--they are not orthogonal and they do not maximize explained variance, but they are much more interpretable.

I'm trying to figure out a way to take these components derived from simulated data, and test how well the work on a "real" data set. My first thinking was to try confirmatory factor analysis--by my thinking, these alternative components derived from theory are essentially latent variables that should underlie the observed real data, and using CFA to test that hypothesis should possibly be a valid test to see if that's the case.

However, this brings up two issues: First, this line of reasoning is not an application of factor analysis I've seen anywhere in literature, so I feel shaky on my theoretical footing. Second, I can't find any example application where CFA has been run in the way I'm trying to run it, where I go into the analysis with two pre-computed factor-weight vectors.

I've tried using the sem package in R, with the following model specification

# First components derived from training/expert analysis

f1 -> x1 , NA ,  0.17
f1 -> x2 , NA ,  0.19
f1 -> x3 , NA ,  0.20

# Second components derived from training/expert analysis

f2 -> x1 , NA ,  0.08
f2 -> x2 , NA , -0.03
f2 -> x4 , NA ,  0.20

# Variances -- Not sure about setting latent variances to 1...
f1 <-> f1 , NA , 1
f2 <-> f2 , NA , 1

x2 <-> x2 , sigma.x2 , NA
x1 <-> x1 , sigma.x1 , NA
x3 <-> x3 , sigma.x3 , NA
x4 <-> x4 , sigma.x4 , NA


When I try to run with this model, I get the following results:

> x.cor<-cor(real.data)
> cfa.model<-specifyModel("sem_model.r")
> cfa.fit<-sem(cfa.model,syn.cor,1500)
Warning message:
In eval(expr, envir, enclos) :
Could not compute QR decomposition of Hessian.
Optimization probably did not converge.
> summary(cfa.fit)
Error in summary.objectiveML(cfa.fit) :
coefficient covariances cannot be computed
In vcov.sem(object, robust = robust, analytic = analytic.se) :
singular Hessian: model is probably underidentified.


Am I doing completely the wrong analysis for what I'm trying to test? If not, what am I doing wrong in my implementation? If so, what is an alternative?

• A quick estimate of the applicability of your expert variables should be possible by looking at their angle with the plane created the 1st and 2nd PCA components. Regarding your implementation of the CFA... You do not really have much free parameters. You are now only estimating the variance of the $x_i$. What happens if you have the variance or regression coefficients for f1 and f2 as free parameters? – Sextus Empiricus Aug 22 '17 at 20:54

To deal with your questions directly (and then to offer some thoughts on the circumstances surrounding your questions):

Am I doing completely the wrong analysis for what I'm trying to test?

If not, what am I doing wrong in my implementation?

This part of your question is, I suspect, a bit off-topic for CV, since it's really about your sem syntax. I'm pretty unfamiliar with sem, but in lavaan (which is the package I know), you would want syntax that would look something like below, where you set up a baseline model where your factor loadings are freely estimated (I've just used a simple one-factor three variable example, but I think you'll get the gist), and then compare against it the fit of a model where you constrain the loadings to the particular values you are interested in (both models would have their latent variance fixed to 1, via the std.lv = T option in the model fitting syntax):

baseline = '
f1 =~ x1 + x2 + x3
'

alt.model = '
f1 =~ .17*x1 + .19*x2 + .20*x4
'

anova(baseline.fit, alt.model.fit)


Other Thoughts