[R] generate fourth vector based on known correlations
Greg Snow
Greg.Snow at intermountainmail.org
Wed Sep 26 23:10:41 CEST 2007
Here is one approach (someone with better linear algebra skills may be
able to shorten this, some steps could be combined, but I wanted to show
each step):
# create x1-x3
x1 <- rnorm(100, 50, 3)
x2 <- rnorm(100) + x1/5
x3 <- rnorm(100) + x2/5
# find current correlations
cor1 <- cor( cbind(x1,x2,x3) )
cor1
# create 1st version of z
z <- rnorm(100)
# combine in a matrix
m1 <- cbind( x1, x2, x3, z )
# center and scale
m2 <- scale(m1)
# find cholesky decomp
c1 <- chol(var(m2))
# force to be independent
m3 <- m2 %*% solve(c1)
zapsmall(cor(m3)) # check
# create new correlation matrix:
cor2 <- cbind( rbind( cor1, z=c(.5,.3,.1) ), z=c(.5,.3,.1,1) )
# create new matrix
m4 <- m3 %*% chol(cor2)
# uncenter and unscale
m5 <- sweep( m4, 2, attr(m2, 'scaled:scale'), '*')
m5 <- sweep( m5, 2, attr(m2, 'scaled:center'), '+')
# check to see if we succeeded
all.equal( cbind(x1,x2,x3), m5[, 1:3] )
all.equal( cor(m5)[1:3,1:3], cor1 )
all.equal( cor(m5), cor2 )
cor(m5)
Hope this helps,
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at intermountainmail.org
(801) 408-8111
> -----Original Message-----
> From: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] On Behalf Of bkelcey at umich.edu
> Sent: Wednesday, September 26, 2007 2:25 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] generate fourth vector based on known correlations
>
>
>
> I am trying to generate a fourth vector,z, given three known
> and fixed vectors, x1,x2,x3 with corresponding known and
> fixed correlations with themeselves and with z. That is, all
> correlations are known and prespecified. How can I do this?
> Thank you,
> ben
>
> ______________________________________________
> 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.
>
More information about the R-help
mailing list