HI Christian,
Thanks so much for the detailed explanation! I look forward to the new release
of nsprcomp package! At the meantime, I will use the function below for
calculation of "adjusted" standard deviation. I have 2 more questions, hope you
can shed some lights on:
1). Assume now I can calculate these "adjusted" standard deviation from sparse
PCA, should the percent variation explained by each sparse PC be calculated
using the sum of all these "adjusted" variance (i.e. square of the "adjusted"
standard deviation) as the denominator (then these percent variation explained
will always add up to 1 if all sparse PCs are counted, or using the sum of the
PC variances estimated by REGULAR PCA as the denominator (then, adding up all
PCs may not be equal to 1)?
2). How do you choose the 2 important parameters in nsprcomp(), ncomp and k? If
for example, my regular PCA showed that I need 20 PCs to account for 80% of the
variation in my dataset, does it mean I should set ncomp=20? And then what
about any rules setting the value of "k"?
3). Would you recommend nscumcomp() or nsprcomp() in general?
Thanks so much again,
John
________________________________
From: Christian Sigg <r-h...@sigg-iten.ch>
Cc: r-help@r-project.org
Sent: Thursday, September 5, 2013 2:43 PM
Subject: Re: [R] sparse PCA using nsprcomp package
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
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
