[R] Generation of correlated variables

Bert Gunter gunter.berton at gene.com
Thu Mar 15 22:22:00 CET 2012


Rui:

Try reading it again. Orthogonal polynomials are generated (subject to
the caveats regarding machine precision stated therein). Note,
especially, the "raw" argument.

Cheers,
Bert

On Thu, Mar 15, 2012 at 1:06 PM, Rui Barradas <rui1174 at sapo.pt> wrote:
> Hello, again.
>
>>
>> This is not really working. Here what I have for the moment.
>>
>
> Right, I've read only half of the description line of function 'poly'. The
> other half states that
>
> "These are all orthogonal to the constant polynomial of degree 0."
>
> But not pairwise orthogonal.
>
> You can look for packages implementing Gram-Schmidt orthogonalization
>
> library(sos)
> r1 <- findFn('Gram')
> r2 <- findFn('Schmidt')
> r1 & r2
>
> or use a (time consuming) trick: form a matrix with the first column as x1,
> and the others as the
> residuals, which are orthogonal, of the regressions of that column with all
> the previous ones.
> Example:
>
> orthogonal <- function(x){
>        x <- cbind(1, x)
>        nc <- ncol(x)
>        res <- matrix(NA, nrow(x), nc - 1)
>        res[, 1] <- x[, 2]
>        if(nc > 2){
>                for(j in 3:nc)
>                        res[, j - 1] <- lm.fit(x[, 1:(j - 1)], x[, j])$residuals
>        }
>        res
> }
>
> set.seed(1)
> X <- matrix(1:10, ncol=1)
> X <- cbind(X, X + rnorm(10), X + rnorm(10))
> cor(X)
>          [,1]      [,2]      [,3]
> [1,] 1.0000000 0.9726364 0.9485793
> [2,] 0.9726364 1.0000000 0.8917848
> [3,] 0.9485793 0.8917848 1.0000000
>
> Y <- orthogonal(X)
> cor(Y)
>              [,1]          [,2]          [,3]
> [1,]  1.000000e+00  2.804232e-17 -3.758236e-17
> [2,]  2.804232e-17  1.000000e+00 -1.492894e-17
> [3,] -3.758236e-17 -1.492894e-17  1.000000e+00
>
> This works but I would use a Gram-Schmidt algorithm.
>
> Rui Barradas
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Generation-of-correlated-variables-tp4475799p4476257.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at 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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm



More information about the R-help mailing list