# How to weight principal componets by variance?

In accordance with Kaisers rule, only those PCs that had an eigenvalue greater than one were retained for analysis. The PCs that had eigenvalues greater than one were then rotated using a varimax rotation to improve their interpretation and maximize the dispersion of loadings across PCs. These rotated PC scores were weighted by variance and then used to reconstruct the original observations.

I run PCA on 4 variables using the `prcomp` library. All variables were normalized to have a mean of zero and a standard deviation of one (z-score) before the PCA.

``````prc <- prcomp(adpt.pca)
``````

The results are below:

``````Standard deviations (1, .., p=4):
[1] 1.4053803 1.0682221 0.8730050 0.3488124

Rotation (n x k) = (4 x 4):
PC1       PC2         PC3         PC4
income    -0.6183326 0.3817256 -0.18174425 -0.66250989
ndvi       0.1800165 0.7300110  0.65526401  0.07284918
cs_dist    0.3476116 0.5522182 -0.73234784  0.19464803
education -0.6814873 0.1281581 -0.03556312  0.71964281
``````

Then I selected the PCs with eigenvalues > 1 and I performed a `varimax` rotation, like so:

``````varimax2 <- varimax(prc\$rotation[, 1:2])
``````

The results of the varimax rotation:

``````\$loadings

PC1    PC2
income    -0.716  0.122
ndvi      -0.107  0.744
cs_dist    0.115  0.642
education -0.680 -0.137

PC1  PC2
Proportion Var 0.25 0.25
Cumulative Var 0.25 0.50

\$rotmat
[,1]      [,2]
[1,]  0.9269670 0.3751429
[2,] -0.3751429 0.9269670
``````

How can I weight these rotated PC scores by variance in order to reconstruct the original observations?

The complete code:

``````library(data.table)
library(dplyr)

wd <- "path/"

# glimpse(mydt)

# Z-score normalize
adpt.pca <- as.data.frame(scale(mydt), center = TRUE, scale = TRUE)

varimax2 <- varimax(prc\$rotation[, 1:2])
``````

Below a sample of the original dataset (before the Z-score standarization) with 20 rows:

``````dput(mydt)
structure(list(income = c(0.0001063, 0.000106, 6.72e-05, 7.97e-05,
0.0001197, 4.09e-05, 5.17e-05, 0.0001092, 8.62e-05, 7.27e-05,
0.0001034, 0.0001159, 7.24e-05, 9.17e-05, 8.06e-05, 0.0001049,
8.15e-05, 9.05e-05, 0.0001063, 5.99e-05), ndvi = c(0.434779405593872,
0.519024193286896, 0.484442293643951, 0.358367592096329, 0.613705396652222,
0.508738815784454, 0.705485105514526, 0.454894632101059, 0.396738857030869,
0.408085465431213, 0.425091296434402, 0.360570818185806, 0.455742985010147,
0.44114676117897, 0.498669385910034, 0.404618799686432, 0.51068776845932,
0.295410215854645, 0.606453955173492, 0.46584877371788), cs_dist = c(1515.64929199219,
3037.51879882812, 2663.20043945312, 1761.39184570312, 344.697448730469,
252.047805786133, 5528.3486328125, 2387.2802734375, 2771.0546875,
877.851745605469, 1342.23034667969, 3318.9130859375, 1075.06188964844,
5190.70166015625, 739.960021972656, 4005.1572265625, 684.494079589844,
426.935241699219, 1222.70263671875, 2597.5166015625), education = c(0.0001015,
9.71e-05, 6.14e-05, 8.47e-05, 9.97e-05, 5.29e-05, 4.74e-05, 0.0001464,
0.0001042, 7.53e-05, 0.0001143, 9.4e-05, 6.57e-05, 5.52e-05,
7.98e-05, 9.5e-05, 6.98e-05, 9.64e-05, 0.0001063, 6.43e-05)), row.names = c(NA,
-20L), class = c("data.table", "data.frame"), na.action = structure(c(`2402` = 2174L,
`2404` = 2176L), class = "omit"), .internal.selfref = <pointer: 0x0000023dc7801200>)
``````

Note: The results I showed were the outcome when using my complete dataset and not the 20 rows I shared. If needed I can share the entire dataset.

R 4.4.0, RStudio 2024.04.2 Build 764, Windows 11.

I'm not sure yet if this approach is the correct one, so if anyone has any insights it would be nice to comment.

``````# z-score normalization
mydt_normalized <- as.data.frame(scale(mydt))

# PCA
pca_result <- PCA(mydt_normalized, graph = FALSE)

# extract PCs with eigenvalues > 1, this step is to identify the desired # of PCs
eigenvalues <- pca_result\$eig[, "eigenvalue"]
n_components <- sum(eigenvalues > 1)

# PCA with the selected number of components
pca_selected <- PCA(mydt_normalized, ncp = n_components, graph = FALSE)

# varimax rotation

# calculation of rotated PC scores

# weight the rotated PC scores by variance
variance_weights <- eigenvalues[1:n_components] / sum(eigenvalues[1:n_components])
weighted_scores <- sweep(rotated_scores, 2, variance_weights, "*")

# reconstruction the original observations

# Calculate the mean and standard deviation of the original data
original_means <- colMeans(mydt)
original_sds <- apply(mydt, 2, sd)

# Convert the reconstructed data back to original scale
reconstructed_original <- sweep(reconstructed_data, 2, original_sds, "*")
reconstructed_original <- sweep(reconstructed_original, 2, original_means, "+")

# option to display full numbers
options(scipen = 999)

# comparison of the original vs reconstructed obs
# the original observations
mydt[1, ]
income      ndvi  cs_dist education
<num>     <num>    <num>     <num>
1: 0.0001063 0.4347794 1515.649 0.0001015

# the reconstructed observations
income            ndvi          cs_dist          education
<num>           <num>         <num>            <num>
1: 0.0001071987   0.4290230290  1864.3736031922  0.0001081663
``````

I think the term "loadings" may be ambiguous, and could represent either unit-length eigenvectors or these scaled by the variances, which I think are the square roots of the eigenvalues. So the "loadings" (and so scores) may already be "weighted by variance" in a sense; dividing by the sum of the variances would — in usual usage of "weighting" — be the only difference between "scaling" and "weighting". You'd want to know whether "loadings" refers to unit-length eigenvectors or variance-scaled eigenvectors by the particular R functions you're applying.