R help - Sep 2013 - sparse PCA using nsprcomp package

Hi all, I am using nsprcomp() from nsprcomp package to run sparse PCA. The
output is very much like regular PCA by prcomp() in that it provides
"sdev" for standard deviation of principle components (PC).


For regular PCA by prcomp(), we can easily calculate the percent of total
variance explained by the first k PCs by using cumsum(obj$sdev^2) because these
PCs are independent of each other so you can simply add up the variance of these
PCs. For sparse PCA, as far as I understand, the generated PCs are not
independent of each other anymore, so you can not simply add up variances to
calculate percentage of variance explained by the first k PCs. For example, in
the package of elasticnet where spca() also performs sparse PCA, one of the
output from spca() is "pev" for percent explained variation which is
based on so-called "adjusted" variance that adjusted for the fact that
these variances of PCs are not independent anymore.

My question is for nsprcomp, how can I calculate percent explained variation by
using "sdev" when I know these PCs are not independent of each other?

Thanks!

John
	[[alternative HTML version deleted]]

Hi John

I am currently traveling and have sporadic net access, I therefore can only
answer briefly. It's also quite late, I hope what follows still makes
sense...
> For regular PCA by prcomp(), we can easily calculate the percent of total
variance explained by the first k PCs by using cumsum(obj$sdev^2) because these
PCs are independent of each other so you can simply add up the variance of these
PCs. For sparse PCA, as far as I understand, the generated PCs are not
independent of each other anymore, so you can not simply add up variances to
calculate percentage of variance explained by the first k PCs. For example, in
the package of elasticnet where spca() also performs sparse PCA, one of the
output from spca() is "pev" for percent explained variation which is
based on so-called "adjusted" variance that adjusted for the fact that
these variances of PCs are not independent anymore.
You are correct that measuring explained variance is more involved with sparse
(or non-negative) PCA, because the principal axes no longer correspond to
eigenvectors of the covariance matrix, and are usually not even orthogonal.

The next update for the 'nsprcomp' package is almost done, and one of
the changes will concern the reported standard deviations. In the current
version (0.3), the standard deviations are computed from the scores matrix X*W,
where X is the data matrix and W is the (pseudo-)rotation matrix consisting of
the sparse loadings. Computing variance this way has the advantage that
'sdev' is consistent with the scores matrix, but it has the disadvantage
that some of the explained variance is counted more than once because of the
non-orthogonality of the principal axes. One of the symptoms of this counting is
that the variance of a later component can actually exceed the variance of an
earlier component, which is not possible in regular PCA.

In the new version of the package, 'sdev' will report the _additional_
standard deviation of each component, i.e. the variance not explained by the
previous components. Given a basis of the space spanned by the previous PAs, the
variance of the PC is computed after projecting the current PA to the
ortho-complement space of the basis. This procedure reverts back to standard PCA
if no sparsity or non-negativity constraints are enforced on the PAs.
> My question is for nsprcomp, how can I calculate percent explained
variation by using "sdev" when I know these PCs are not independent of
each other?
The new version of the package will do it for you. Until then, you can use
something like the following function

asdev <- function(X, W) {
    nc <- ncol(W)
    sdev <- numeric(nc)
    Q <- qr.Q(qr(W))
    Xp <- X
    for (cc in seq_len(nc)) {
        sdev[cc] <- sd(Xp%*%W[ , cc])
        Xp <- Xp - Xp%*%Q[ , cc]%*%t(Q[ , cc])   
    }
    return(sdev)
}

to compute the additional variances for given X and W.

The package documentation will explain the above in some more detail, and I will
also have a small blog post which compares the 'nsprcomp' and
'spca' routine from the 'elasticnet' package on the
'marty' data from the EMA package.

Best regards
Christian