[R] R-code to generate random rotation matrix for rotation testing

[Martin Maechler]

>>>>> 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