search

GrossSBM / missSBM

An R package for adjusting Stochastic Block Models from networks data sampled under various missing data conditions
http://grosssbm.github.io/missSBM
GNU General Public License v3.0
12 stars 2 forks source link

Issues with missSBM::smooth #33

Closed TabouyT closed 5 years ago

TabouyT commented 5 years ago

The smooth function doesn't always improve the estimation, it makes it even worst and it shouldn't be ... Example with R code in JSS paper (.Rnw code)... See the ggplot attached with this issue

Rplot01.pdf

jchiquet commented 5 years ago

@TabouyT You are not using the last version of the package in the vignette, so it is hard to tell...

I experienced the same issue, so I decided to get back to the old version of ICL smoothing, which consists in choosing the best model according to the best ICL (not the best vBound).*

jchiquet commented 5 years ago

@TabouyT I changed a bit your code in JSS and recompile everything : it seems to work fine with the last version. ICL_both_example.pdf

TabouyT commented 5 years ago

Copy that ! It's moving too fast ... ;)

Demiperimetre commented 5 years ago

Comparison with blockmodels under a full sampling not consistent anymore... :{

Demiperimetre commented 5 years ago

R code : `r library(blockmodels) mod = BM_bernoulli("SBM_sym",beligerent_adjacency) mod$estimate() table(apply(mod$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[3]]$fittedSBM$memberships)

`

Demiperimetre commented 5 years ago

I give you my ticket that it is because of ICL instead of J

jchiquet commented 5 years ago

You can keep your ticket : it is juste because of the initialization ("hierarchical" in place au "spectral") : The following piece of code is consistent with blockmodels

collection_sbm_full <- 
  missSBM::estimate(
    sampledNet  = prepare_data(beligerent_adjacency), 
    vBlocks     = vBlocks, "node",
    clusterInit = "spectral"
  )
smooth(collection_sbm_full, "both")

I get

> table(apply(mod$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[3]]$fittedSBM$memberships)   
     1  2  3
  1  0  0 35
  2 42  0  0
  3  0 15  0

SO shall we always initialize with spectral clustering ? Hard to say...

jchiquet commented 5 years ago

I am closing this because it is not due to smoothing, but initialization

jchiquet commented 5 years ago

Well,

Il semblerait que le smoothing de l'ICL soit de retour.

Pour conclure, on fait ? Soit on utilise le J, mais on ne peut pas véritablement parler de smoothing de l'ICL dans ce cas. Soit on laisse le choix à l'utilisateur de lisser selon l'ICL ou la borne variationelle

D'autres idées ? Une opinion ? @Demiperimetre @TabouyT

Demiperimetre commented 5 years ago

...

TabouyT commented 5 years ago

Compliqué …

Le 11 avr. 2019 à 23:09, Pierre Barbillon notifications@github.com a écrit :

...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-482314930, or mute the thread https://github.com/notifications/unsubscribe-auth/AbRNCnInKq-T-sQDD0ugyLEOzwG2yQF0ks5vf6RsgaJpZM4ck-AY.

Demiperimetre commented 5 years ago

J'ai fait tourner le code du JSS. J'ai l'impression qu'il n'y a rien de catastrophique dans la mesure où le minimum est bien marqué ! Je serai d'avis à laisser le J s'il n'y a pas de trucs plus catastrophiques que ça quitte à faire une remarque dans le texte pour dire que c'est l'entropie qui peut fausser un lisssage... Je vais essayer d'explorer les cas où il y a des montagnes russes pour comprendre s'il n'y a pas un truc bizarre à ce niveau

ICLs

Demiperimetre commented 5 years ago

J'ai pushé sur le git document/papier-JSS les résultats sauvegardés et qq plots qui ont l'air d'indiquer que le problème de l'ICL vient de l'espérance du SBM.

Demiperimetre commented 5 years ago

it is called strange.R

Demiperimetre commented 5 years ago

a new smoothing proposition ??

jchiquet commented 5 years ago

merci, je regarde ça ce matin (12/04/2019)

jchiquet commented 5 years ago

Je pense avoir réconcilié la @Demiperimetre (vignette) et @TabouyT (JSS) en remplçant la fonction de split dans le forward par un sample sample (plutôt qu'un classif hiérarchique).

Demiperimetre commented 5 years ago

La classification en full observed est toujours cohérente avec BM !

TabouyT commented 5 years ago

Bon ... Le nouveau smoothing fait n'importe quoi sur les données de la blogosphère ... L'ICL est dégradé après 2 tours de lissage et à une mine affreuse par rapport à ne à lisser ...

Demiperimetre commented 5 years ago

Ça marche ? Sinon @TabouyT, tu peux extraire un Rdata avec la classif du smoothing pour vérifier ce qu'il se passe

Demiperimetre commented 5 years ago

on est repassé à ICL plutôt que J à nombre de classes fixées ?

jchiquet commented 5 years ago

@Demiperimetre

on est repassé à ICL plutôt que J à nombre de classes fixées ?

Oui... Faudrait tester voir si blockmodels fait toujours pareil. Je pense (j'espère) que les inconsistences étaient plutôt dues à la méthode de découpage des groupes (maintenant random, avant avec clustering hiérarchique).

Demiperimetre commented 5 years ago

et non c'est à nouveau le même problème. De plus, j'ai des warnings inquiétants de la fonction estimate Warning messages:nference for model with 6 blocks. 1: did not converge in 100 iterations 2: Quick-TRANSfer stage steps exceeded maximum (= 4150)

Est ce un autre problème ?

jchiquet commented 5 years ago

Warning messages:nference for model with 6 blocks. 1: did not converge in 100 iterations 2: Quick-TRANSfer stage steps exceeded maximum (= 4150)

Ce sont les k-means associés au spectral clustering qui n'ont pas convergé à l'initialisation

Demiperimetre commented 5 years ago

c'est étrange non alors que je suis sur les données complètes... On n'avait jamais eu ce message avant sur les mêmes données ?

TabouyT commented 5 years ago

Jamais de mon coté …

Le 18 avr. 2019 à 19:00, Pierre Barbillon notifications@github.com a écrit :

c'est étrange non alors que je suis sur les données complètes... On n'avait jamais eu ce message avant sur les mêmes données ?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-484592809, or mute the thread https://github.com/notifications/unsubscribe-auth/AG2E2CQCO66GTALTTBPCLM3PRCSL7ANCNFSM4HET4AMA.

jchiquet commented 5 years ago

Tu ne veux pas mettre ton bout de code ici stp ?

J'ai p-e changé le nombre d'itération max ou min dans le k-means

Demiperimetre commented 5 years ago

`r

data("war") library(igraph) B2 = as.matrix(get.adjacency(war$beligerent))

vBlocks <- 1:6 B2 = prepare_data(B2) collection_sbm_full <- missSBM::estimate(B2,vBlocks = vBlocks, "node") smooth(collection_sbm_full, "both") plot(collection_sbm_full$ICL, type = "l")

`

jchiquet commented 5 years ago

Il ne te plait pas l'ICL dans ce cas ? Moi je le trouve beau.

Demiperimetre commented 5 years ago

oui mais la classif est déconnante par rapport à blockmodels

Demiperimetre commented 5 years ago

`r library(blockmodels) mod2b = BM_bernoulli("SBM_sym",B2) mod2b$estimate()

table(apply(mod2b$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[3]]$fittedSBM$memberships)

`

TabouyT commented 5 years ago

De mon coté, enfin missSBM trouve le même nombre de classes que block models sans données manquantes.

Le 18 avr. 2019 à 19:07, Pierre Barbillon notifications@github.com a écrit :

`r library(blockmodels) mod2b = BM_bernoulli("SBM_sym",B2) mod2b$estimate()

table(apply(mod2b$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[3]]$fittedSBM$memberships)

`

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-484595019, or mute the thread https://github.com/notifications/unsubscribe-auth/AG2E2CTFSPXIZB6QXVHGUTLPRCTGTANCNFSM4HET4AMA.

Demiperimetre commented 5 years ago 
rm(list=ls())
library(igraph)
library(missSBM)
library(blockmodels)

data("war")
B2 = as.matrix(get.adjacency(war$beligerent))

vBlocks <- 1:6
B3 = prepare_data(B2)
collection_sbm_full <-  missSBM::estimate(B3,vBlocks = vBlocks, "node",clusterInit = "hierarchical")
smooth(collection_sbm_full, "both")
plot(collection_sbm_full$ICL, type = "l")

mod2b = BM_bernoulli("SBM_sym",B2)
mod2b$estimate()

# initialize from blockmodels
collection_sbmBM3 <-  missSBM::estimate(B3,vBlocks = 3, "node",clusterInit = apply(mod2b$memberships[[3]]$Z,1,which.max))
collection_sbmBM3$bestModel$fittedSBM$vICL

# comparing ICLs
mod2b$ICL[3]*-2
collection_sbm_full$bestModel$fittedSBM$vICL
collection_sbmBM3$bestModel$fittedSBM$vICL

# comparing clustering
Q=3
table(apply(mod2b$memberships[[Q]]$Z,1,which.max),collection_sbm_full$models[[Q]]$fittedSBM$memberships)
table(apply(mod2b$memberships[[Q]]$Z,1,which.max),collection_sbmBM3$bestModel$fittedSBM$memberships)
Demiperimetre commented 5 years ago

si on initialize avec blocmodel on est très proche mais l'icl est de 2 points moins bon. Je pense que ça peut être dû aux bornages des tau dans le VEM. À un moment on pouvait le bloquer dans les controles, j'ai l'impression que ce n'est plus possible

Demiperimetre commented 5 years ago 

# icl
iclalamano = function(pi,alpha,Z,adj)
{
  adj0 = 1-adj
  diag(adj0) = 0
  vA = 0
  for (i in 2:nrow(adj))
    for (j in 1:(i-1))
      for (q in 1:ncol(Z))
        for (l in 1:ncol(Z))
          vA = vA + Z[i,q]*Z[j,l]*(log(pi[q,l])*adj[i,j]+log(1-pi[q,l])*adj0[i,j])
  return(vA + sum(Z%*%log(alpha)))
}

pismooth = collection_sbm_full$bestModel$fittedSBM$connectParam
alphasmooth = collection_sbm_full$bestModel$fittedSBM$mixtureParam
Zsmooth = collection_sbm_full$bestModel$fittedSBM$blocks

pi3 = collection_sbmBM3$bestModel$fittedSBM$connectParam
alpha3 = collection_sbmBM3$bestModel$fittedSBM$mixtureParam
Z3 = collection_sbmBM3$bestModel$fittedSBM$blocks

piBM = mod2b$model_parameters[[3]]$pi
alphaBM = colMeans(mod2b$memberships[[3]]$Z)
ZBM = mod2b$memberships[[3]]$Z

-2*iclalamano(pismooth,alphasmooth,Zsmooth,B2) + log(nrow(B2)) * (ncol(Zsmooth)-1) + log(nrow(B2)*(nrow(B2)-1)/2)* (ncol(Zsmooth)+1)*ncol(Zsmooth)/2
-2*iclalamano(pi3,alpha3,Z3,B2) + log(nrow(B2)) * (ncol(Z3)-1) + log(nrow(B2)*(nrow(B2)-1)/2)* (ncol(Z3)+1)*ncol(Z3)/2
-2*iclalamano(piBM,alphaBM,ZBM,B2) + log(nrow(B2)) * (ncol(ZBM)-1) + log(nrow(B2)*(nrow(B2)-1)/2)* (ncol(ZBM)+1)*ncol(ZBM)/2

quand j'ajoute un seuil sur les valeurs de tau, les valeurs d'icl calculées à la main sont très différentes de celles données dans la classe

jchiquet commented 5 years ago

quand j'ajoute un seuil sur les valeurs de tau, les valeurs d'icl calculées à la main sont très différentes de celles données dans la classe

@Demiperimetre ce qui veut dire ? Tu vois une manière de changer le code qui solutionnerait le problème ?

Demiperimetre commented 5 years ago

non qu'il se passe qq chose d'étrange entre l'icl qu'il calcule et celui que je calcule à la main. Cela veut peut être dire qu'il n'utilise pas ce qu'il faut lorsqu'il compare les modèles ?

jchiquet commented 5 years ago

non qu'il se passe qq chose d'étrange entre l'icl qu'il calcule et celui que je calcule à la main. Cela veut peut être dire qu'il n'utilise pas ce qu'il faut lorsqu'il compare les modèles ?

C'est qui "il" ?

Demiperimetre commented 5 years ago

l'algo de smoothing

Demiperimetre commented 5 years ago 

vBlocks <- 1:6
B3 = prepare_data(B2)
collection_sbm_full <-  missSBM::estimate(B3,vBlocks = vBlocks, "node",clusterInit = "spectral")
collection_sbm_full$ICL
collection_sbm_full$bestModel$fittedSBM$vICL
pibrut = collection_sbm_full$bestModel$fittedSBM$connectParam
alphabrut = collection_sbm_full$bestModel$fittedSBM$mixtureParam
Zbrut = collection_sbm_full$bestModel$fittedSBM$blocks
-2*iclalamano(pibrut,alphabrut,Zbrut,B2) + log(nrow(B2)) * (ncol(Zbrut)-1) + log(nrow(B2)*(nrow(B2)-1)/2) * (ncol(Zbrut)+1)*ncol(Zbrut)/2

toujours le même problème. Le truc que je n'arrive pas à comprendre dans le code c'est à quel moment, on calcule les Vexpec, on dirait qu'il n'utilise pas les valeurs de pi, alpha, Z qu'il renvoie...

jchiquet commented 5 years ago

Bon,

je ne sais pas ce qu'il s'était passé, mais manifestement une des corrections que j'avais faite dans le calcul de la borne variationnelle n'avait pas été pushé. Sans doute un conflit mal résolu.

Je trouve désormais pareil que "ICL à la Mano".

Et j'ai viré le code de update_block, l'ancien apparemment faisait le job.

jchiquet commented 5 years ago

@TabouyT @Demiperimetre Merci de faire retourner tout avec le dernier code pour qu'on en finisse...

jchiquet commented 5 years ago

Bon,

Le calcule de l'ICL est désormais correct. On ne trouve pas comme blockmodel à nombre fixe de classe pour l'exemple de guerre mais il s'agit sans doute d'une question d'initialisation. On pourrait envisager d'utiliser le kmeans de blockmodels pour initialiser, c'est p-e ça le soucis.

En tout cas je reviens au smoothing de l'ICL : j'ai essayé toutes les configurations (tau, ICL, vBound) et celle qui me paraît la plus juste vis à vis de l'exemple de JSS et celle sans le filtre sur les tau et avec le smoothing de l'ICL.

Et dans ce cas au moins le code suivant est cohérent entre blockmodel et nous. Je regarderai cette histoire d'initialisation avec les kmeans de blockmodel plutôt ue ceux de R, c'est p-e ça le soucis

data("war")
library(igraph)
B2 = as.matrix(get.adjacency(war$beligerent))

library(blockmodels)
mod2b = BM_bernoulli("SBM_sym",B2)
mod2b$estimate()

B3 = prepare_data(B2)
# collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node")
collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node", clusterInit = list(apply(mod2b$memberships[[3]]$Z,1,which.max)))

print(table(apply(mod2b$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[1]]$fittedSBM$memberships))
TabouyT commented 5 years ago

Les dernières modifications marche bien de mon coté, j’ai un exemple satisfaisant.

T.

Le 27 avr. 2019 à 10:34, Julien Chiquet notifications@github.com a écrit :

Bon,

Le calcule de l'ICL est désormais correct. On ne trouve pas comme blockmodel à nombre fixe de classe pour l'exemple de guerre mais il s'agit sans doute d'une question d'initialisation. On pourrait envisager d'utiliser le kmeans de blockmodels pour initialiser, c'est p-e ça le soucis.

En tout cas je reviens au smoothing de l'ICL. Et au moins le code suivant est cohérent

data("war") library(igraph) B2 = as.matrix(get.adjacency(war$beligerent))

library(blockmodels) mod2b = BM_bernoulli("SBM_sym",B2) mod2b$estimate()

B3 = prepare_data(B2)

collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node")

collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node", clusterInit = list(apply(mod2b$memberships[[3]]$Z,1,which.max)))

print(table(apply(mod2b$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[1]]$fittedSBM$memberships)) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-487267382, or mute the thread https://github.com/notifications/unsubscribe-auth/AG2E2CTQFND5CA66MVMSDE3PSQFY7ANCNFSM4HET4AMA.

Demiperimetre commented 5 years ago

De mon côté, j'obtiens les mêmes résultats. Lorsqu'on lance le smoothing de missSBM on tombe sur un ICL plus faible de 2 points pour une classification totalement différente de blockmodels. Ce qu'il y a de bizarre est que blockmodel donne un ICL de 1267, missSBM avec l'initialisation SBM 1273 et en lissant à partir d'un spectral 1271. Je ne comprends pas d'où vient cette différence de 4-5 points, j'avais supposé que c'était dû au seuillage sur les taus.

Mais tant pis, on a déjà passé assez de temps là-dessus. Si ça marche pour l'exemple du JSS c'est bien. Pour la vignette je vais analyser avec la classif de missSBM.

+++ P.

Le sam. 27 avr. 2019 à 10:26, Tabouy Timothée notifications@github.com a écrit :

Les dernières modifications marche bien de mon coté, j’ai un exemple satisfaisant.

T.

Le 27 avr. 2019 à 10:34, Julien Chiquet notifications@github.com a écrit :

Bon,

Le calcule de l'ICL est désormais correct. On ne trouve pas comme blockmodel à nombre fixe de classe pour l'exemple de guerre mais il s'agit sans doute d'une question d'initialisation. On pourrait envisager d'utiliser le kmeans de blockmodels pour initialiser, c'est p-e ça le soucis.

En tout cas je reviens au smoothing de l'ICL. Et au moins le code suivant est cohérent

data("war") library(igraph) B2 = as.matrix(get.adjacency(war$beligerent))

library(blockmodels) mod2b = BM_bernoulli("SBM_sym",B2) mod2b$estimate()

B3 = prepare_data(B2)

collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node")

collection_sbm_full <- missSBM::estimate(B3, vBlocks = 3, "node", clusterInit = list(apply(mod2b$memberships[[3]]$Z,1,which.max)))

print(table(apply(mod2b$memberships[[3]]$Z,1,which.max),collection_sbm_full$models[[1]]$fittedSBM$memberships)) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/jchiquet/missSBM/issues/33#issuecomment-487267382>, or mute the thread < https://github.com/notifications/unsubscribe-auth/AG2E2CTQFND5CA66MVMSDE3PSQFY7ANCNFSM4HET4AMA .

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-487290424, or mute the thread https://github.com/notifications/unsubscribe-auth/AHP4KVT6Z3FE5KQPVFIPNL3PSRPB5ANCNFSM4HET4AMA .

-- Pierre Barbillon https://www6.inra.fr/mia-paris/Equipes/Membres/Pierre-Barbillon

jchiquet commented 5 years ago

Ok, mais gardons cette issue ouverte... Ça serait bien qu'un jour on fasse exactement comme blockmodels su cet exemple... P-e un problème d'optim, de point fixe... Je ne sais pas.

TabouyT commented 5 years ago

Ok, je propose quand même d'arrêter de se triturer l'esprit pour peu, c'est sans doute comme dit julien des petits réglages qui un coup donneront l'avantage à l'un puis un autre coup à l'autre. On avance.

T.

Le lun. 29 avr. 2019 à 07:25, Julien Chiquet notifications@github.com a écrit :

Ok, mais gardons cette issue ouverte... Ça serait bien qu'un jour on fasse exactement comme blockmodels su cet exemple... P-e un problème d'optim, de point fixe... Je ne sais pas.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jchiquet/missSBM/issues/33#issuecomment-487454982, or mute the thread https://github.com/notifications/unsubscribe-auth/AG2E2CWPNSGVOXC5S7GV7GDPS2BFBANCNFSM4HET4AMA .

jchiquet commented 5 years ago

Ok, donc ça serait cool qu'on ait une proposition complète pour la vignette à mon retour de vacances que je puisse soumettre au CRAN ...

Demiperimetre commented 5 years ago

:s