####
## This code accompany the paper
## Gollini I., and Murphy T.B. (2016)
## " Joint Modelling of Multiple Network Views "
## Journal of Computational and Graphical Statistics 25(1), 246-265
## arXiv:1301.3759
####
## author: Gollini I., <isabella.gollini@ucd.ie>
## date: last update of this code October 7th 2018
library(lvm4net)
D <- 2
########################
########################
####### Girls Datasets #######
########################
########################
library(RSiena) # load the datasets
Y1 <- s501
Y2 <- s502
Y3 <- s503
Nb <- nrow(Y1)
Y123 <- list(Y1 = Y1, Y2 = Y2, Y3 = Y3)
########################
########## LSM ##########
########################
modLSM1 <- lsm(Y1, D)
modLSM2 <- lsm(Y2, D)
modLSM3 <- lsm(Y3, D)
########################
########## LSJM ##########
########################
modLSJM123 <- lsjm(Y123, D) # It takes ~ 40 seconds
########################
#######match rotations #######
########################
Z123 <- modLSJM123$EZ
Z1 <- rotXtoY(modLSM1$lsmEZ, Z123)$Xrot
Z2 <- rotXtoY(modLSM2$lsmEZ, Z123)$Xrot
Z3 <- rotXtoY(modLSM3$lsmEZ, Z123)$Xrot
XYlimb <- range(Z123, Z1, Z2, Z3)
namesb <- paste('Network ', 1:3, sep ='')
colPl <- rainbow(Nb,alpha=.8)
########################
######## LSM PLOTS ########
########################
par(mfrow = c(1,3))
plotY(Y1, Ndata = 1, EZ = Z1, VZ = modLSM1$lsmVZ, main = namesb[1], xlim = XYlimb, ylim = XYlimb, colPl = colPl)
plotY(Y2, Ndata = 1, EZ = Z2, VZ = modLSM2$lsmVZ, main = namesb[2], xlim = XYlimb, ylim = XYlimb, colPl = colPl)
plotY(Y3, Ndata = 1, EZ = Z3, VZ = modLSM3$lsmVZ, main = namesb[3], xlim = XYlimb, ylim = XYlimb, colPl = colPl)
Zlsm <- list()
Zlsm[[1]] <- Z1
Zlsm[[2]] <- Z2
Zlsm[[3]] <- Z3
bpbLSM <- boxroc(Y123,
EZ = Zlsm,
xiT = c(modLSM1$xiT, modLSM2$xiT, modLSM3$xiT),
Lroc = 150,
ROC = TRUE,
BOXPLOT = TRUE,
labelsPlot = namesb
)
########################
######## LSJM PLOTS ########
########################
plot(modLSJM123, Y123, drawCB = TRUE, plotZtilde = TRUE, colPl = colPl)
plot(modLSJM123, Y123, drawCB = TRUE, colPl = colPl, main = 'Multiple networks')
bpbLSJM <- boxroc(Y123,
EZ = modLSJM123$lsmEZ,
xiT = modLSJM123$xiT,
Lroc = 150,
ROC = TRUE,
BOXPLOT = TRUE,
labelsPlot = namesb
)
########################
########################
####### PPI Datasets ########
########################
########################
data(PPInet) # provided in the package
Yg <- PPIgen
Yp <- PPIphy
Na <- nrow(Yg)
Ygp <- list()
Ygp[[1]] <- Yg
Ygp[[2]] <- Yp
########################
########## LSM ##########
########################
modLSMg <- lsm(Yg, D)
modLSMp <- lsm(Yp, D)
########################
########## LSJM ##########
########################
modLSJMgp <- lsjm(Ygp, D) # It takes ~ 45 seconds
########################
#######match rotations #######
########################
Zgp <- modLSJMgp$EZ
Zg <- rotXtoY(modLSMg$lsmEZ, Zgp)$Xrot
Zp <- rotXtoY(modLSMp$lsmEZ, Zgp)$Xrot
XYlima <- range(Zgp, Zg, Zp)
namesa <- c('Genetic Interactions', 'Physical Interactions')
colPl <- rainbow(Na,alpha=.8)
########################
######## LSM PLOTS ########
########################
par(mfrow = c(1,2))
plotY(Yg, Ndata = 1, EZ = Zg, VZ = modLSMg$lsmVZ, main = namesa[1], xlim = XYlima, ylim = XYlima, colPl = colPl)
plotY(Yp, Ndata = 1, EZ = Zp, VZ = modLSMp$lsmVZ, main = namesa[2], xlim = XYlima, ylim = XYlima, colPl = colPl)
Zlsmgp <- list()
Zlsmgp[[1]] <- Zg
Zlsmgp[[2]] <- Zp
bpaLSM <- boxroc(Ygp,
EZ = Zlsmgp,
xiT = c(modLSMg$xiT, modLSMp$xiT),
Lroc = 150,
ROC = TRUE,
BOXPLOT = TRUE,
labelsPlot = namesa
)
########################
######## LSJM PLOTS #######
########################
plot(modLSJMgp, Ygp, drawCB = TRUE, plotZtilde = TRUE, mainZtilde = namesa, colPl = colPl)
plot(modLSJMgp, Ygp, drawCB = TRUE, colPl = colPl)
bpaLSJM <- boxroc(Ygp,
EZ = modLSJMgp$lsmEZ,
xiT = modLSJMgp$xiT,
Lroc = 150,
ROC = TRUE,
BOXPLOT = TRUE,
labelsPlot = namesa
)
########################
########################
#### Supplementary Material ####
########################
########################
##########################
# Remove unconnected nodes from Yg #
##########################
cond <- seq(1:nrow(Yg))[rowSums(Yg) + colSums(Yg) == 0]
Yg2 <- Yg[-cond, -cond]
direct <- !all(Yg2 == t(Yg2))
nYg2 <-network(Yg2, directed = direct)
timings <- numeric(3)
names(timings) <- c('latentnet', 'VBLPCM', 'lvm4net')
########################
### Euclidean distance -- mcmc ###
########################
library(latentnet)
timings[1] <- system.time(euclmcmc <- ergmm(nYg2 ~ euclidean(d = D)))[[3]]
#######################
# Euclidean distance -- variational #
#######################
library(VBLPCM)
timings[2] <- system.time(euclvar<- vblpcmfit(vblpcmstart(nYg2, G = 1, d = D), maxiter = 300))[[3]]
##########################
# (Euclidean distance)^2 -- variational #
##########################
library(lvm4net)
timings[3] <- system.time(eucl2var <- lsm(Yg2, D) )[[3]]
########################
#######match rotations #######
########################
Z <- eucl2var$lsmEZ
Zm <- rotXtoY(euclmcmc$mkl$Z,Z)$Xrot
Zv <- rotXtoY(euclvar$V_z,Z)$Xrot
N <- nrow(Yg2)
Att<- rainbow(N,alpha=.8)
########################
######### LSM plots ########
########################
par(mfrow = c(1,3))
plotY(Yg2, Ndata = 1, EZ = Zm, main = 'latentnet', font.main = 1, cex.main = 1, colPl = Att)
plotY(Yg2, Ndata = 1, EZ = Zv, main = 'VBLPCM', font.main = 1, cex.main = 1, colPl = Att)
plotY(Yg2, Ndata = 1, EZ = Z, main = 'lvm4net', font.main = 1, cex.main = 1, colPl = Att)
########################
######### ROC plots ########
########################
Zmv <- list()
Zmv[[1]] <- Zm
Zmv[[2]] <- Zv
Zmv[[3]] <- Z
a <- boxroc(Yg2,
EZ = Zmv,
xiT = c(euclmcmc$mkl$beta, euclvar$V_xi_e, eucl2var$xiT),
Lroc = 150,
ROC = TRUE,
BOXPLOT = TRUE,
labelsPlot = c('latentnet', 'VBLPCM', 'lvm4net'),
powdist = c(1, 1, 2),
main =''
)
res <- cbind(a$AUC, timings)
colnames(res) <- c('AUC', 'Timings in Sec')
res
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.