### Install MR.RGM Package Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Instructions to install the MR.RGM R package from CRAN. ```r install.packages("MR.RGM") ``` -------------------------------- ### RGM Function Example: Simulate Data and Run Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Demonstrates the usage of the RGM function with simulated individual-level data. It includes data simulation for response variables (Y), instrument variables (X), and causal interaction matrices (A, B), followed by calling the RGM function with a 'Threshold' prior. ```r # Model: Y = AY + BX + E # Set seed set.seed(9154) # Number of data points n = 10000 # Number of response variables and number of instrument variables p = 5 k = 6 # 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 D[4, 5] = 1 # Fourth response variable is influenced by the 5th instrument D[5, 6] = 1 # Fifth response variable is influenced by the 6th 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 } } } # Create variance-covariance matrix Sigma = 1 * diag(p) Mult_Mat = solve(diag(p) - A) Variance = Mult_Mat %*% Sigma %*% t(Mult_Mat) # Generate instrument data matrix X = matrix(runif(n * k, 0, 5), 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) } # Print true causal interaction matrices between response variables and between response and instrument variables A #> [,1] [,2] [,3] [,4] [,5] #> [1,] 0.0 -0.1 0.0 0.0 0.1 #> [2,] 0.1 0.0 -0.1 0.1 0.1 #> [3,] 0.0 -0.1 0.0 0.0 0.1 #> [4,] 0.0 -0.1 0.0 0.0 0.0 #> [5,] 0.0 0.1 0.0 0.0 0.0 B #> [,1] [,2] [,3] [,4] [,5] [,6] #> [1,] 1 1 0 0 0 0 #> [2,] 0 0 1 0 0 0 #> [3,] 0 0 0 1 0 0 #> [4,] 0 0 0 0 1 0 #> [5,] 0 0 0 0 0 1 # Apply RGM on individual level data with Threshold prior Output1 = RGM(X = X, Y = Y, D = D, prior = "Threshold") ``` -------------------------------- ### R Simulation: High-Dimensional IVs with PCA Reduction Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Demonstrates an R simulation setup with numerous instrumental variables (IVs) and dimensionality reduction via PCA. It generates data, computes summary statistics, and prepares a compact matrix for RGM analysis, requiring MASS and igraph libraries. ```R # Load necessary libraries library(MASS) library(igraph) #> #> Attaching package: 'igraph' #> The following objects are masked from 'package:stats': #> #> decompose, spectrum #> The following object is masked from 'package:base': #> #> union # Set seed for reproducibility set.seed(9154) # Number of data points n = 10000 # Number of response variables p = 5 # Number of SNPs per response variable num_snps_per_y = 100 # Total number of SNPs k = num_snps_per_y * p # Initialize causal interaction matrix between response variables A = matrix(sample(c(-0.1, 0.1), p^2, replace = TRUE), p, p) diag(A) = 0 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) # Assign values to D matrix using a loop for (run in 1:p) { D[run, ((run - 1) * 100 + 1) : (run * 100)] = 1 } # 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 } } } # Calculate Variance-Covariance matrix Sigma = diag(p) Mult_Mat = solve(diag(p) - A) Variance = Mult_Mat %*% Sigma %*% t(Mult_Mat) # Generate instrument data matrix (X) X = matrix(rnorm(n * k, 0, 1), nrow = n, ncol = k) # Initialize response data matrix (Y) 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 # Perform PCA for each response variable to get top 20 PCs top_snps_list = list() for (i in 1:p) { X_sub = X[, (num_snps_per_y * (i - 1) + 1):(num_snps_per_y * i)] pca = prcomp(X_sub, center = TRUE, scale. = TRUE) top_20_pcs = pca$x[, 1:20] top_snps_list[[i]] = top_20_pcs } # Combine the top PCs from all response variables compact_X = do.call(cbind, top_snps_list) # Calculate summary level data based on compact_X Sxx_compact = t(compact_X) %*% compact_X / n Syx_compact = t(Y) %*% compact_X / n # Create D_New D_New = matrix(0, nrow = p, ncol = 20 * p) # Assign values to D matrix using a loop for (run in 1:p) { D_New[run, ((run - 1) * 20 + 1) : (run * 20)] = 1 } # Apply RGM on summary level data for Spike and Slab Prior using the compact_X matrix Output = RGM(Syy = Syy, Syx = Syx_compact, Sxx = Sxx_compact, D = D_New, n = n, prior = "Spike and Slab") ``` -------------------------------- ### Load MR.RGM Library Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Loads the MR.RGM package into the R workspace for use. ```r library("MR.RGM") ``` -------------------------------- ### Initialize and Plot Subgraph Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Initializes a matrix for a subgraph and sets specific directed edges. It then plots this subgraph using igraph, visualizing the initial causal structure. ```r # Start with a random subgraph Gamma = matrix(0, nrow = p, ncol = p) Gamma[5, 2] = Gamma[3, 5] = Gamma[2, 3] = 1 # Plot the subgraph to get an idea about the causal network plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix(Gamma, mode = "directed")), layout = igraph::layout_in_circle, main = "Subgraph") ``` -------------------------------- ### Execute NetworkMotif Function Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Collects posterior probability arrays (GammaPst) from various outputs and uses them with an initial graph structure (Gamma) to execute the NetworkMotif function. This function likely calculates motif occurrences or probabilities. ```r # Store the GammaPst arrays from outputs GammaPst1 = Output1$GammaPst GammaPst2 = Output2$GammaPst GammaPst3 = Output3$GammaPst # Get the posterior probabilities of Gamma with these GammaPst matrices NetworkMotif(Gamma = Gamma, GammaPst = GammaPst1) NetworkMotif(Gamma = Gamma, GammaPst = GammaPst2) NetworkMotif(Gamma = Gamma, GammaPst = GammaPst3) ``` -------------------------------- ### Calculate Beta and SigmaHat for RGM Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Centralizes the input data Y and X by subtracting column means. Then, it calculates the Sxx matrix and fits linear models (lm) for each pair of response (Y[,i]) and predictor (X[,j]) to estimate the Beta coefficients and SigmaHat values. ```r # 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 } } ``` -------------------------------- ### Plot Log-Likelihoods Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Plots the log-likelihood values over iterations using a line graph. This helps visualize the convergence of the model. Requires a log-likelihood vector (LLPst) from the output object. ```r plot(Output1$LLPst, type = 'l', xlab = "Iterations", ylab = "Log-likelihood") plot(Output2$LLPst, type = 'l', xlab = "Iterations", ylab = "Log-likelihood") plot(Output3$LLPst, type = 'l', xlab = "Iterations", ylab = "Log-likelihood") ``` -------------------------------- ### Apply RGM with Spike and Slab Prior (Beta/SigmaHat) Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Applies the RGM function using calculated Beta and SigmaHat matrices, along with Sxx, and a specified prior ('Spike and Slab'). This estimates the causal interaction matrix (AEst) and network structure (zAEst). ```r Output3 = RGM(Sxx = Sxx, Beta = Beta, SigmaHat = SigmaHat, D = D, n = 10000, prior = "Spike and Slab") ``` -------------------------------- ### Apply RGM with Spike and Slab Prior (Summary Data) Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Applies the RGM function using pre-calculated summary statistics (Syy, Syx, Sxx) and a specified prior ('Spike and Slab'). This estimates the causal interaction matrix (AEst) and network structure (zAEst). ```r Output2 = RGM(Syy = Syy, Syx = Syx, Sxx = Sxx, D = D, n = 10000, prior = "Spike and Slab") ``` -------------------------------- ### Calculate Summary Statistics for RGM Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Calculates covariance matrices (Syy, Syx, Sxx) from input data Y and X, normalized by n. These are used as inputs for the RGM function. ```r Syy = t(Y) %*% Y / n Syx = t(Y) %*% X / n Sxx = t(X) %*% X / n ``` -------------------------------- ### Display Causal Graph Structure Matrix Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Retrieves and displays the estimated graph structure matrix (zBEst) between response and instrument variables from different output objects. This matrix typically indicates the presence (1) or absence (0) of causal links. ```r Output1$zBEst Output2$zBEst Output3$zBEst ``` -------------------------------- ### Display Estimated Causal Interaction Matrices Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Prints the estimated causal interaction matrices, AEst and zAEst, which represent the estimated causal effects and their binary presence, respectively. These matrices are crucial for understanding direct and indirect causal relationships. ```R # Print estimated causal interaction matrices Output$AEst #> [,1] [,2] [,3] [,4] [,5] #> [1,] 0.0000000000 -0.10304239 0.012795897 0.006911765 0.09119095 #> [2,] 0.1102430697 0.00000000 -0.101718794 0.095093965 0.10762267 #> [3,] -0.0039579652 -0.10801784 0.000000000 -0.008388722 0.10255348 #> [4,] -0.0013919956 -0.09076408 0.011536478 0.000000000 0.01075220 #> [5,] -0.0005533104 0.10848253 -0.002160112 -0.011725145 0.00000000 Output$zAEst #> [,1] [,2] [,3] [,4] [,5] #> [1,] 0 1 0 0 1 #> [2,] 1 0 1 1 1 #> [3,] 0 1 0 0 1 #> [4,] 0 1 0 0 0 #> [5,] 0 1 0 0 0 ``` -------------------------------- ### Plot Causal Networks Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Generates and displays two causal network plots: the true causal network and the estimated causal network. This visualization helps in comparing the model's performance against the ground truth. ```R # Create a layout for multiple plots par(mfrow = c(1, 2)) # Plot the true causal network plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix((A != 0) * 1, mode = "directed")), layout = igraph::layout_in_circle, main = "True Causal Network") # Plot the estimated causal network plot(Output$Graph, main = "Estimated Causal Network") ``` ``` -------------------------------- ### Display Causal Interaction Matrix Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Retrieves and displays the estimated causal interaction matrix (BEst) between response and instrument variables from different output objects. This matrix represents the estimated coefficients. ```r Output1$BEst Output2$BEst Output3$BEst ``` -------------------------------- ### Plot Estimated Causal Network Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Plots the estimated causal network graph from an output object. Requires a graph object stored within the output, typically accessed via '$Graph'. ```r plot(Output1$Graph, main = "Estimated Causal Network") ``` -------------------------------- ### Plotting Causal Network Structure Source: https://github.com/bitansa/mr.rgm/blob/master/README.md Defines a helper function to adjust arrowhead size and then plots the true causal network and the estimated causal network using igraph. It converts adjacency matrices to graph objects and applies a circular layout. ```r # Define a function to create smaller arrowheads smaller_arrowheads <- function(graph) { igraph::E(graph)$arrow.size = 0.60 # Adjust the arrow size value as needed return(graph) } # Create a layout for multiple plots par(mfrow = c(1, 2)) # Plot the true causal network plot(smaller_arrowheads(igraph::graph_from_adjacency_matrix((A != 0) * 1, mode = "directed")), layout = igraph::layout_in_circle, main = "True Causal Network") ``` === COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.