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-help@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]]

