up vote 3 down vote favorite

I'm simulating a multivariate dataset using a modified form of the mvrnorm() function (from the MASS package in R).

The problem is that I'm getting some negative values after the eigenvalue transformation (I assume) because I have many negative correlations and a number of small means.

Is there a special way to deal with this phenomenon? My idea is to just add some factor to all the datapoints (such as the 1st quartile) because I don't care so much what the exact means are, as long as the correlation structure remains intact.

A simple example:

require(MASS)     # For mvrnorm()
require(Matrix) # For nearPD()

corr <- diag(5)
corr[5,1] <- .5
corr[1,5] <- .5
corr[4,2] <- -.5
corr[2,4] <- -.5

set.seed(1000)
mm <- mvrnorm(n=10, mu=rep(1,5), Sigma=nearPD(corr, corr=TRUE)$mat, 
         empirical=TRUE)

As you can see, mm has nonpositive examples. Since I want to model physical measurements, this makes no sense.

edit2: multivariate sampling from log-normal distribution still results in negative examples:

mvrlnorm <- function (n = 1, mu, Sigma, tol = 1e-06, empirical = TRUE) {
    require(Matrix)
    p <- length(mu)
    if (!all(dim(Sigma) == c(p, p)))
        stop("incompatible arguments")
    eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
    ev <- eS$values
    if (!all(ev >= -tol * abs(ev[1L]))) 
        stop("'Sigma' is not positive definite")
      # HERE BE log-normal distribution
    X <- matrix(rlnorm(p * n), nrow=n)
    if (empirical) {
        X <- scale(X, TRUE , FALSE)
        X <- X %*% svd(X, nu = nrow(X), nv = ncol(X))$v
        X <- scale(X, FALSE, TRUE)
    }
    retval <- eS$vectors %*% diag(sqrt(ev), length(ev)) %*% t(eS$vectors)
    retval <- X %*% retval
    retval <- sweep(retval, 2, mu, "+")
    X <- retval


    nm <- names(mu)
    if (is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
        nm <- dn[[1L]]
    dimnames(X) <- list(nm, NULL)
    if (n == 1) 
        drop(X)
    else t(X)
}
share|cite|improve this question
  • 4
    A reproducible example would help. – chl♦ Mar 27 '12 at 20:19
  • To add to chl's remark: it is not clear what you mean by an "eigenvalue transformation"; one cannot tell whether you are referring to negative eigenvalues or negative random values; and--most importantly--the purpose and constraints of the simulation are ambiguous, because if you can arbitrarily add "some factor" to the simulated data--which changes the means--then exactly what statistical characteristics are you attempting to reproduce with this simulation? – whuber♦ Mar 28 '12 at 15:27
  • 1
    Added an example... sorry about that. The point of the simulation is to see if a new correlation metric for a certain kind of data can accurately recover the input correlation matrix (hence, the exact value of the means don't matter much). – zzk Mar 28 '12 at 16:32
  • 1
    What is the basis for choosing $1$ in rep(1,5) instead of, say, $10$, which would make your problem go away? – whuber♦ Mar 28 '12 at 16:47
  • in my actual code, I set the means by sampling from a log-normal distribution (which is a reasonable approx of the actual physical phenomenon - population means of a community of organisms). It's not unusual for the mean of a variable to be close to 1. – zzk Mar 28 '12 at 16:52
up vote 2 down vote

A multivariate normal distribution on $\mathbb{R}^p$ has its support equal to the whole $\mathbb{R}^p$ unless the covariance matrix does not have full rank. Therefore, there always is a positive probability to observe negative components when generating $$ X \sim \mathcal{N}_p(\mu,\Sigma) $$ Thus, to answer your questions:

  1. The issue has nothing to do with the mvnorm() function. It is doing what it is supposed to do.
  2. If you are imposing positivity on your distribution, you cannot use a normal distribution. Use instead a distribution restricted to $\mathbb{R}_+^p$.
share|cite|improve this answer
  • I've edited the source code and replaced the normal distribution with a log normal and I still observe negative components. Want to see the code? – zzk Mar 28 '12 at 16:55
  • 1
    Indeed, I would be most interested in seeing negative log normal realisations! – Xi'an Mar 28 '12 at 16:59
  • added to OP. Running with same seed and correlation matrix results in negative realizations. – zzk Mar 28 '12 at 17:04
  • 1
    @zzk: your code consists in transforming an iid lognormal vector $X$ by linear transforms, i.e. $Y=\mu+A X$. The resulting vector $Y$ is no longer log-normal, since only the location families like the normal distribution remain invariant under linear transforms. And there is nothing paradoxical in observing some components of $Y$ to be negative. – Xi'an Mar 28 '12 at 17:07
  • 1
    @zzk: exactly my point above. While it works for the normal distribution, it does not work for the log-normal distribution! This is of course an interesting problem, how to find a joint distribution on $\mathbb{R}^p_+$ with a given covariance matrix, however I do not know of a closed form solution. – Xi'an Mar 29 '12 at 5:56

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.