# ---------------------------------------------------------
# Example 1:
# Run RGM based on individual level data with Threshold prior based on the model Y = AY + BX + E
# Data Generation
set.seed(9154)
# Number of data points
n = 10000
# Number of response variables and number of instrument variables
p = 3
k = 4
# Initialize causal interaction matrix between response variables
A = matrix(sample(c(-0.1, 0.1), p^2, replace = TRUE), p, p)
# Diagonal entries of A matrix will always be 0
diag(A) = 0
# Make the network sparse
A[sample(which(A != 0), length(which(A != 0)) / 2)] = 0
# Create D matrix (Indicator matrix where each row corresponds to a response variable
# and each column corresponds to an instrument variable)
D = matrix(0, nrow = p, ncol = k)
# Manually assign values to D matrix
D[1, 1:2] = 1 # First response variable is influenced by the first 2 instruments
D[2, 3] = 1 # Second response variable is influenced by the 3rd instrument
D[3, 4] = 1 # Third response variable is influenced by the 4th instrument
# Initialize B matrix
B = matrix(0, p, k) # Initialize B matrix with zeros
# Calculate B matrix based on D matrix
for (i in 1:p) {
for (j in 1:k) {
if (D[i, j] == 1) {
B[i, j] = 1 # Set B[i, j] to 1 if D[i, j] is 1
}
}
}
# Define Sigma matrix
Sigma = diag(p)
# Compute Mult_Mat
Mult_Mat = solve(diag(p) - A)
# Calculate Variance
Variance = Mult_Mat %*% Sigma %*% t(Mult_Mat)
# Generate instrument data matrix
X = matrix(rnorm(n * k, 0, 1), nrow = n, ncol = k)
# Initialize response data matrix
Y = matrix(0, nrow = n, ncol = p)
# Generate response data matrix based on instrument data matrix
for (i in 1:n) {
Y[i, ] = MASS::mvrnorm(n = 1, Mult_Mat %*% B %*% X[i, ], Variance)
}
# Define a function to create smaller arrowheads
smaller_arrowheads = function(graph) {
igraph::E(graph)$arrow.size = 1 # Adjust the arrow size value as needed
return(graph)
}
# Print true causal interaction matrices between response variables
# and between response and instrument variables
A
B
# Plot the true graph structure between response variables
plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix((A != 0) * 1,
mode = "directed")), layout = igraph::layout_in_circle, main = "True Graph")
# Apply RGM on individual level data for Threshold Prior
Output = RGM(X = X, Y = Y, D = D, prior = "Threshold")
# Get the graph structure between response variables
Output$zAEst
# Get the estimated causal strength matrix between response variables
Output$AEst
# Get the graph structure between response and instrument variables
Output$zBEst
# Get the estimated causal strength matrix between response and instrument variables
Output$BEst
# Get the estimated causal network between responses
Output$Graph
# Plot the estimated causal network
plot(Output$Graph, main = "Estimated Causal Network")
# Plot posterior log-likelihood
plot(Output$LLPst, type = 'l', xlab = "Number of Iterations", ylab = "Log-likelihood")
# -----------------------------------------------------------------
# Example 2:
# Run RGM based on Syy, Syx and Sxx with Spike and Slab prior based on the model Y = AY + BX + E
# Data Generation
set.seed(9154)
# Number of data points
n = 10000
# Number of response variables and number of instrument variables
p = 3
k = 4
# Initialize causal interaction matrix between response variables
A = matrix(sample(c(-0.1, 0.1), p^2, replace = TRUE), p, p)
# Diagonal entries of A matrix will always be 0
diag(A) = 0
# Make the network sparse
A[sample(which(A!=0), length(which(A!=0))/2)] = 0
# Create D matrix (Indicator matrix where each row corresponds to a response variable
# and each column corresponds to an instrument variable)
D = matrix(0, nrow = p, ncol = k)
# Manually assign values to D matrix
D[1, 1:2] = 1 # First response variable is influenced by the first 2 instruments
D[2, 3] = 1 # Second response variable is influenced by the 3rd instrument
D[3, 4] = 1 # Third response variable is influenced by the 4th instrument
# Initialize B matrix
B = matrix(0, p, k) # Initialize B matrix with zeros
# Calculate B matrix based on D matrix
for (i in 1:p) {
for (j in 1:k) {
if (D[i, j] == 1) {
B[i, j] = 1 # Set B[i, j] to 1 if D[i, j] is 1
}
}
}
Sigma = diag(p)
Mult_Mat = solve(diag(p) - A)
Variance = Mult_Mat %*% Sigma %*% t(Mult_Mat)
# Generate instrument data matrix
X = matrix(rnorm(n * k, 0, 1), nrow = n, ncol = k)
# Initialize response data matrix
Y = matrix(0, nrow = n, ncol = p)
# Generate response data matrix based on instrument data matrix
for (i in 1:n) {
Y[i, ] = MASS::mvrnorm(n = 1, Mult_Mat %*% B %*% X[i, ], Variance)
}
# Calculate summary level data
Syy = t(Y) %*% Y / n
Syx = t(Y) %*% X / n
Sxx = t(X) %*% X / n
# Print true causal interaction matrices between response variables
# and between response and instrument variables
A
B
# Plot the true graph structure between response variables
plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix(((A != 0) * 1),
mode = "directed")), layout = igraph::layout_in_circle, main = "True Graph")
# Apply RGM on summary level data for Spike and Slab Prior
Output = RGM(Syy = Syy, Syx = Syx, Sxx = Sxx,
D = D, n = 10000, prior = "Spike and Slab")
# Get the graph structure between response variables
Output$zAEst
# Get the estimated causal strength matrix between response variables
Output$AEst
# Get the graph structure between response and instrument variables
Output$zBEst
# Get the estimated causal strength matrix between response and instrument variables
Output$BEst
# Get the estimated causal network between responses
Output$Graph
# Plot the estimated causal network
plot(Output$Graph, main = "Estimated Causal Network")
# Plot posterior log-likelihood
plot(Output$LLPst, type = 'l', xlab = "Number of Iterations", ylab = "Log-likelihood")
# -----------------------------------------------------------------
# Example 3:
# Run RGM based on Sxx, Beta and SigmaHat with Spike and Slab prior
# based on the model Y = AY + BX + E
# Data Generation
set.seed(9154)
# Number of datapoints
n = 10000
# Number of response variables and number of instrument variables
p = 3
k = 4
# Initialize causal interaction matrix between response variables
A = matrix(sample(c(-0.1, 0.1), p^2, replace = TRUE), p, p)
# Diagonal entries of A matrix will always be 0
diag(A) = 0
# Make the network sparse
A[sample(which(A!=0), length(which(A!=0))/2)] = 0
# Create D matrix (Indicator matrix where each row corresponds to a response variable
# and each column corresponds to an instrument variable)
D = matrix(0, nrow = p, ncol = k)
# Manually assign values to D matrix
D[1, 1:2] = 1 # First response variable is influenced by the first 2 instruments
D[2, 3] = 1 # Second response variable is influenced by the 3rd instrument
D[3, 4] = 1 # Third response variable is influenced by the 4th instrument
# Initialize B matrix
B = matrix(0, p, k) # Initialize B matrix with zeros
# Calculate B matrix based on D matrix
for (i in 1:p) {
for (j in 1:k) {
if (D[i, j] == 1) {
B[i, j] = 1 # Set B[i, j] to 1 if D[i, j] is 1
}
}
}
Sigma = diag(p)
Mult_Mat = solve(diag(p) - A)
Variance = Mult_Mat %*% Sigma %*% t(Mult_Mat)
# Generate DNA expressions
X = matrix(rnorm(n * k, 0, 1), nrow = n, ncol = k)
# Initialize response data matrix
Y = matrix(0, nrow = n, ncol = p)
# Generate response data matrix based on instrument data matrix
for (i in 1:n) {
Y[i, ] = MASS::mvrnorm(n = 1, Mult_Mat %*% B %*% X[i, ], Variance)
}
# Centralize Data
Y = t(t(Y) - colMeans(Y))
X = t(t(X) - colMeans(X))
# Calculate Sxx
Sxx = t(X) %*% X / n
# Generate Beta matrix and SigmaHat
Beta = matrix(0, nrow = p, ncol = k)
SigmaHat = matrix(0, nrow = p, ncol = k)
for (i in 1:p) {
for (j in 1:k) {
fit = lm(Y[, i] ~ X[, j])
Beta[i, j] = fit$coefficients[2]
SigmaHat[i, j] = sum(fit$residuals^2) / n
}
}
# Print true causal interaction matrices between response variables
# and between response and instrument variables
A
B
# Plot the true graph structure between response variables
plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix(((A != 0) * 1),
mode = "directed")), layout = igraph::layout_in_circle, main = "True Graph")
# Apply RGM based on Sxx, Beta and SigmaHat for Spike and Slab Prior
Output = RGM(Sxx = Sxx, Beta = Beta, SigmaHat = SigmaHat,
D = D, n = 10000, prior = "Spike and Slab")
# Get the graph structure between response variables
Output$zAEst
# Get the estimated causal strength matrix between response variables
Output$AEst
# Get the graph structure between response and instrument variables
Output$zBEst
# Get the estimated causal strength matrix between response and instrument variables
Output$BEst
# Get the estimated causal network between responses
Output$Graph
# Plot the estimated causal network
plot(Output$Graph, main = "Estimated Causal Network")
# Plot posterior log-likelihood
plot(Output$LLPst, type = 'l', xlab = "Number of Iterations", ylab = "Log-likelihood")
Run the code above in your browser using DataLab