Correlated random draws with graph structured correlation

Question

I have a problem where I have a graph structure, such that some nodes are connected to other nodes i.e. we have an adjacency matrix of size n*n with a 1 corresponding to a connection and 0 otherwise.

Each of the n nodes in the graph have associated with them two parameters $\alpha$, $\beta$ both distributed jointly normally with correlation $\rho$.

Lets say that for the $i^{th}$ node the parameter values are $\alpha_i$ and $\beta_i$.

Now, I can easily get correlated random draws for each node separately from a bivariate normal distribution.

However, I want that the draws of connected nodes in the graph should be correlated.

Lets say that nodes $i$ and $j$ are connected via an edge in the graph (corresponds to 1 in the adjacency matrix). Then I want that $\alpha_i$ should also be correlated with $\alpha_j$ and $\beta_j$ and similarly $\beta_i$ should be correlated with $\alpha_j$ and $\beta_j$. And as mentioned earlier, $\alpha_i$ and $\beta_i$ are also correlated.

Any thoughts on how to incorporate the graph structured covariance into the random draws?

Answer 1

In general, if you want to draw a $N \times 1$ random vector $x$ with a multivariate normal distribution with mean zero and $N \times N$ variance matrix $\Sigma$, then you do the following:

1. Draw $N$ independent and identically distributed $N(0,1)$ random variables, and stack them up in a vector $y$.
2. Calculate the Cholesky decomposition of $\Sigma$ (the CC' form). I think all statistical programming languages have a function to do this. R has it in the function chol.
3. Multiply $x=Cy$
4. Now, $x$ is distributed normal with mean zero and variance $\Sigma$.

The other part of your question is how to construct $\Sigma$. You did not say anything about the variances of the $\alpha,\beta$, so I am going to assume that each $\alpha$ needs to have variance $\sigma^2_\alpha$ and each $\beta$ needs to have variance $\sigma^2_\beta$.

First, let's be clear about how we are stacking up the $\alpha$s and $\beta$s to make our vector $x$: \begin{align} x = \left( \begin{array}{l} \alpha_1 \\ \beta_1 \\ \alpha_2\\ \beta_2\\ \vdots \\ \alpha_N\\ \beta_N \end{array} \right) \end{align}

Since you want the correlation between each $\alpha,\beta$ pair to be $\rho$, that means that the correlation matrix will start out (before we deal with the adjacency correlations) looking like: \begin{align} \left[ \begin{array}{l l l l l l l l} 1 & \rho & 0 & 0 & 0 & \cdots & 0 & 0\\ \rho & 1 & 0 & 0 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \rho & 0 & \cdots & 0 & 0\\ 0 & 0 & \rho & 1 & 0 & \cdots & 0 & 0\\ \vdots &&&\vdots&&&& \vdots\\ 0 & 0 & 0 & 0 & 0 & \cdots & 1 & \rho\\ 0 & 0 & 0 & 0 & 0 & \cdots & \rho & 1\\ \end{array} \right] \end{align}

You don't say anything about how you want the adjacent nodes' variables to be correlated. So, I'll assume you want the correlation between two adjacent nodes' $\alpha$s to be $\rho_\alpha$, the correlation between two adjacent nodes' $\beta$s to be $\rho_\beta$ and the "cross-correlation" between the $\alpha$ and $\beta$ of two adjacent nodes to be $\rho_{\alpha\beta}$.

Let's imagine that nodes 1 and 2 are adjacent to one another and that nodes 1 and N are adjacent to one another. Then, we want the correlation matrix to look like: \begin{align} \left[ \begin{array}{l l l l l l l l} 1 & \rho & \rho_\alpha & \rho_{\alpha\beta} & 0 & \cdots & \rho_\alpha & \rho_{\alpha\beta}\\ \rho & 1 & \rho_{\alpha\beta} & \rho_\beta & 0 & \cdots & \rho_{\alpha\beta} & \rho_\beta\\ \rho_\alpha & \rho_{\alpha\beta} & 1 & \rho & 0 & \cdots & 0 & 0\\ \rho_{\alpha\beta} & \rho_\beta & \rho & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & \cdots & 0 & 0\\ \vdots &&&\vdots&&&& \vdots\\ 0 & 0 & 0 & 0 & 0 & \cdots & 0 & 0\\ \rho_\alpha & \rho_{\alpha\beta} & 0 & 0 & 0 & \cdots & 1 & \rho\\ \rho_{\alpha\beta} & \rho_\beta & 0 & 0 & 0 & \cdots & \rho & 1\\ \end{array} \right] \end{align}

To get this, we add, to the starting correlation matrix, the matrix: \begin{align} A \otimes \left[\begin{array}{l l} \rho_\alpha & \rho_{\alpha\beta}\\ \rho_{\alpha\beta} & \rho_\beta \end{array} \right] \end{align}

In that expression, $\otimes$ means kroenecker product, and $A$ is the adjacency matrix (I'm assuming symmetric).

That's almost it. One problem may arise. We have not done anything to assure that the correlation matrix we have constructed is positive definite. If there are a lot of adjacencies and if the correlations we have assigned have large absolute values, then there is a chance that the matrix will not be positive definite. You must check for this before calculating the Choleski decomposition, and if the matrix is not positive definite, you have to make the correlations smaller in absolute value.

So, the steps in doing the draw are:

1. Choose the correlations, $\rho$, $\rho_\alpha$, $\rho_\beta$, $\rho_{\alpha\beta}$
2. Choose the variances of $\alpha$ and $\beta$: $\sigma^2_\alpha$, $\sigma^2_\beta$
3. If you don't want $\alpha$ and $\beta$ to have zero means, then choose their means: $\mu_\alpha,\mu_\beta$.
4. Construct the "starting" correlation matrix above
5. Add the kroenecker product matrix to the starting correlation matrix to get the final correlation matrix
6. Check that the final correlation matrix is positive definite. If not, lower the correlations in step 1 and try again.
7. Take the choleski decompostion of the final correlation matrix
8. Draw $2N$ normal (0,1) random variables and stack them up in a vector $y$
9. Multiply the vector $y$ by the choleski decomposition
10. Multiply the $1^{\text{st}}$, $3^\text{rd}$, etc elements of $y$ by $\sqrt{\sigma^2_\alpha}$
11. Multiply the $2^\text{nd}$, $4^\text{th}$, etc elements of $y$ by $\sqrt{\sigma^2}_\beta$
12. Add $\mu_\alpha$ to the $1^{\text{st}}$, $3^\text{rd}$, etc elements of $y$
13. Add $\mu_\beta$ to the $2^\text{nd}$, $4^\text{th}$, etc elements of $y$

I attach an R program which does all this:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/123018/correlated-random-draws-with-graph-structured-correlation

# This program draws four pairs of random variables.  Each is a pair with members named alpha, beta.  
# Each pair of random variable is located at a node in a network.  There is a 4 by 4 adjacency 
# matrix, A, which shows which nodes are adjacent to one another.  All random variables are 
# normal 
# 
# Parameters:
# 
# rho - correlation within nodes between alpha and beta
# rho_a - correlation between alphas of adjacent nodes
# rho_b - correlation between betas of adjacent nodes
# rho_ab - correlation between the alpha at one and the beta at the other adjacent nodes
# 
# mu_a  - mean of alpha
# mu_b  - mean of beta
# sig_a - standard deviation of alpha
# sig_b - standard deviation of beta
# 
# William B Vogt

library(matrixcalc)

set.seed(12344321)


# Set the parameters

rho    <- 0.2
rho_a  <- 0.3
rho_b  <- 0.1
rho_ab <- 0.1

mu_a   <- 0
mu_b   <- 0
sig_a  <- 2
sig_b  <- 3

A <- matrix(c(0,1,1,1,1,0,0,0,1,0,0,0,1,0,0,0),nrow=4,ncol=4)
print(A)

### Construct Sigma, the desired correlation matrix
# first, get the within-node correlations right
Sigma <- diag(4)%x%matrix(c(1,rho,rho,1),nrow=2,ncol=2)
# now, fill in the between-node correlations
Sigma <- Sigma + A%x%matrix(c(rho_a,rho_ab,rho_ab,rho_b),nrow=2,ncol=2)

# check that Sigma is positive definite
is.positive.definite(Sigma)

# factor it
C <- chol(Sigma)  # This function works so that Sigma = C'C

# Draw and transform N(0,1) variates

# We need four nodes times 2 variables, or 8 variables:
x <- rnorm(8,mean=0,sd=1)
y <- t(C)%*%x

# Finally, extract alpha and beta and fix up their means and vars
alpha <- y[seq(from=1,to=7,by=2)]
beta  <- y[seq(from=2,to=8,by=2)]
alpha <- alpha*sig_a + mu_a
beta  <- beta*sig_b + mu_b

# Here is one draw on our network of four nodes:
cbind(alpha,beta)







# Just to check, let's do this draw 10000 times to make sure we are getting the correlations we want
x <- matrix(rnorm(80000,mean=0,sd=1),nrow=8,ncol=10000)
y <- t(C)%*%x
y <- t(y)

# If all is well, these two matrixes should be basically identical
cor(y)
Sigma