Looks like you have some numerical precision issues. Why not use the svd function directly? (See below.)
-tgs x <- read.table( textConnection( "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"), header=F, as.is=TRUE) X<-as.matrix(x[,-1]) V <- t( eigen( t(X) %*% X , T )$vectors ) U <- eigen( X %*% t( X ) , T )$vectors D <- diag( sqrt( eigen( X %*% t( X ) , T )$values ) ) X - U %*% D %*% t(V) svd(X)$v - V svd(X)$u - U svd(X)$d - sqrt( eigen( X %*% t( X ) , T )$values ) X-svd(X)$u %*% diag( svd(X)$d ) %*% t( svd(X)$v ) On Fri, May 21, 2010 at 2:25 PM, Julia El-Sayed Moustafa < jelsa...@imperial.ac.uk> wrote: > > 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: > > adjusted_dummyset > [,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=adjusted_dummyset %*% t(adjusted_dummyset) > > 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=t(adjusted_dummyset) %*% adjusted_dummyset > > 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. > > ______________________________________________ > 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. > [[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.
