Clarification regarding rmvnorm in R

Question

Apologies for my limited knowledge of Random Variables. I am working on a problem which requires me to generate a correlation matrix with the absolute pairwise Pearson correlation should be close to 0.4 (off diagonal elements of the correlation matrix). To clarify further let me provide a concrete example. I am following the steps below:

Step 1: Get the sample correlation matrix between 5 products, 
call it sample_corr matrix (just to get the correlation sign +ve or -ve)


> sample_corr
       a          b          c          d          e
a  1.0000000 -0.1907169  1.0000000 -0.5872416  0.8346421
b -0.1907169  1.0000000 -0.1907169  0.6511824 -0.2783893
c  1.0000000 -0.1907169  1.0000000 -0.5872416  0.8346421
d -0.5872416  0.6511824 -0.5872416  1.0000000 -0.6638346
e  0.8346421 -0.2783893  0.8346421 -0.6638346  1.0000000

Now I need to discard the sample values and substitute them by number 0.4 (with maintaining the sign from the above sample matrix).

Step 2: replace values >0 and <0 with 0.4 and -0.4 respectively
sample_corr[sample_corr>0.0]<-0.4
sample_corr[sample_corr< 0]<-(-0.4)
diag(sample_corr)<-1

This replacement of the numbers is not a correct approach as it fails the positive-semi-definite test

> is.positive.semi.definite(sample_corr)
[1] FALSE

Now, despite that error, I use this matrix to generate normal random numbers using rmvnorm with "svd" option (assuming standard deviation as 1)

gen_norm<-rmvnorm(5000,mean = rep(0,5),sigma=sample_corr,method = "svd")

The obtained correlation structure is as follows:

> cor(gen_norm)
           [,1]       [,2]       [,3]       [,4]       [,5]
[1,]  1.0000000 -0.2068077 -0.3214124 -0.2837120 -0.3949594
[2,] -0.2068077  1.0000000 -0.2800788 -0.3087180  0.3548876
[3,] -0.3214124 -0.2800788  1.0000000 -0.2260995  0.4010528
[4,] -0.2837120 -0.3087180 -0.2260995  1.0000000 -0.3444349
[5,] -0.3949594  0.3548876  0.4010528 -0.3444349  1.0000000

Which is not very close to desired correlation but is acceptable to what I was looking for.

My questions are: 1. Why don't I get an error while I am using rmvnorm() to generate the random numbers ? 2. How the correlation structure is corrected ?

Answer 1

Instead of trying to explain the behavior of the fix that rmvnorm() applies, you could instead try and repair the matrix yourself. Here's one way of doing so that appears to work better than the built-in repair within rmvnorm(), at least for this specific case (which is a pretty weird case if you ask me):

################### generate data #######################

# randomly generate a symmetric matrix
C1<-matrix(.8*(runif(25,0,1)>.5)-.4,5,5)
for(j in 1:5){
  for(n in j:5){
    C1[j,n]<-C1[n,j]
  }
}
diag(C1)<-1

# check to see if the random matrix is not positive semi-definite 
# so we can repair it (if its positive semi-definite, there's 
# nothing to repair)
eigen(C1)$values

##  e.g.
#    > eigen(C1)$values
#    [1]  1.9489125  1.4000000  1.4000000  0.6000000 -0.3489125

#################### repair the matrix ###################

# Perform a rank one update with the minimum distance necessary to make the matrix PSM
D<-eigen(C1)$values
V<-eigen(C1)$vectors
C2<- C1+V[,5]%*%t(V[,5])*(.Machine$double.eps-min(D[5],0))

# fix the diagonal
s<- diag(1/sqrt(diag(C2)))
C2<-s%*%C2%*%s

Sidenote: Check for yourself that if you apply this repair procedure on a matrix whose eigenvalues are all non-negative, C1 will equal C2.

######################### compare ########################

C1
C2

#    > C1
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.4000000 -0.4000000 -0.4000000  0.4000000
#    [2,]  0.4000000  1.0000000  0.4000000  0.4000000 -0.4000000
#    [3,] -0.4000000  0.4000000  1.0000000 -0.4000000 -0.4000000
#    [4,] -0.4000000  0.4000000 -0.4000000  1.0000000  0.4000000
#    [5,]  0.4000000 -0.4000000 -0.4000000  0.4000000  1.0000000
#    
#    > C2
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.2861488 -0.3186038 -0.2861488  0.3186038
#    [2,]  0.2861488  1.0000000  0.3186038  0.2861488 -0.3186038
#    [3,] -0.3186038  0.3186038  1.0000000 -0.3186038 -0.4240045
#    [4,] -0.2861488  0.2861488 -0.3186038  1.0000000  0.3186038
#    [5,]  0.3186038 -0.3186038 -0.4240045  0.3186038  1.0000000

# check to see if the procedure you did before will return similar 
# results with the repaired martix
cor(rmvnorm(5000,mean = rep(0,5),sigma=C2,method = "svd"))
#               [,1]       [,2]       [,3]       [,4]       [,5]
#    [1,]  1.0000000  0.2872494 -0.3294268 -0.2792972  0.3085973
#    [2,]  0.2872494  1.0000000  0.3042871  0.3004519 -0.3306211
#    [3,] -0.3294268  0.3042871  1.0000000 -0.3246630 -0.4447972
#    [4,] -0.2792972  0.3004519 -0.3246630  1.0000000  0.3149295
#    [5,]  0.3085973 -0.3306211 -0.4447972  0.3149295  1.0000000

Answer 2

I thought this is worth sharing to people looking for similar solution as I. I found a function called Nearest Positive Definite Matrix nearPD in the Matrixpackage of R. It is what I was looking for in addition to Emily's solution. Here is an excellent blog I found by Nick Higham which explains the function. Lastly, here is the link for the function discreption.

Answer 3

(mvtnorm version 1.1-0)

As of today (April 2020), The function rmvnorm with method = "svd" amounts to

    s. <- svd(sigma)
    if (!all(s.$d >= -sqrt(.Machine$double.eps) * abs(s.$d[1]))) {
        warning("sigma is numerically not positive semidefinite")
    }
    R <- t(s.$v %*% (t(s.$u) * sqrt(pmax(s.$d, 0))))

    retval <- matrix(rnorm(n * ncol(sigma)), nrow = n, byrow = !pre0.9_9994) %*% R
    retval <- sweep(retval, 2, mean, "+")
    colnames(retval) <- names(mean)
    retval

So basically, if I am not mistaken, this means that all negative eigenvalues of sigma are set to zero when its matrix square root (denoted R) is computed.