Building association network and sampling in multivariate normal distribution

We choose a linear network for the support of the precision matrix (inverse of the covariance matrix).

p <- 6
g <- make_lattice(length = p,dim = 1)
V(g)$name <- LETTERS[1:p]
#g is an object of class 'igraph'
class(g)

## [1] "igraph"

print(g)

## IGRAPH 3000480 UN-- 6 5 -- Lattice graph
## + attr: name (g/c), dimvector (g/n), nei (g/n), mutual (g/l), circular
## | (g/l), name (v/c)
## + edges from 3000480 (vertex names):
## [1] A--B B--C C--D D--E E--F

#using igraph plot
plot(g)

#a nicer representation with ggnet2
ggnet(g,label = TRUE,color = "grey")

We add edge weights that are non-zero values of the precision matrix (inverse of the variance-covariance matrix). These are association strengths.

#Adding weights
E(g)$weight <- c(1,1,0.2,1,1)
#Get the adjacency matrix of g
A <- get.adjacency(g,attr = "weight")
#adding variance (diagonal term)
diag(A) = 2

Careful: Remember that precision matrix coefficients and partial correlation coefficients (the ones of interest for interpretation) are of opposite signs.

A = as.matrix(A)
p_cor_theo <- - A/(sqrt(diag(A) %*% t(diag(A))))
diag(p_cor_theo) = 1
print(A)

##   A B   C   D E F
## A 2 1 0.0 0.0 0 0
## B 1 2 1.0 0.0 0 0
## C 0 1 2.0 0.2 0 0
## D 0 0 0.2 2.0 1 0
## E 0 0 0.0 1.0 2 1
## F 0 0 0.0 0.0 1 2

print(p_cor_theo)

##      A    B    C    D    E    F
## A  1.0 -0.5  0.0  0.0  0.0  0.0
## B -0.5  1.0 -0.5  0.0  0.0  0.0
## C  0.0 -0.5  1.0 -0.1  0.0  0.0
## D  0.0  0.0 -0.1  1.0 -0.5  0.0
## E  0.0  0.0  0.0 -0.5  1.0 -0.5
## F  0.0  0.0  0.0  0.0 -0.5  1.0

We now represent the partial correlation matrix as a network using ggnet2

#building igraph objects
g_pcor <- graph_from_adjacency_matrix(abs(p_cor_theo),weighted = TRUE)
# convert it to network format
g_pcor = asNetwork(g_pcor)
ggnet2(g_pcor,label = TRUE,color = "grey",edge.size = "weight",edge.color = "blue")

## Warning in ggnet2(g_pcor, label = TRUE, color = "grey", edge.size = "weight", :
## ggnet2 does not know how to handle self-loops

We now draw samples in a multivariate normal distribution of precision matrix A.

#choose mean
mu <- rep(1,p)
X <- rmvnp(n=1000, mu, Omega = as.matrix(A))

We now compute the sample correlation matrix

S_cor = cor(X)
print(S_cor)

##             A           B           C           D           E           F
## A  1.00000000 -0.59186377  0.29186181 -0.01512471  0.03738444 -0.06132183
## B -0.59186377  1.00000000 -0.53195106  0.05881456 -0.06283624  0.08135608
## C  0.29186181 -0.53195106  1.00000000 -0.11885419  0.06431436 -0.08810876
## D -0.01512471  0.05881456 -0.11885419  1.00000000 -0.59125126  0.38895696
## E  0.03738444 -0.06283624  0.06431436 -0.59125126  1.00000000 -0.58708085
## F -0.06132183  0.08135608 -0.08810876  0.38895696 -0.58708085  1.00000000

The sign of correlation between neighbors is correct but some other correlations are spurious (e.g. A-C). We can build a sample correlation network using a threshold.

thrs <- 0.35
S_cor_thrs <- S_cor
S_cor_thrs[which(abs(S_cor)<thrs)] = 0

#represent marginal correlation network
diag(S_cor_thrs) = 0
g_cor_thrs <- graph_from_adjacency_matrix(abs(S_cor_thrs),weighted = TRUE)
g_cor_thrs = asNetwork(g_cor_thrs)
ggnet2(g_cor_thrs,label = TRUE,edge.size = "weight",edge.color = "blue")

The marginal correlation network is disconnected with this threshold value. What happens if we decrease the threshold ?

