Poisson log-normal model for fish species: network inference
S. Robin
nov 2021
Data
- \(n =89\) stations (sites) in the Barents’ sea
- \(m = 30\) fish species
- \(d = 4\) environmental covariates: latitude, longitude, depth, temperature
rm(list=ls()); set.seed(1)
library(PLNmodels)## This is package 'PLNmodels' version 1.2.2## Use future::plan(multicore/multisession) to speed up PLNPCA/PLNmixture/stability_selection.data(barents)
X <- as.matrix(cbind(barents$Latitude, barents$Longitude, barents$Depth, barents$Temperature))
colnames(X) <- c('Latitude', 'Longitude', 'Depth', 'Temperature')
Y <- as.matrix(barents$Abundance)Scaling covariates (for numerical reasons)
X <- scale(X)
print(head(X))## Latitude Longitude Depth Temperature
## [1,] -1.645660 -1.271961723 0.3006760 1.976666
## [2,] -1.483765 -0.993522216 0.8001587 1.785028
## [3,] -1.277716 -0.721765258 -0.5317952 1.497572
## [4,] -1.520559 -0.503468685 -0.3804368 1.689210
## [5,] -1.336587 -0.004505089 0.8304304 1.401753
## [6,] -1.366023 0.213791485 0.2249968 1.689210d <- ncol(X)
Y <- Y[, which(colMeans(Y)>1)]
n <- nrow(Y); m <- ncol(Y)
print(m)## [1] 18print(head(Y[, 1:10]))## Re_hi An_de Hi_pl Me_ae Ra_ra Mi_po Ar_at Ma_vi Bo_sa Cl_ha
## 1 0 0 31 108 0 325 0 2 0 4
## 2 0 0 4 110 0 349 0 0 0 0
## 3 0 0 27 788 0 6 0 0 0 0
## 4 0 0 13 295 0 2 0 3 0 0
## 5 0 0 23 13 2 240 0 16 0 3
## 6 1 0 20 97 0 0 0 0 0 15Removing too rare species
X <- scale(X)
print(head(X))## Latitude Longitude Depth Temperature
## [1,] -1.645660 -1.271961723 0.3006760 1.976666
## [2,] -1.483765 -0.993522216 0.8001587 1.785028
## [3,] -1.277716 -0.721765258 -0.5317952 1.497572
## [4,] -1.520559 -0.503468685 -0.3804368 1.689210
## [5,] -1.336587 -0.004505089 0.8304304 1.401753
## [6,] -1.366023 0.213791485 0.2249968 1.689210d <- ncol(X)
Y <- Y[, which(colMeans(Y)>1)]
n <- nrow(Y); m <- ncol(Y)
print(m)## [1] 18print(head(Y[, 1:10]))## Re_hi An_de Hi_pl Me_ae Ra_ra Mi_po Ar_at Ma_vi Bo_sa Cl_ha
## 1 0 0 31 108 0 325 0 2 0 4
## 2 0 0 4 110 0 349 0 0 0 0
## 3 0 0 27 788 0 6 0 0 0 0
## 4 0 0 13 295 0 2 0 3 0 0
## 5 0 0 23 13 2 240 0 16 0 3
## 6 1 0 20 97 0 0 0 0 0 15PLN model
Diagonal covariance matrix
No interaction: \(\Sigma\) (and therefore \(\Omega\)) diagonal
plnDiag <- PLN(Y ~ X, control=PLN_param(trace=0, covariance='diagonal'))
plnDiagCrit <- c(unlist(plnDiag$criteria), c(n_edges=0, EBIC=plnDiag$criteria$BIC, pen_loglik=plnDiag$criteria$loglik, density=0))
print(plnDiagCrit)## nb_param loglik BIC ICL n_edges EBIC pen_loglik
## 108.000 -3987.140 -4229.526 -5448.262 0.000 -4229.526 -3987.140
## density
## 0.000PlotCov(plnDiag)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
Full covariance matrix
All interactions: \(\Sigma\) (and possibly \(\Omega\)) full
plnFull <- PLN(Y ~ X, control=PLN_param(trace=0, covariance='full'))
plnFullCrit <- c(unlist(plnFull$criteria), c(n_edges=m*(m-1)/2, EBIC=plnFull$criteria$BIC, pen_loglik=plnFull$criteria$loglik, density=1))
print(plnFullCrit)## nb_param loglik BIC ICL n_edges EBIC pen_loglik
## 261.000 -3859.482 -4445.249 -5482.684 153.000 -4445.249 -3859.482
## density
## 1.000PlotCov(plnFull)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
Abiotic effects
cluster <- kmeans(t(plnFull$model_par$B[-1, ]), centers=4, iter.max=1e2, nstart=1e2)$cluster
par(mfrow=c(1, 2))
image.plot(1:d, 1:m, plnFull$model_par$B[-1, ], main='abiotic effects', xlab='covariates', ylab='species')
text(1:d, rep(2, d), labels=substr(colnames(X), 1, 4), col=1)
image.plot(1:d, 1:m, plnFull$model_par$B[-1, order(cluster)], main='species clustering', xlab='covariates', ylab='species')
abline(h=.5+cumsum(table(cluster)), col=1, lwd=2)
text(1:d, rep(2, d), labels=substr(colnames(X), 1, 4), col=1)
Reordering the covariance matrix
PlotCov(plnFull, cluster=cluster)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
Network inference
Few interactions: \(\Omega\) sparse
plnNet <- PLNnetwork(Y ~ X, control=PLNnetwork_param(trace=0))
plot(plnNet)
Extending the grid for the penalty parameter \(\lambda\)
lambda <- max(plnNet$penalties)^seq(-2, 1, length.out=50)
plnNet <- PLNnetwork(Y ~ X, penalties=lambda, control=PLNnetwork_param(trace=0))
plot(plnNet)
grid
A first network
netBic <- plnNet$getBestModel(crit='BIC')
print(netBic$criteria)## nb_param loglik BIC ICL n_edges EBIC pen_loglik
## 1 119 -4048.801 -4315.874 -5933.665 11 -4330.353 -4051.978
## density
## 1 0.06790123PlotCov(netBic)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
par(mfrow=c(1, 2))
plot(netBic)## Warning: vertex attribute frame.color contains NAs. Replacing with default
## value blackPlotNet(netBic, color=cluster)
Stability selection
Resampling procedure to assess edge robustness
stability_selection(plnNet)##
## Stability Selection for PLNnetwork:
## subsampling: ++++++++++++++++++++netStab <- plnNet$getBestModel(crit='StARS')
print(netStab$criteria)## nb_param loglik BIC ICL n_edges EBIC pen_loglik
## 1 116 -4078.822 -4339.163 -6000.742 8 -4350.967 -4082.154
## density
## 1 0.04938272PlotCov(netStab)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
par(mfrow=c(1, 2))
plot(netStab)## Warning: vertex attribute frame.color contains NAs. Replacing with default
## value blackPlotNet(netStab, color=cluster)
Model comparison
print(rbind(plnDiag = plnDiagCrit,
plnFull = plnFullCrit,
netBic = netBic$criteria,
netStab = netStab$criteria))## nb_param loglik BIC ICL n_edges EBIC pen_loglik
## plnDiag 108 -3987.140 -4229.526 -5448.262 0 -4229.526 -3987.140
## plnFull 261 -3859.482 -4445.249 -5482.684 153 -4445.249 -3859.482
## netBic 119 -4048.801 -4315.874 -5933.665 11 -4330.353 -4051.978
## netStab 116 -4078.822 -4339.163 -6000.742 8 -4350.967 -4082.154
## density
## plnDiag 0.00000000
## plnFull 1.00000000
## netBic 0.06790123
## netStab 0.04938272Comparing estimated abiotic effects
abiotic = cbind(plnDiag = as.vector(plnDiag$model_par$B[-1, ]),
plnFull = as.vector(plnFull$model_par$B[-1, ]),
netBic = as.vector(netBic$model_par$B[-1, ]),
netStab = as.vector(netStab$model_par$B[-1, ]))
plot(as.data.frame(abiotic), pch=20)
Not accounting for environmental covariates
PLN model
plnNoCovar <- PLN(Y ~ 1, control=PLN_param(trace=0))
PlotCov(plnNoCovar)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
PlotCov(plnNoCovar, cluster=cluster)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
plnNoCovarCrit <- c(unlist(plnNoCovar$criteria),
c(n_edges=m*(m-1)/2, EBIC=plnNoCovar$criteria$BIC,
pen_loglik=plnNoCovar$criteria$loglik, density=0))Network inference
lambda <- max(plnNet$penalties)^seq(-4, 1, length.out=50)
plnNetNoCovar <- PLNnetwork(Y ~ 1, penalties=lambda, control=PLNnetwork_param(trace=0))
plot(plnNetNoCovar)
netNoCovarBic <- plnNetNoCovar$getBestModel(crit='BIC')
PlotCov(netNoCovarBic, cluster=cluster)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
stability_selection(plnNetNoCovar)##
## Stability Selection for PLNnetwork:
## subsampling: ++++++++++++++++++++netNoCovarStab <- plnNetNoCovar$getBestModel(crit='StARS')
PlotCov(netNoCovarStab, cluster=cluster)## Warning in plot.window(...): "color.scale" n'est pas un paramètre graphique## Warning in plot.xy(xy, type, ...): "color.scale" n'est pas un paramètre
## graphique## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique
## Warning in axis(side = side, at = at, labels = labels, ...): "color.scale"
## n'est pas un paramètre graphique## Warning in box(...): "color.scale" n'est pas un paramètre graphique## Warning in title(...): "color.scale" n'est pas un paramètre graphique
Network comparison
print(rbind(netBic = netBic$criteria,
netNoCovarBic = netNoCovarBic$criteria,
netStab = netStab$criteria,
netNoCovarStab = netNoCovarStab$criteria))## nb_param loglik BIC ICL n_edges EBIC
## netBic 119 -4048.801 -4315.874 -5933.665 11 -4330.353
## netNoCovarBic 158 -4139.988 -4494.591 -5959.710 122 -4508.402
## netStab 116 -4078.822 -4339.163 -6000.742 8 -4350.967
## netNoCovarStab 187 -4108.460 -4528.147 -5871.380 151 -4529.141
## pen_loglik density
## netBic -4051.978 0.06790123
## netNoCovarBic -4141.198 0.75308642
## netStab -4082.154 0.04938272
## netNoCovarStab -4108.544 0.93209877