<jupyter_start><jupyter_text>While analyzing the PCA plots for protein levels two problems have arisen:- the loadings plot did not show any variables that define PC2  - in order to select variables that explain both PC1 and PC2, a way to extract the contributions of variables of each PC (but not of joint PC1 and PC2) was needed- it was not obvious how are the contributions defined, given different terminologies in useThis notebook addresses these questions using the protein levels data.  Loading the data<jupyter_code>path = 'data/clean/protein/levels.csv'
protein_levels = read_csv(path, index_col=2).drop(['soma_id', 'target_full_name', 'entrez_gene_symbol'], axis=1)
def z_score(v):
    return (v - v.mean()) / v.std()
from numpy import log10
log_matrix = protein_levels.applymap(log10)
z_log_matrix = log_matrix.apply(z_score).apply(z_score, axis=1)
%%R -i z_log_matrix -o pca
pca = prcomp(t(z_log_matrix), scale=T)<jupyter_output><empty_output><jupyter_text>Analysis  A quick look at the distribution of contributions to PC1/2<jupyter_code>%%R
pc_contributions_distribution = function(pca_result, axes, choice='var', labels=FALSE) {
    g = fviz_pca_contrib(pca, choice=choice, axes=axes)
    if(!labels)
        g = g + theme(
            axis.text.x=element_blank(),
            axis.ticks.x=element_blank()
        )
    g
}
%%R -w 300 -h 300 -u px
# "Reference line corresponds to the expected value if the contribution where uniform"
pc_contributions_distribution(pca, axes=1)<jupyter_output><empty_output><jupyter_text>The cut-off point for PC1 is much more to the right (PC1 has a heavy tail), thus more variables contribute to this dimension (compared to PC2, below). Dr Tim noted that similar pattern appears when the data are not centered, as then the first PC is essentially swapping the dataset into the origin. While the z_log_matrix might not be centered on zero (due to the order of z-scores) the `prcomp` function does center it anyway (when scale=T).<jupyter_code>%%R -w 300 -h 300 -u px
pc_contributions_distribution(pca, axes=2)<jupyter_output><empty_output><jupyter_text>While PC2 contributions are much more clearly defined by a small number of proteins, these dimensions do not account for enough variation to be shown on the original biplot.  Getting the top contributors for each of the PCs Loadings are named "rotations" in `prcomp`; columns correspond to the eigenvectors.<jupyter_code>%%R -o loadings -o stdev
loadings = as.data.frame(pca$rotation)
stdev = as.data.frame(pca$sdev)
loadings.head()
from helpers.unsupervised import top_variables
show_source(top_variables)
def top_variables(loadings, pc, n=20):
    return loadings[pc].abs().nlargest(n)<jupyter_output><empty_output><jupyter_text>Quick check to verify that my results are same as of the R-package:<jupyter_code>top_variables(loadings, 'PC1', n=5) # I do not scale by variance - do not need it
%%R -w 250 -h 200 -u px
fviz_pca_contrib(pca, choice="var", axes=1, top=5)<jupyter_output><empty_output><jupyter_text>How are the contributions are defined? I recreate the calculation of contributions step-by-step for benefit of the reader (including me during the write-up): 1) The loadings are normalized (scaled to unit vector):<jupyter_code>from numpy import isclose

assert all(
    isclose(pc_squared, 1) # equal 1 or very, very close
    # (rtol=1e-05, atol=1e-08) - to allow for rounding errors
    for pc_squared in loadings.pow(2).sum()
)<jupyter_output><empty_output><jupyter_text>2) "sdev" gives us the square root of eigenvalues, which can be used to calculate the variance explained by each principal component:<jupyter_code>variance_by_pc = stdev.pow(2)
variance_explained = variance_by_pc / variance_by_pc.sum()
variance_explained.columns = ['variance_explained']
variance_explained.index = 'PC' + variance_explained.index
variance_explained.head()<jupyter_output><empty_output><jupyter_text>In documentation the contributions are defined as: $\text{contrib} = \frac{\text{variable.cos}^2 \cdot 100}{\text{total cos}^2 \text{of the component}}$ As the loadings are cosines of rotation of variables required to form the PCs:<jupyter_code>absolute_contributions = loadings.pow(2) * variance_explained.variance_explained
relative_contributions = 100 * absolute_contributions / absolute_contributions.sum() 
relative_contributions['PC1'].nlargest()<jupyter_output><empty_output><jupyter_text>Which agrees with the values on the plot above. **But how are the contributions for more than one PC calculated?** These are just means of contributions (i.e. means of squared loadings, weighted by variance explained by each of the principal components):<jupyter_code>pc_1_2_contrib = absolute_contributions[['PC1', 'PC2']].mean(axis=1)
(100 * pc_1_2_contrib / pc_1_2_contrib.sum()).nlargest()<jupyter_output><empty_output><jupyter_text>Please, compare with this plot:<jupyter_code>%%R -w 250 -h 200 -u px
fviz_pca_contrib(pca, choice="var", axes=1:2, top=5)<jupyter_output><empty_output><jupyter_text>This is of course different from getting a mean of relative contributions (which would assume equal importance of each of the principal components):<jupyter_code>relative_contributions[['PC1', 'PC2']].mean(axis=1).nlargest()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09381541389153188, 'fraction_numerical': 0.017126546146527116, 'mean_word_length': 5.970822281167108, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '309063', 'license': 'MIT', 'n_tokens_mistral': 1721, 'n_tokens_neox': 1627, 'n_words': 548, 'path': 'exploration/Notes_on_PCA_with_prcomp_and_factoextra.ipynb', 'repo_name': 'krassowski/meningitis-integration', 'size': '5255'}