thrs = 0.3
S_cor_thrs <- S_cor
S_cor_thrs[which(abs(S_cor)<thrs)] = 0
diag(S_cor_thrs) = 0
#representing the network with ggnet
g_cor_thrs <- graph_from_adjacency_matrix(abs(S_cor_thrs),weighted = TRUE)
g_cor_thrs = asNetwork(g_cor_thrs)
ggnet2(g_cor_thrs,label = TRUE,edge.size = "weight",edge.color = "blue")

Network is still disconnected ! Moreover we have spurious associations. The two highly correlated components have strong marginal correlations inside them.

Inference with the graphical lasso

We saw that marginal correlations do not recover the original association network. We now aim at recovering the correct partial correlation network using graphical lasso.

S <- cov(X)
#applying graphical lasso with fixed penalty parameter
mod1 <- glasso(S,rho = 0.05)
pcor_estimated <- - mod1$wi/sqrt(diag(mod1$wi) %*% t(diag(mod1$wi)))
diag(pcor_estimated) = 1
rownames(pcor_estimated) <- colnames(pcor_estimated) <- colnames(S)
#representing the partial correlation network
g_pcor <- graph_from_adjacency_matrix(abs(pcor_estimated),weighted = TRUE)
g_pcor = asNetwork(g_pcor)
ggnet2(g_pcor,label = TRUE,edge.size = "weight",edge.color = "blue")

## Warning in ggnet2(g_pcor, label = TRUE, edge.size = "weight", edge.color =
## "blue"): ggnet2 does not know how to handle self-loops

We now run glasso algorithm with a range of penalty parameters to represent coefficient path.

rho_seq <- seq(0.0001,1,length.out = 100)
wi_melted_list <- list()
#for loop over penalty parameters
for(rho in rho_seq){
mod_loc <- glasso(S,rho)
wi_loc <- mod_loc$wi
rownames(wi_loc) <- colnames(wi_loc) <- colnames(S)
#creating data frame for plotting
wi_melted <- melt(wi_loc)
wi_melted <- wi_melted[which(wi_melted[,1] != wi_melted[,2]),]
wi_melted <- cbind(wi_melted,rep(rho,ncol(wi_melted)),paste0(wi_melted[,1],"|",wi_melted[,2]))
colnames(wi_melted) = c("Node1","Node2","Coeff","Penalty","Edge")
wi_melted_list = c(wi_melted_list,list(wi_melted))
}
#merging results
wi_melted_df = do.call(rbind,wi_melted_list)

#representing coefficient path using ggplot2
ggplot(wi_melted_df,aes(x = Penalty, y = Coeff, group = Edge, colour = Edge)) + 
  geom_line(show.legend = FALSE) +  coord_trans(x="log10") + theme_bw()

Decreasing rho increase the number of non-null coefficients, that is the number of edges of the partial correlation network. How to choose rho ? Several methods, implemented in R, provide automatic selection of penalty parameter. ‘GLASSOO’ package selects the penalty parameter using cross-validation.

# cross validation for lam
mod2 = GLASSO(X,nlam = 100, trace = "none")
plot(mod2)

mod3 = GLASSO(X,lam = exp(-1.773))
mod3$Omega

##           [,1]      [,2]      [,3]         [,4]      [,5]         [,6]
## [1,] 1.5541584 0.5416053 0.0000000  0.000000000 0.0000000  0.000000000
## [2,] 0.5416053 1.3271850 0.4403733  0.000000000 0.0000000  0.000000000
## [3,] 0.0000000 0.4403733 1.5639287  0.000000000 0.0000000  0.000000000
## [4,] 0.0000000 0.0000000 0.0000000  1.594146899 0.5342061 -0.005093197
## [5,] 0.0000000 0.0000000 0.0000000  0.534206141 1.2978359  0.526863102
## [6,] 0.0000000 0.0000000 0.0000000 -0.005093197 0.5268631  1.593558126

Multivariate Poisson-lognormal model

In this section: - we draw sample in a Poisson-lognormal law - we estimate partial correlations using ‘PLNmodels’ R package - we introduce covariable (abiotic environment) and see how it can introduce confounding effects

p <- 5
n <- 100

g <- igraph::make_star(p,mode = "undirected")
g = igraph::add_edges(g,c(2,3,3,4,4,5,5,2))
V(g)$name <- LETTERS[1:p]
E(g)$weight <- c(2,2,2,2,-0.5,-0.5,-0.5,-0.5)
A <- get.adjacency(g,attr = "weight")
diag(A) = 5
#same mean on each component
mu <- rep(1,p)
Z <- rmvnp(n, mu, Omega = as.matrix(A))

