Training School BOTTOMS-UP
Marc Ohlmann, Laboratoire d’Écologie Alpine
2022-09-30
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"
## 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
## 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
## [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()
