# [R] Data reconstruction following PCA using Eigen function

Julia El-Sayed Moustafa jelsayed at imperial.ac.uk
Fri May 21 20:25:05 CEST 2010

```Hi all,

As a molecular biologist by training, I'm fairly new to R (and statistics!),
and was hoping for some advice. First of all, I'd like to apologise if my
question is more methodological rather than relating to a specific R
function. I've done my best to search both in the forum and elsewhere but
can't seem to find an answer which works in practice.

I am carrying out principal component analysis in R using the Eigen
function, and have applied this to my dataset and can identify patterns in
my data which result from batch effects. What I'd ultimately like to do is
use PCA as a noise filter to remove the batch effects from my data, and the
way I would like to try and do this is by reconstructing the data, minus the
components (or that part of them) which correlate with batch.

As I have an extremely large dataset, as a first step I am attempting simply
to apply the Eigen function to a dummy dataset consisting of only 10
individuals and ten variables derived from my originial data, then
reconstruct the original dummy data matrix, just to see if I can get the
basics of the method to work first. Here is the dummy matrix (adjusted by
variable mean-subtraction and division by the variable standard deviation
applied to each data point) I am using:

[,1]        [,2]        [,3]        [,4]       [,5]
[,6]        [,7]       [,8]       [,9]      [,10]
Sample1   0.7329881  0.76912670  2.45906143 -0.06411602  1.2427801
0.3785717  2.34508664  1.1043552 -0.1883830  0.6503095
Sample2  -2.0446131  1.72783245 -0.40965941 -0.21713655 -0.2386781
-0.1944390 -0.25181541  0.8882743  0.8404783  0.2531209
Sample3   1.7575174  0.03851425 -0.87537424  1.82811160  0.8342636
-0.8155942 -0.90893068 -0.5529098 -1.6586626  0.1761717
Sample4   0.2731749  0.24045167 -0.39913821 -0.48525909  0.1448994
-2.0173360 -0.01073639 -0.3219478 -1.1536431 -0.2521545
Sample5  -0.2014241 -1.04646151  0.28101160  0.74348390  0.1738312
-0.8431262 -0.08842512  1.2909658  1.2013136  0.6706926
Sample6   0.1743534 -1.70657357 -0.09170187 -0.55605031 -0.2940946
1.4525891 -0.39068509 -0.3373913 -0.0533732  0.9658389
Sample7  -0.8533191 -0.34438091 -1.23890437 -0.77360636  0.5926479
0.7742632 -1.12515017 -0.5720099  0.2243808  0.5420693
Sample8   0.4176988 -0.35906123 -0.07190644  0.90045123 -1.0621902
0.2693762 -0.38033715  0.6267548  0.4767652  0.3012347
Sample9   0.1088066 -0.32197951  0.46665158 -1.72560781  0.7375796
0.2794331  1.00171777 -0.1087306  1.2519195 -0.8848459
Sample10 -0.3651829  1.00253167 -0.12004007  0.34972942 -2.1310389
0.7162622 -0.19072440 -2.0173607 -0.9407956 -2.4224372

I first multiply the matrix by the transposed matrix:

> xxt
Sample1    Sample2    Sample3    Sample4    Sample5    Sample6
Sample7    Sample8    Sample9   Sample10
Sample1  16.045163 -1.1367278 -2.5390033 -1.4762231  1.0158736 -1.8408484
-5.8176751 -1.5163240  3.5304834 -6.2647180
Sample2  -1.136728  9.0985066 -5.2175085 -0.8332460  0.7980575 -3.3607995
1.6342076 -0.3098849 -0.3462389 -0.3263661
Sample3  -2.539003 -5.2175085 12.4738749  3.7747189 -0.9562913 -1.3252885
-0.9175020  0.5849437 -6.0796087  0.2016647
Sample4  -1.476223 -0.8332460  3.7747189  6.1161097 -1.0232149 -3.0984138
-1.1213965 -1.9014901 -1.0503243  0.6134802
Sample5   1.015874  0.7980575 -0.9562913 -1.0232149  6.0758629  0.2183710
-0.9466710  2.1466445 -0.2626433 -7.0659890
Sample6  -1.840848 -3.3607995 -1.3252885 -3.0984138  0.2183710  6.6590617
3.0772441  1.0977948  0.3980588 -1.8251802
Sample7  -5.817675  1.6342076 -0.9175020 -1.1213965 -0.9466710  3.0772441
5.8681617 -0.9215294  0.1646922 -1.0195316
Sample8  -1.516324 -0.3098849  0.5849437 -1.9014901  2.1466445  1.0977948
-0.9215294  3.1757171 -2.2533118 -0.1025599
Sample9   3.530483 -0.3462389 -6.0796087 -1.0503243 -0.2626433  0.3980588
0.1646922 -2.2533118  7.2986178 -1.3997251
Sample10 -6.264718 -0.3263661  0.2016647  0.6134802 -7.0659890 -1.8251802
-1.0195316 -0.1025599 -1.3997251 17.1889249

Then I apply the Eigen function and save the eigenvectors and eigenvalues
eigen(xxt,TRUE)
\$values
[1]  2.693535e+01  1.898845e+01  1.668814e+01  1.184849e+01  8.136377e+00
4.893302e+00  2.007293e+00  4.994285e-01  3.173378e-03 -1.652399e-15

\$vectors
[,1]        [,2]        [,3]        [,4]        [,5]
[,6]         [,7]        [,8]        [,9]      [,10]
[1,]  0.60510924  0.13494101 -0.54092789 -0.11000076  0.11721749
0.37004560 -0.106993925 -0.10503855  0.19436075 -0.3162278
[2,]  0.05581480 -0.39453663  0.09160465  0.67089969  0.08088382
0.32755098 -0.009124933  0.11721316 -0.39379401 -0.3162278
[3,] -0.28666669  0.71121151 -0.01000162  0.03014755 -0.06261445
0.22739798  0.422156670  0.02330019 -0.27677054 -0.3162278
[4,] -0.13469379  0.23709870 -0.15611654  0.27760975 -0.50201617
-0.31821931 -0.584748137  0.12405010  0.11661762 -0.3162278
[5,]  0.26143053  0.12386029  0.29769330  0.19250271  0.27617639
-0.49754797  0.084826842 -0.59908097  0.02670441 -0.3162278
[6,]  0.01918933 -0.04796651  0.35560411 -0.57607725  0.07402744
0.14338749 -0.455866694 -0.01271614 -0.45276432 -0.3162278
[7,] -0.11887886 -0.18220751  0.42748177 -0.07360955 -0.34327294
0.38589551  0.142875991 -0.16600492  0.59142740 -0.3162278
[8,] -0.02725610  0.06425762  0.15301980 -0.01332073  0.53079711
-0.18890512  0.012099700  0.65880817  0.34630946 -0.3162278
[9,]  0.24706215 -0.31886828 -0.13294196 -0.25401667 -0.42308063
-0.38135353  0.482282448  0.24225089 -0.19843310 -0.3162278
[10,] -0.62111060 -0.32779019 -0.48541562 -0.14413474  0.25188193
-0.06825162  0.012492037 -0.28278192  0.04634233 -0.3162278

I also calculate:

> xtx
[,1]       [,2]       [,3]       [,4]       [,5]        [,6]
[,7]       [,8]       [,9]       [,10]
[1,]  9.0000000 -3.1796310  2.0417037  3.9514142  2.7275520 -1.66571377
1.5629736 -0.9104748 -4.8506341  0.68480373
[2,] -3.1796310  9.0000000  1.0981249  0.5494749 -1.2662064 -1.89318349
2.1005566 -0.5052055 -1.7947356 -3.90496189
[3,]  2.0417037  1.0981249  9.0000000 -1.1689542  2.3836861  1.22540724
8.5924385  4.2130304  1.8322104  0.72643120
[4,]  3.9514142  0.5494749 -1.1689542  9.0000000 -1.7132579 -2.51679155
-2.8679239  0.5181511 -3.9536139  0.84095707
[5,]  2.7275520 -1.2662064  2.3836861 -1.7132579  9.0000000 -2.17714376
3.1966738  4.1903147  0.7937017  5.20170437
[6,] -1.6657138 -1.8931835  1.2254072 -2.5167916 -2.1771438  9.00000000
0.3764605 -1.9821441  2.3330854 -0.08204933
[7,]  1.5629736  2.1005566  8.5924385 -2.8679239  3.1966738  0.37646053
9.0000000  3.5708620  1.7809104 -0.28160124
[8,] -0.9104748 -0.5052055  4.2130304  0.5181511  4.1903147 -1.98214408
3.5708620  9.0000000  5.3281685  6.32863320
[9,] -4.8506341 -1.7947356  1.8322104 -3.9536139  0.7937017  2.33308543
1.7809104  5.3281685  9.0000000  2.27959505
[10,]  0.6848037 -3.9049619  0.7264312  0.8409571  5.2017044 -0.08204933
-0.2816012  6.3286332  2.2795951  9.00000000
>
and apply the eigen function to this:

> eigen(xtx,TRUE)
\$values
[1] 2.693535e+01 1.898845e+01 1.668814e+01 1.184849e+01 8.136377e+00
4.893302e+00 2.007293e+00 4.994285e-01 3.173378e-03 6.009067e-16

\$vectors
[,1]        [,2]         [,3]        [,4]         [,5]
[,6]        [,7]        [,8]        [,9]       [,10]
[1,] -0.01604003  0.56323526 -0.190946669 -0.40105857 -0.002543387
0.17417910  0.25192284 -0.59870554  0.18090811 -0.01463338
[2,]  0.08466679 -0.16695878 -0.455330886  0.54953906 -0.002740213
-0.47379253  0.19091718 -0.36818489  0.23997046  0.03943000
[3,] -0.42016626 -0.01672674 -0.438168969 -0.16986221 -0.259624942
0.03774746 -0.20059192  0.33694834  0.50122476 -0.35847979
[4,]  0.17380688  0.46248616 -0.009158497  0.19699882 -0.641525880
-0.08643216  0.32002024  0.38196270 -0.08822512  0.20466364
[5,] -0.38230998  0.27820625  0.061021059  0.01793892  0.556873562
-0.31199982  0.48447946  0.33729239 -0.05542533 -0.11569124
[6,] -0.01078626 -0.35623684  0.086516360 -0.57830505 -0.317806785
-0.56648878  0.22620425 -0.09749859 -0.18851088 -0.11373558
[7,] -0.41357826 -0.06924232 -0.495909637 -0.11281244  0.008420391
0.07292844 -0.09886884  0.01075021 -0.44335225  0.59469674
[8,] -0.48705511  0.07929178  0.158922061  0.33799785 -0.279059362
0.07881225 -0.04227138 -0.30875079 -0.45445853 -0.47881039
[9,] -0.33349167 -0.40387477  0.287225004  0.07874711 -0.163416914
0.37364320  0.52438989 -0.09014761  0.33696991  0.27201980
[10,] -0.34649980  0.24943468  0.446371805  0.03960183 -0.072567359
-0.40852699 -0.43011632 -0.14042978  0.30891048  0.38025985

I create a diagonal matrix containing the square  roots of the eigenvalues
along the diagonal:

> diagonalmatrix_eigenvaluesqrts
V1   V2   V3    V4    V5    V6    V7    V8    V9    V10
1  5.61 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.0000
2  0.00 1.41 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.0000
3  0.00 0.00 1.12 0.000 0.000 0.000 0.000 0.000 0.000 0.0000
4  0.00 0.00 0.00 0.941 0.000 0.000 0.000 0.000 0.000 0.0000
5  0.00 0.00 0.00 0.000 0.727 0.000 0.000 0.000 0.000 0.0000
6  0.00 0.00 0.00 0.000 0.000 0.471 0.000 0.000 0.000 0.0000
7  0.00 0.00 0.00 0.000 0.000 0.000 0.296 0.000 0.000 0.0000
8  0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.246 0.000 0.0000
9  0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.159 0.0000
10 0.00 0.00 0.00 0.000 0.000 0.000 0.000 0.000 0.000 0.0127

I have been operating under the assumption that the original data matrix
could be reconstructed by multiplying the eigenvectors from xxt X a diagonal
matrix with the square roots of the eigenvalues along the diagonal X the
eigenvectors from xtx as follows:

Reconstructed_dummyset=eigenvectors_xxt %*%
as.matrix(diagonalmatrix_eigenvaluesqrts) %*% t(eigenvectors_xtx)
> Reconstructed_dummyset
[,1]         [,2]        [,3]       [,4]        [,5]
[,6]        [,7]        [,8]         [,9]      [,10]
[1,]  0.253193770  0.402535986 -1.14743281  0.5698820 -1.31593686
-0.23278020 -1.10482209 -1.78241789 -1.344859283 -1.4554367
[2,] -0.593623956  0.320035080 -0.30357878 -0.1156343 -0.25992539
-0.25866729 -0.17538542  0.05041073  0.221261845 -0.2675080
[3,]  0.620290857 -0.322491313  0.63113668  0.2551444  0.90025477
-0.36701269  0.61029139  0.90661927  0.226875492  0.6933749
[4,]  0.045394154  0.135401863  0.48897212  0.2772398  0.14530592
-0.12141324  0.33943264  0.53647586  0.002607645  0.4395889
[5,] -0.007416873  0.213143606 -0.90810973  0.2103079 -0.34091303
-0.06529543 -0.82817426 -0.61633082 -0.543192842 -0.2177753
[6,]  0.068080740 -0.533029067 -0.14658587 -0.2022178 -0.09846926
0.29932974 -0.12798021 -0.14611192 -0.016267532  0.1053571
[7,] -0.113178997 -0.309552826  0.18447270 -0.1195559  0.01775153
0.15282100  0.02644529  0.41110216  0.623380590  0.3368958
[8,] -0.076800910 -0.117830736 -0.03085105 -0.1712106  0.39036065
-0.10570034 -0.05526268 -0.08333122 -0.028936234  0.1511139
[9,] -0.187036903  0.212383109 -0.41884937  0.2704928 -0.69211148
0.50295414 -0.45885438 -0.75166730 -0.301317053 -0.6536756
[10,] -0.008314191 -0.002179248  1.66522298 -0.9826679  1.25832941
0.20033204  1.75042615  1.49448058  1.149522805  0.8527934

However as you can see I've clearly got my wires crossed somewhere!

I'd really appreciate any advice you could offer... Thanks!

Julia
--
View this message in context: http://r.789695.n4.nabble.com/Data-reconstruction-following-PCA-using-Eigen-function-tp2226535p2226535.html
Sent from the R help mailing list archive at Nabble.com.

```