Martin Krautschke
2010-Dec-27 21:47 UTC
[R] R-code to generate random rotation matrix for rotation testing
Dear list, I am looking for an implementation of random rotation matrix generation in R to do a rotation test: I want to use the matrices to create random multivariate normal matrices with common covariance structure and mean based on an observed data matrix. The rRotationMatrix-function in the mixAK-package is an option, but as far as I can tell I need to draw rotation matrices with determinant -1 as well. Roast and Romer in the limma-bioconductor package appear to have implemented something similar, which seems not to be general enough for my purposes, however. Inspired by the code in the ffmanova-rotationtest function I thought of the following, but it appears to me that there only the covariance, not the mean, is preserved: ##### # a given Y has independent, multivariate normal rows library(mvtnorm) Y <- rmvnorm(4,mean=1:10,sigma=diag(1:10)+3) # Generation of a set of random matrices Z for (i in 1:10) { # R is random matrix of independent standard-normal entries R <- matrix(rnorm(16),ncol=4) R <- qr.Q(qr(R, LAPACK = TRUE)) # Z shall be a random matrix with the same mean and covariance structure as Y Z <- crossprod(R,Y) } ##### A suggestion for the procedure exists (in Dorum et al. http://www.bepress.com/sagmb/vol8/iss1/art34/ , end of chapter 2.1), but a hint to a (fast) implementation would be greatly appreciated. Best regards and a happy new year, Martin Krautschke ----------------------- Martin Krautschke Student at University of Vienna -- Sicherer, schneller und einfacher. Die aktuellen Internet-Browser - jetzt kostenlos herunterladen! http://portal.gmx.net/de/go/atbrowser
Martin Maechler
2010-Dec-29 11:34 UTC
[R] R-code to generate random rotation matrix for rotation testing
>>>>> Martin Krautschke <m.krautschke at gmx.at> >>>>> on Mon, 27 Dec 2010 22:47:26 +0100 writes:> I am looking for an implementation of random rotation > matrix generation in R to do a rotation test: I want to > use the matrices to create random multivariate normal > matrices with common covariance structure and mean based > on an observed data matrix. > The rRotationMatrix-function in the mixAK-package is an > option, but as far as I can tell I need to draw rotation > matrices with determinant -1 as well. Roast and Romer in > the limma-bioconductor package appear to have implemented > something similar, which seems not to be general enough > for my purposes, however. Inspired by the code in the > ffmanova-rotationtest function I thought of the following, > but it appears to me that there only the covariance, not > the mean, is preserved: > ##### > # a given Y has independent, multivariate normal rows > library(mvtnorm) > Y <- rmvnorm(4,mean=1:10,sigma=diag(1:10)+3) > # Generation of a set of random matrices Z > for (i in 1:10) { > # R is random matrix of independent standard-normal entries > R <- matrix(rnorm(16),ncol=4) > R <- qr.Q(qr(R, LAPACK = TRUE)) > # Z shall be a random matrix with the same mean and covariance structure as Y > Z <- crossprod(R,Y) > } > ##### > A suggestion for the procedure exists (in Dorum et al. http://www.bepress.com/sagmb/vol8/iss1/art34/ , end of chapter 2.1), but a hint to a (fast) implementation would be greatly appreciated. > Best regards and a happy new year, > Martin Krautschke > ----------------------- > Martin Krautschke > Student at University of Vienna and this is not a home work problem? Just in case, I don't give you the complete solution in R, but in words : Think geometrically: Rotation in the above sense only preserves the mean when that is the zero vector. Consequently: Your procedure must rather be 1) Y0 <- Y - mY 2) Z0 <- Q' %*% Y0 3) Z <- Z0 + mY and to make this work with data matrices Y, Z, the mean vector mY must either be a matrix with constant columns or the result of as.vector()ing such a matrix. Regards, Martin Maechler, ETH Zurich