We represent the partial correlation network built from the precision matrix.

A = as.matrix(A)
p_cor_theo <- - A/(sqrt(diag(A) %*% t(diag(A))))
#removing self-loops
diag(p_cor_theo) = 0
#increasing edge width
p_cor_theo = 3*p_cor_theo
#representing the partial correlation network
g_pcor <- graph_from_adjacency_matrix(abs(p_cor_theo),weighted = TRUE,mode = "undirected")
g_pcor = asNetwork(g_pcor)
ggnet2(g_pcor,label = TRUE,edge.size = "weight",edge.color = c("red","red","red","red","blue","blue","blue","blue"))

#sampling in the Poisson
Y <- apply(Z,1:2,function(k) rpois(n = 1,lambda = exp(k)))
rownames(Y) <- paste0("s",1:n)
covariates <- data.frame(env = rep(1,n))
rownames(covariates) <- paste0("s",1:n)

#overdispersion: variance is greater than mean
mean(Y[,1])

## [1] 4.99

var(Y[,1])

## [1] 59.46455

#prepare data for PLNmodels
PLN_data <- prepare_data(Y,covariates)
#similar syntax than R GLM
network_models <- PLNnetwork(Y ~ 1 , data = PLN_data)

## 
##  Initialization...
##  Adjusting 30 PLN with sparse inverse covariance estimation
##  Joint optimization alternating gradient descent and graphical-lasso
##  sparsifying penalty = 1.160399 
    sparsifying penalty = 1.071826 
    sparsifying penalty = 0.9900148 
    sparsifying penalty = 0.9144478 
    sparsifying penalty = 0.8446488 
    sparsifying penalty = 0.7801775 
    sparsifying penalty = 0.7206272 
    sparsifying penalty = 0.6656224 
    sparsifying penalty = 0.614816 
    sparsifying penalty = 0.5678876 
    sparsifying penalty = 0.5245412 
    sparsifying penalty = 0.4845034 
    sparsifying penalty = 0.4475217 
    sparsifying penalty = 0.4133628 
    sparsifying penalty = 0.3818111 
    sparsifying penalty = 0.3526678 
    sparsifying penalty = 0.325749 
    sparsifying penalty = 0.3008848 
    sparsifying penalty = 0.2779186 
    sparsifying penalty = 0.2567053 
    sparsifying penalty = 0.2371112 
    sparsifying penalty = 0.2190127 
    sparsifying penalty = 0.2022956 
    sparsifying penalty = 0.1868546 
    sparsifying penalty = 0.1725921 
    sparsifying penalty = 0.1594183 
    sparsifying penalty = 0.1472501 
    sparsifying penalty = 0.1360106 
    sparsifying penalty = 0.125629 
    sparsifying penalty = 0.1160399 
##  Post-treatments
##  DONE!

#search for optimal penalty value
plot(network_models, "diagnostic")

As for the graphical lasso, we represent the coefficient path along the penalty parameter.

coefficient_path(network_models, corr = TRUE) %>% 
  ggplot(aes(x = Penalty, y = Coeff, group = Edge, colour = Edge)) + 
  geom_line(show.legend = FALSE) +  coord_trans(x="log10") + theme_bw()

PLNmodels allows selecting optimal penalty parameter according to an information criterion: the Bayesian Information Criterion (BIC)

model_BIC <- getBestModel(network_models, "BIC")  
#PLNmodels has a method to represent the estimated partial correlation network
model_BIC$plot_network()

Adding covariates in the PLN model

We now introduce a covariable that represents abiotic environment and influences the mean of the normal layer.

## Building mean matrix
#env 1: we assume that A and B are abundant
mu_1 <- c(5,5,1,1,1)
#env 1: we assume that C, D and E are abundant
mu_2 <- c(1,1,5,5,5)

#sampling Z values in env1 and env2
Z1 <- rmvnp(n/2, mu_1, Omega = as.matrix(A))
Z2 <- rmvnp(n/2, mu_2, Omega = as.matrix(A))
Z <- rbind(Z1,Z2)

Y <- apply(Z,1:2,function(k) rpois(n = 1,lambda = exp(k)))
PLN_data <- prepare_data(Y,covariates)

