首页 文章

不可能从这个相关矩阵创建相关变量?

提问于
浏览
1

我想生成由相关矩阵指定的相关变量 .

首先,我生成相关矩阵:

require(psych)
require(Matrix)

cor.table <- matrix( sample( c(0.9,-0.9) , 2500 , prob = c( 0.8 , 0.2 ) , repl = TRUE ) , 50 , 50 )
k=1
while (k<=length(cor.table[1,])){
    cor.table[1,k]<-0.55
    k=k+1
    }
k=1
while (k<=length(cor.table[,1])){
    cor.table[k,1]<-0.55
    k=k+1
    }   
ind<-lower.tri(cor.table)
cor.table[ind]<-t(cor.table)[ind]
diag(cor.table) <- 1

该相关矩阵不一致,因此不可能进行特征值分解 . 为了使它一致,我使用nearPD:

c<-nearPD(cor.table)

完成后,我生成相关变量:

fit<-principal(c, nfactors=50,rotate="none")
fit$loadings
loadings<-matrix(fit$loadings[1:50, 1:50],nrow=50,ncol=50,byrow=F)
loadings

cases <- t(replicate(50, rnorm(10)) ) 
multivar <- loadings %*% cases
T_multivar <- t(multivar)
var<-as.data.frame(T_multivar)
cor(var)

然而,由此产生的相关性远非我最初指定的任何内容 .

是不可能创建这样的相关性或我做错了什么?

从格雷格·斯诺的评论中可以清楚地看到问题是我的初始相关矩阵是不合理的 .

那么问题是如何使矩阵合理化 . 目标是:

  • 49个变量中的每一个都应该与第一个变量相关联.5 .

  • ~40个变量应该具有相互高的> .8相关性

  • 剩下的~9个变量应该相互之间具有低或负的相关性 .

整个要求不可能吗?

2 回答

  • 1

    尝试使用MASS包中的 mvrnorm 函数,而不是尝试自己构造变量 .

    **编辑

    这是一个正定矩阵(因此它作为相关矩阵)并且接近您的标准,您可以从那里调整值(所有Eigen值都需要为正数,因此您可以看到更改数字的方式如何影响事情):

    cor.mat <- matrix(0.2,nrow=50, ncol=50)
    cor.mat[1,] <- cor.mat[,1] <- 0.55
    cor.mat[2:41,2:41] <- 0.9
    cor.mat[42:50, 42:50] <- 0.25
    diag(cor.mat) <- 1
    
    eigen(cor.mat)$values
    
  • 2

    根据您的上述规范进行的一些数值实验表明,生成的矩阵永远不会(从来没有 - 很好,几乎没有......)是正定的,但它也不会与具有这些值的PD相差很远(使 lcor 低于负值将几乎肯定会让事情变得更糟......)

    rmat <- function(n=49,nhcor=40,hcor=0.8,lcor=0) {
        m <- matrix(lcor,n,n)  ## fill matrix with 'lcor'
        ## select high-cor variables
        hcorpos <- sample(n,size=nhcor,replace=FALSE)
        ## make all of these highly correlated
        m[hcorpos,hcorpos] <- hcor                
        ## compute min real part of eigenvalues
        min(Re(eigen(m,only.values=TRUE)$values))
    }
    set.seed(101)
    r <- replicate(1000,rmat())
    ## NEVER pos definite
    max(r)
    ## [1] -1.069413e-15
    par(las=1,bty="l")
    png("eighist.png")
    hist(log10(abs(r)),breaks=50,col="gray",main="")
    dev.off()
    

    enter image description here

相关问题