## Warning in common_samples(counts, covariates): There are no matching names in the count matrix and the covariates data.frame.
## Function will proceed assuming:
## - samples are in the same order;
## - samples are rows of the abundance matrix.

We now estimate a PLN model without taking into account the covariable.

#build the network model
network_models <- PLNnetwork(Y ~ 1 , data = as.data.frame(cbind(Y,covariates)))

## 
##  Initialization...
##  Adjusting 30 PLN with sparse inverse covariance estimation
##  Joint optimization alternating gradient descent and graphical-lasso
##  sparsifying penalty = 4.725373 
    sparsifying penalty = 4.364689 
    sparsifying penalty = 4.031536 
    sparsifying penalty = 3.723812 
    sparsifying penalty = 3.439577 
    sparsifying penalty = 3.177037 
    sparsifying penalty = 2.934537 
    sparsifying penalty = 2.710546 
    sparsifying penalty = 2.503652 
    sparsifying penalty = 2.312551 
    sparsifying penalty = 2.136036 
    sparsifying penalty = 1.972994 
    sparsifying penalty = 1.822397 
    sparsifying penalty = 1.683295 
    sparsifying penalty = 1.55481 
    sparsifying penalty = 1.436133 
    sparsifying penalty = 1.326514 
    sparsifying penalty = 1.225263 
    sparsifying penalty = 1.131739 
    sparsifying penalty = 1.045355 
    sparsifying penalty = 0.9655636 
    sparsifying penalty = 0.891863 
    sparsifying penalty = 0.8237878 
    sparsifying penalty = 0.7609088 
    sparsifying penalty = 0.7028293 
    sparsifying penalty = 0.6491829 
    sparsifying penalty = 0.5996314 
    sparsifying penalty = 0.553862 
    sparsifying penalty = 0.5115862 
    sparsifying penalty = 0.4725373 
##  Post-treatments
##  DONE!

plot(network_models, "diagnostic")

model_BIC <- getBestModel(network_models, "BIC")  
model_BIC$plot_network()

We detect a spurious positive association between C, D and E. This is due to the fact that C, D and E are abundant in env2. Similarly, we detect a spurious positive association between A and B. We now estimate a PLN model again but we take into account the covariable to hopefully remove spurious associations.

# Adding covariates
covariates <- data.frame(env = c(rep(1,n/2),rep(2,n/2)))

network_models <- PLNnetwork(Y ~ env , data = as.data.frame(cbind(Y,covariates)))

## 
##  Initialization...
##  Adjusting 30 PLN with sparse inverse covariance estimation
##  Joint optimization alternating gradient descent and graphical-lasso
##  sparsifying penalty = 0.9875057 
    sparsifying penalty = 0.9121302 
    sparsifying penalty = 0.8425081 
    sparsifying penalty = 0.7782002 
    sparsifying penalty = 0.7188008 
    sparsifying penalty = 0.6639354 
    sparsifying penalty = 0.6132577 
    sparsifying penalty = 0.5664483 
    sparsifying penalty = 0.5232118 
    sparsifying penalty = 0.4832755 
    sparsifying penalty = 0.4463875 
    sparsifying penalty = 0.4123151 
    sparsifying penalty = 0.3808434 
    sparsifying penalty = 0.351774 
    sparsifying penalty = 0.3249234 
    sparsifying penalty = 0.3001223 
    sparsifying penalty = 0.2772142 
    sparsifying penalty = 0.2560547 
    sparsifying penalty = 0.2365102 
    sparsifying penalty = 0.2184576 
    sparsifying penalty = 0.2017829 
    sparsifying penalty = 0.186381 
    sparsifying penalty = 0.1721547 
    sparsifying penalty = 0.1590143 
    sparsifying penalty = 0.1468769 
    sparsifying penalty = 0.1356659 
    sparsifying penalty = 0.1253106 
    sparsifying penalty = 0.1157458 
    sparsifying penalty = 0.106911 
    sparsifying penalty = 0.09875057 
##  Post-treatments
##  DONE!

plot(network_models, "diagnostic")

model_BIC <- getBestModel(network_models, "BIC")  
model_BIC$plot_network()

Training School BOTTOMS-UP

Aim

Gaussian graphical model

Building association network and sampling in multivariate normal distribution

Inference with the graphical lasso

Multivariate Poisson-lognormal model

Adding covariates in the PLN model