Timothée Tabouy

Pierre Barbillon

Julien Chiquet

Variational Inference for Stochastic Block Models from Sampled Data

This paper deals with non-observed dyads during the sampling of a network and consecutive issues in the inference of the Stochastic Block Model (SBM). We review sampling designs and recover Missing At Random (MAR) and Not Missing At Random (NMAR) conditions for the SBM. We introduce variants of the variational EM algorithm for inferring the SBM under various sampling designs (MAR and NMAR) all available as an R package. Model selection criteria based on Integrated Classification Likelihood are derived for selecting both the number of blocks and the sampling design. We investigate the accuracy and the range of applicability of these algorithms with simulations. We explore two real-world networks from ethnology (seed circulation network) and biology (protein-protein interaction network), where the interpretations considerably depends on the sampling designs considered. Stochastic Block Model • Variational inference • Missing data • Sampled network

Introduction

Networks arise in many fields of application for providing an intuitive way to represent interactions between entities. In this paper, a network is composed by a fixed set of nodes, and an interaction between a pair of nodes (dyad) is called an edge. We consider undirected binary networks with no loop, which can be represented by symmetric adjacency matrices filled with zeros and ones.

Various statistical models exist for depicting the probability distribution of the adjacency matrix (see, e.g. [START_REF] Goldenberg | A survey of statistical network models[END_REF][START_REF] Snijders | Statistical models for social networks[END_REF], for a survey). A highly desirable feature is their capability to describe the heterogeneity of real-world networks. In this perspective, the family of models endowed with a latent structure (reviewed in [START_REF] Matias | Modeling heterogeneity in random graphs through latent space models: a selective review[END_REF] offers a natural way to introduce heterogeneity. Within this family the Stochastic Block Model (in short SBM, see Frank and[START_REF] Frank | Cluster inference by using transitivity indices in empirical graphs[END_REF][START_REF] Holland | Stochastic blockmodels: First steps[END_REF] describes a broad variety of network topologies by positing a latent structure (or a clustering) on the nodes, then making the probability distribution of the adjacency matrix dependent on this latent structure. In order to estimate SBMs, Bayesian approaches were first developed [START_REF] Snijders | Statistical models for social networks[END_REF][START_REF] Snijders | Estimation and prediction for stochastic blockmodels for graphs with latent block structure[END_REF]Nowicki and[START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] prior to variational approaches [START_REF] Daudin | A mixture model for random graphs[END_REF][START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF]. On the theoretical side, [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] study the conditions for identifiability and the consistency of the variational estimators; [START_REF] Bickel | Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels[END_REF] prove their asymptotic normality. Several generalizations are possible such as weighted or directed variants [START_REF] Mariadassou | Uncovering latent structure in valued graphs: A variational approach[END_REF], mixed-membership and overlapping SBM [START_REF] Airoldi | Mixed membership stochastic blockmodels[END_REF][START_REF] Latouche | Overlapping stochastic block models with application to the french political blogosphere[END_REF], degree-corrected SBM [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF], dynamic SBM [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], or multiplex SBM [START_REF] Barbillon | Stochastic block models for multiplex networks: an application to networks of researchers[END_REF].

This paper deals with inference in the SBM when the network is not fully observed. We consider cases where all the nodes are observed but information regarding the presence/absence of an edge is missing for some dyads. In other words the adjacency matrix contains missing values, a situation often met with real-world networks. For instance in social sciences, network data consists in interactions between individuals: the set of individuals is fixed, possibly known from a census. Information about the presence/absence of an edge is only available when at least one of the two individuals is available for an interview, otherwise it is missing. See [START_REF] Thompson | Model-based estimation with link-tracing sampling designs[END_REF], [START_REF] Thompson | Adaptive Sampling[END_REF], [START_REF] Kolaczyk | Statistical analysis of network data, methods and models[END_REF] and [START_REF] Handcock | Modeling social networks from sampled data[END_REF] for a review of network sampling techniques. Even though some papers deal with SBM inference under missing data condition [START_REF] Aicher | Learning latent block structure in weighted networks[END_REF][START_REF] Vinayak | Graph clustering with missing data: Convex algorithms and analysis[END_REF], the sampling mechanism responsible for the missing values is overlooked in the inference, contrary to the approach developed in our paper.

Our contributions. A typology of sampling designs is introduced in Section 2.2. We adapt the theory developed in [START_REF] Rubin | Inference and missing data[END_REF] and [START_REF] Little | Statistical analysis with missing data[END_REF] to the SBM by splitting the sampling designs into the three usual classes of missing data: i) Missing Completely At Random (MCAR), where the sampling does not depend on the data, neither on the observed nor on the unobserved part of the network.

ii) Missing At Random (MAR), where the probability of being sampled is independent on the value of the missing data. For network data, the sampling does not depend on the presence/absence of an edge of an unobserved (or missing) dyad. MCAR is a particular case of MAR.

iii) Not Missing At Random (NMAR), where the sampling scheme is guided by unobserved dyads in some way.

Section 2.3 introduces several examples of sampling designs (MAR and NMAR) for which we derive conditions for identifiability of the SBM parameters. Estimation of the SBM in the MAR cases can be handled with the Variational EM (VEM) of [START_REF] Daudin | A mixture model for random graphs[END_REF] by conducting the inference only on the observed part of the network (Section 3.1). NMAR is more difficult to deal with as the sampling design must be taken into account in the inference. We introduce in Section 3.2 a general variational algorithm [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF] to deal with NMAR cases when the sampling design relies on a probability distribution which is explicitly known1 . Our variational approach is based on a double mean-field approximation applied to the latent distribution of the clustering and to the distribution of the missing dyads. We implement VEM algorithms that produce unbiased estimators for three natural NMAR sampling designs: a dyad-centered strategy, a node-centered strategy, and a block-centered strategy. We also derive an Integrated Classification Likelihood criterion (ICL, [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] for selecting the number of blocks. Although it is not possible to distinguish whether the sampling is MAR or NMAR [START_REF] Molenberghs | Every missing not at random model got a missing at random counterpart with equal fit[END_REF], the ICL can also be used to select which sampling design is the best fit for the data.

In Section 4.2 we show the good performance of our VEM algorithms on simulations for both MAR (Section 4.1) and NMAR conditions. Finally we investigate two very different real-world networks with missing values, namely a Kenyan seed exchange network (Section 5.1), and a protein-protein interaction (PPI) network (Section 5.2).

Related works.

In the few papers dealing with missing data for networks, the sampling design is rarely discussed. Even if not explicitly stated they all assume MAR conditions. [START_REF] Aicher | Learning latent block structure in weighted networks[END_REF] propose a weighted SBM modeling simultaneously the presence/absence of an edge and its weight. Missing data are handled by dropping the corresponding terms in the likelihood and the inference is conducted by a variational algorithm. In [START_REF] Vincent | Estimating the size and distribution of networked populations with snowball sampling[END_REF] a Bayesian augmentation procedure is introduced to estimate simultaneously the size of the population and the clustering when the sampling design is a one-wave snowball. Apart from the SBM, the exponential random graph model has been studied in the MAR setting in [START_REF] Handcock | Modeling social networks from sampled data[END_REF].

The matrix completion literature brings additional insights since SBM inference can be seen as a low-rank matrix estimation. [START_REF] Vinayak | Graph clustering with missing data: Convex algorithms and analysis[END_REF] introduce a convex program for the matrix completion problem where the underlying matrix has a simple affiliation structure defined via an SBM. The entries are sampled independently with the same probability, corresponding to a MAR case. In [START_REF] Davenport | 1-bit matrix completion[END_REF] the case of noisy 1-bit observations is studied and a likelihood-based strategy is developed with theoretical justifications ensuring good matrix completion. [START_REF] Chatterjee | Matrix estimation by universal singular value thresholding[END_REF] proves strong results for large matrices with noisy entries estimation, by means of a universal singular value thresholding.

Another related question is when the status of some dyads (absence/presence) is not clear in errorfully observed graph. Such uncertainties can be taken into account [START_REF] Priebe | Statistical inference on errorfully observed graphs[END_REF][START_REF] Balachandran | On the propagation of low-rate measurement error to subgraph counts in large networks[END_REF]. The latter reference studies the error propagation made by using estimators computed on observed sub-graphs, in order to estimate the number of existing edges in the real underlying graph.

Statistical framework

Stochastic Block Model

In an SBM, nodes from a set N {1, . . . , n} are distributed among a set Q {1, . . . , Q} of hidden blocks that model the latent structure of the graph. The blocks are described by the latent random vectors Z i q = (Z i1 , . . . , Z iQ ) i∈N with multinomial distribution M(1, α = (α 1 , . . . , α Q )). The probability of an edge between any dyad in D N × N only depends on the blocks the two nodes belong to. Hence, the presence of an edge between i and j, indicated by the binary variable Y ij , is independent on the other edges conditionally on the latent blocks:

Y ij | Z iq = 1, Z j = 1 ∼ ind B(π q ), ∀(i, j) ∈ D, ∀(q, ) ∈ Q × Q,
where B stands for the Bernoulli distribution. In the following, π = (π q ) (q, )∈Q×Q is the

Q × Q matrix of connectivity probabilities, Y = (Y ij ) (i,j
)∈D is the n × n adjacency matrix of the random graph, Z = (Z iq ) i∈N ,q∈Q is the n × Q matrix of the latent blocks and θ = (α, π) are the unknown parameters. In the undirected binary case, Y ij = Y ji for all (i, j) ∈ D and Y ii = 0 for all i ∈ N . Similarly, π q = π q for all (q, ) ∈ Q × Q.

Sampled data in the SBM framework

The sampled data is an n × n matrix with entries in {0, 1, NA}. It corresponds to the adjacency matrix Y where unobserved dyads have been replaced by NA's. More formally, let R be the n × n sampling matrix recording the data sampled during this process, such that R ij = 1 if Y ij is observed and 0 otherwise; also define

D o = {(i, j) : R ij = 1}, D m = {(i, j) : R ij = 0}, Y o = {Y ij : (i, j) ∈ D o } and Y m = {Y ij : (i, j) ∈ D m } to
denote the sets of variables respectively associated with the observed and missing data.

The number of nodes n is assumed to be known. The sampling design is the description of the stochastic process that generates R. It is assumed that the network exists before the sampling design acts upon it. Moreover, the sampling design is fully characterized by the conditional distribution p ψ (R|Y ), the parameters of which are such that ψ and θ live in a product space Θ × Ψ. Hence the joint probability density function of the observed data satisfies

p θ,ψ (Y o , R) = p θ (Y o , Y m , Z)p ψ (R|Y o , Y m , Z)dY m dZ. (1)
Simplifications may occur in (1) depending on the sampling design, leading to the three usual types of missingness (MCAR, MAR and NMAR). This typology depends on the relations between the adjacency matrix Y , the latent structure Z and the sampling R, so that the missingness is characterized by four directed acyclic graphs displayed in Figure 1.

Z Y R Z Y R (a) (b) Z Y R Z Y R (c) (d)
Figure 1: DAGs of relationships between Y, Z and R in the framework of missing data for SBM. DAG where R is a parent node are not reviewed since the network exists before the sampling design acts upon it. The systematic edge between Z and Y is due to the definition of the SBM. Note that the DAG (b) may correspond to MAR or NMAR samplings.

On the basis of these DAGs, the sampling design is

MCAR if R |= (Y m , Z, Y o ), MAR if R |= (Y m , Z) | Y o
, and NMAR otherwise. We derive Proposition 1 from these definitions.

Proposition 1. If the sampling is MCAR or MAR then i) arg max θ p θ,ψ (Y o , R) = arg max θ p θ (Y o )
for any ψ such that p θ,ψ (Y o , R) = 0 and ii) the sampling design necessary satisfies DAG (a) or (b).

Proof.

To prove i), if R satisfies MAR conditions, then p ψ (R|Y o , Y m , Z) = p ψ (R|Y o ).
Moreover, θ and ψ lie in a product space so that (1) factorizes into

p θ,ψ (Y o , R) = p θ (Y o )p ψ (R|Y o ).
This corresponds to the ignorability condition of [START_REF] Rubin | Inference and missing data[END_REF] and [START_REF] Handcock | Modeling social networks from sampled data[END_REF]. The proof of ii) is postponed to the supplementary materials.

Sampling design examples

MAR examples

Definition 1 (Random-dyad sampling). Each dyad (i, j) ∈ D has the same probability P(R ij = 1) = ρ to be observed independently of the others.

This design is trivially MCAR because each dyad is sampled with the same probability ρ which does not depend on Y .

Definition 2 (Star and snowball sampling). The star sampling consists in selecting uniformly a set of nodes, then observing corresponding rows of matrix Y . Snowball sampling is initialized by a star sampling which gives a first "wave" of nodes. The second wave is composed by the neighbors of the first. Successive waves can then be obtained. The final set of observed dyads corresponds to all dyads involving at least one of these nodes.

These two designs are node-centered and MAR. Indeed, selecting nodes independently in star sampling or in the first wave of snowball sampling corresponds to MCAR sampling. Successive waves are then MAR since they are built on the basis of the previously observed part of Y . Expressions of the corresponding distributions p ψ (R|Y o ) are given in [START_REF] Handcock | Modeling social networks from sampled data[END_REF].

Identifiability of random-dyad and star sampling designs. Since random-dyad and star samplings are MCAR, the identifiability is assessed in two steps by proving the identifiability of, first, the sampling parameter ψ = ρ and second, the SBM parameters θ = (α, π) given ρ. Our proofs, postponed to the supplementary materials, follow [START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF] who established the identifiability of the SBM without missing data.

Proposition 2. The sampling parameter ρ > 0 of random-dyad (resp. star) sampling is identifiable w.r.t. the sampling distribution.

Theorem 1. Let n ≥ 2Q and assume that for any 1 ≤ q ≤ Q, ρ > 0, α q > 0 and that the coordinates of πα are pairwise distinct. Then, under random-dyad (resp. star) sampling, SBM parameters are identifiable w.r.t. the distribution of the observed part of the SBM up to label switching.

NMAR examples

Definition 3 (Double standard sampling). Let ρ 1 , ρ 0 ∈ [0, 1]. Double standard sampling consists in observing dyads with probabilities

P(R ij = 1|Y ij = 1) = ρ 1 , P(R ij = 1|Y ij = 0) = ρ 0 . (2) Denote S o = (i,j)∈D o Y ij , So = (i,j)∈D o (1 -Y ij )
and similarly for S m , Sm . In this dyad-centered sampling design satisfying DAG (b), the log-likelihood is

log p ψ (R|Y ) = S o log ρ 1 + So log ρ 0 +S m log(1-ρ 1 )+ Sm log(1-ρ 0 ), with ψ = (ρ 0 , ρ 1 ). (3)
Definition 4 (Star sampling based on degrees -Star degree sampling). Star degree sampling consists in observing all dyads corresponding to nodes selected with probabilities {ρ 1 , . . . , ρ n } such that

ρ i = logistic(a + bD i ) for all i ∈ N where (a, b) ∈ R 2 , D i = j Y ij and logistic(x) = (1 + e -x ) -1 .
In this node-centered sampling design satisfying DAG (b), the log-likelihood is

log p ψ (R|Y ) = i∈N o log ρ i + i∈N m log(1 -ρ i ), with ψ = (a, b). (4)
Definition 5 (Class sampling). Class sampling consists in observing all dyads corresponding to nodes selected with probabilities {ρ 1 , . . . , ρ Q } such that ρ q = P(i

∈ N o | Z iq = 1) for all (i, q) ∈ N × Q.
In this node-centered sampling design satisfying DAG (d), the log-likelihood is

log p ψ (R|Z) = i∈N o q∈Q Z iq log ρ q + i∈N m q∈Q Z iq log(1-ρ q ), with ψ = (ρ 1 , . . . , ρ Q ). (5)
Identifiability of class sampling. Theorem 2 establishes the identifiability of the SBM sampled under NMAR class sampling design (see the supplementary materials for the proof). Note that the identifiability of the sampling parameters ψ = (ρ 1 , . . . , ρ Q ) and of the SBM parameters must be proved jointly because of the dependence between the network and the sampling. It is worth mentioning that both α q and ρ q are identifiable and not only their product. Although somewhat counter-intuitive, this fact is supported by the inference algorithm for class sampling in Section 3.3, which weights the recovery of the latent clusters by taking the unbalanced sampling into account.

Theorem 2. Let n ≥ 2Q and assume that for any 1 ≤ q ≤ Q, ρ q > 0, α q > 0, and that the coordinates of o = πα and t

= ( Q k=1 π 1k ρ k α k , . . . , Q k=1 π Qk ρ k α k ) are pairwise distinct.
Then, under class sampling, SBM and class sampling parameters are identifiable w.r.t. the distributions of the SBM and the sampling up to label switching.

Variational Inference

Derivations of the practical variational algorithms considerably change depending on the missing data condition at play. We start by MAR to gently introduce the variational principle for SBM, then develop algorithms in a series of NMAR conditions

MAR inference

By Proposition 1 part (i), inference in the MAR case is conducted on Y o . The EM algorithm is unfeasible since it requires the evaluation of the conditional mean of the complete log-likelihood

E Z|Y o [log p θ (Y o , Z)]
which is intractable when Y comes from an SBM. The variational approach circumvents this limitation by maximizing a lower bound of the loglikelihood based on an approximation pτ of the true conditional distribution p θ (Z|Y o ),

log p θ (Y o ) ≥ J τ,θ (Y o ) log(p θ (Y o )) -KL[p τ (Z)||p θ (Z|Y o )], = E pτ [log(p θ (Y o , Z))] -E pτ [log pτ (Z)], (6)
where τ are some variational parameters and KL is the Kullback-Leibler divergence. The approximated distribution is chosen so that the integration over the latent variables simplifies by factorization. Recall from Section 2.1 that the latent vectors Z i q = (Z i1 , . . . , Z iQ ) i∈N are independent with a multinomial prior distribution. Thus, in order to factorize the likelihood in a convenient way, a natural variational counterpart to p θ (Z|Y o ) is pτ (Z) = i∈N m(Z i q; τ i ), where τ i = (τ i1 , . . . , τ iQ ), and m(•; τ i ) is the multinomial probability density function with parameters τ i . The VEM sketched in Algorithm 1 consists in alternatively maximizing J w.r.t. τ = {τ 1 , . . . , τ n } (the variational E-step) and w.r.t. θ (the M-step). The two maximization problems are solved straightforwardly following [START_REF] Daudin | A mixture model for random graphs[END_REF]:

1. The parameters θ = (α, π) maximizing J θ (Y o ) when τ is held fixed are αq = i∈N o τiq card (N o ) , πq = (i,j)∈D o τiq τj Y ij (i,j)∈D o τiq τj .
2. The variational parameters τ maximizing J τ (Y o ) when θ is held fixed are obtained with the following fixed point relation:

τiq ∝ α q   (i,j)∈D o ∈Q b(Y ij ; π q ) τj   ,
where b(x, π) = π x (1 -π) 1-x the Bernoulli probability density function.

Algorithm 1: Variational EM for MAR inference in SBM Initialization: Set up τ (0) with some clustering algorithm repeat

θ (h+1) = arg max θ J Y o ; τ (h) , θ M-step τ (h+1) = arg max τ J Y o ; τ, θ (h+1) variational E-step until θ (h+1) -θ (h) < ε
Algorithm 1 generates a sequence {τ (h) , θ (h) ; h 0} with increasing J(Y o ; τ (h) , θ (h) ). Since there is no guarantee for convergence to the global maximum, we run the algorithm from several different initializations to finally retain the best solution.

Model selection of the number of blocks. The Integrated Classification Likelihood (ICL) criterion of [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] is relevant for latent variable models where the likelihood -and thus BIC -is intractable. [START_REF] Daudin | A mixture model for random graphs[END_REF] derive a variational ICL for the SBM which we adapt to missing data conditions: if θ = arg max log p θ (Y o , Z) then

ICL(Q) = -2E pτ log p θ(Y o , Z; Q) + Q(Q + 1) 2 log card (D o ) + (Q -1) log card (N o ) .
Note that each dyad is only counted once since we work with symmetric networks.

NMAR inference: the general case

In contrast to the MAR case, conducting inference on the observed dyads only may bias the estimates in the NMAR case. In fact, all observed data (including the sampling matrix R in addition to Y o ) must be taken into account. The likelihood of the observed data is thus log p θ,ψ (Y o , R) and the corresponding completed likelihood has the following decomposition:

log p θ,ψ (Y o , R, Y m , Z) = log p ψ (R|Y o , Y m , Z) + log p θ (Y o , Y m , Z), (7)
where an explicit form of p ψ (R|Y o , Y m , Z) requires further specification of the sampling. The joint distribution p θ (Y o , Y m , Z) has a form similar to the MAR case. Now, the approximation is required both on latent blocks Z and missing dyads Y m to approximate p θ (Z, Y m |Y o ). We suggest a variational distribution where complete independence is forced on Z and Y m , using a multinomial (resp. Bernoulli ) distribution for Z (resp. for Y m ):

pτ,ν (Z, Y m ) = pτ (Z) pν (Y m ) = i∈N m(Z i• ; τ i ) (i,j)∈D m b(Y ij ; ν ij ), (8)
where τ and ν = {ν ij , (i, j) ∈ D m } are two sets of variational parameters respectively associated with Z and Y m . This leads to the following lower bound for log p θ,ψ (Y o , R):

J τ,ν,θ,ψ (Y o , R) = E pτ,ν [log p θ,ψ (Y o , R, Y m , Z)] -E pτ,ν [log pτ,ν (Z, Y m )] .
By means of Decomposition ( 7) of the completed log-likelihood, variational approximation (8) and entropies of multinomial and Bernoulli distributions, one has

J τ,ν,θ,ψ (Y o , R) = E pτ,ν [log p ψ (R|Y o , Y m , Z)] + (i,j)∈D o (q, )∈Q 2 τ iq τ j log b(Y ij , π q ) + (i,j)∈D m (q, )∈Q 2 τ iq τ j log b(ν ij , π q ) + i∈N q∈Q τ iq log(α q /τ iq ) - (i,j)∈D m ν ij log(ν ij ) + (1 -ν ij ) log(1 -ν ij ). (9)
In ( 9), E pτ,ν [log p ψ (R|Y o , Y m , Z)] can be integrated over the variational distribution pτ,ν (Z, Y m ), as expected. The practical computations depend on the sampling design.

The general VEM algorithm used to maximize ( 9) is sketched in its main lines in Algorithm 2. Both the E-step and the M-step split into two parts: the maximization must be performed on the SBM parameters θ and the sampling design parameters ψ respectively. The variational E-step is performed on the parameters τ of the latent block Z and on the parameters ν of the missing data Y m .

Algorithm 2: Variational EM for NMAR inference in SBM

Initialisation: set up τ (0) , ν (0) and ψ (0) repeat θ (h+1) = arg max θ J Y o , R; τ (h) , ν (h) , ψ (h) , θ M-step a) ψ (h+1) = arg max ψ J Y o , R; τ (h) , ν (h) , ψ, θ (h+1) M-step b) τ (h+1) = arg max τ J Y o , R; τ, ν (h) , ψ (h+1) , θ (h+1) VE-step a) ν (h+1) = arg max ν J Y o , R; τ (h+1) , ν, ψ (h+1) , θ (h+1) VE-step b) until θ (h+1) -θ (h) < ε
Interestingly, resolution of the two steps concerned with the optimization of the parameters related with the SBM -that is to say, θ and τ -can be stated almost independently of any further specification of the sampling design.

Proposition 3. Consider the lower bound J τ,ν,θ,ψ (Y o , R) given by (9).

1. The parameters θ = (α, π) maximizing (9) when all others are held fixed are

αq = 1 n i∈N τiq , πq = (i,j)∈D o τiq τj Y ij + (i,j)∈D m τiq τj νij (i,j)∈D τiq τj .
2. The optimal τ in (9) when all other parameters are held fixed verifies

τiq ∝ λ iq α q   (i,j)∈D o ∈Q b(Y ij ; π q ) τj     (i,j)∈D m ∈Q b(ν ij ; π q ) τj  
with λ iq a simple constant depending on the sampling design.

Proof. These results are simply obtained by differentiation of (9).

The two steps concerned with ψ and ν are specific to the sampling designs used to describe R. Further details are provided below for the designs presented in Section 2.3.2.

NMAR: specificities related to the choice of the sampling

In light of Figure 1, NMAR conditions specified by DAGs (b), (c) or (d) induce different simplifications for the conditional distribution of the sampling design R:

DAG (b) p ψ (R|Y o , Y m , Z) = p ψ (R|Y o , Y m ), DAG (c) p ψ (R|Y o , Y m , Z) = p ψ (R|Y o , Y m , Z), DAG (d) p ψ (R|Y o , Y m , Z) = p ψ (R|Z).
This induces different evaluations of E pτ,ν [log p θ,ψ (Y o , R, Y m , Z)] in the lower bound (9) for double standard sampling, star degree sampling and class sampling. We obtain below explicit formulas of ψ and ν by differentiation of the corresponding variational lower bounds. The computations are tedious but straightforward and thus eluded in the following. Double-standard sampling. Let s m = (i,j)∈D m ν ij , sm = (i,j)∈D m (1 -ν ij ) be the variational counterparts of S m and Sm . From (3) we have

E p log p ψ (R|Y ) = S o log ρ 1 + So log ρ 0 + s m log(1 -ρ 1 ) + sm log(1 -ρ 0 ).
Proposition 4 (double standard sampling).

1. The parameters ψ = (ρ 0 , ρ 1 ) maximizing (9) when all others are held fixed are

ρ0 = So So + sm , ρ1 = S o S o + s m . ( 10
)
2. The optimal ν in (9) when all other parameters are held fixed are

νij = logistic   log 1 -ρ 1 1 -ρ 0 + (q, )∈Q 2 τ iq τ j log π q 1 -π q   . Moreover, λ iq = 1 ∀(i, q) ∈ N × Q for optimization of τ in Proposition 3.b).
Class sampling. According to (5) we have

E p log p ψ (R|Y ) = i∈N o q∈Q τ iq log(ρ q ) + i∈N m q∈Q τ iq log(1 -ρ q ).
Proposition 5 (class sampling).

1. The parameters ψ = (ρ 1 ...ρ Q ) maximizing (9) when all others are held fixed are

ρq = i∈N o τ iq i∈N τ iq . ( 11
)
2. The optimal ν in (9) when all other parameters are held fixed verify

νij = logistic   (q, )∈Q 2 τ iq τ j log π q 1 -π q   . Moreover λ iq = ρ 1 {i∈N o } q (1 -ρ q )
1 {i∈N m } for optimization of τ in Proposition 3.b).

Star degree sampling. From Expression (4) of the likelihood, one has

E p log p ψ (R|Y ) = - i∈N m a + b Di + i∈N E p -log(1 + e -(a+bD i ) ) , where Di = E p [D i ] = i∈N m ν ij + i∈N o Y ij is the approximation of the degrees. Because E p -log(1 + e -(a+bD i )
) has no explicit form, an additional variational approximation is needed [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF]. This technique was recently used in random graph framework [START_REF] Latouche | Goodness of fit of logistic models for random graphs[END_REF]. It relies on the following approximation of the logistic function:

g(x) ≥ g(ζ) + x -ζ 2 + h(ζ)(x 2 -ζ 2 ), h(ζ) = -1 2ζ logistic(ζ) - 1 2 (12)
for all (x, ζ) ∈ R × R + . This leads to a lower bound of the initial lower bound:

log p θ,ψ (Y o , R) ≥ J τ,ν,θ,ψ (Y o , R) ≥ J τ,ν,ζ,θ,ψ (Y o , R), (13)
with ζ = (ζ i , i ∈ N ) such that ζ i > 0 is an additional set of variational parameters used to approximate -log(1+e -x ). The second lower bound J τ,ν,ζ,θ,ψ is derived from Equation ( 12) and given in the supplementary materials for completeness . At the end of the day, we have an additional set of variational parameters to optimize, and a corresponding additional step in Algorithm 2. Expression of all the parameters specific to star degree sampling by differentiating J τ,ν,ζ,θ,ψ .

Proposition 6 (star degree sampling).

Let D 2 i = E p D i 2 and D- k = Dk -ν k . 1. The parameters ψ = (a, b) maximizing J τ,ν,ζ,θ,ψ (Y o , R) when others are held fixed are b = 2 n 2 -card (N m ) i∈N (h(ζ i ) Di ) -1 2 i∈N Di -i∈N m Di × i∈N h(ζ i ) 2 i∈N (h(ζ i ) D 2 i ) × i∈N h(ζ i ) -2 i∈N h(ζ i ) Di 2 , â = - b i∈N h(ζ i ) Di + n 2 -card (N m ) i∈N h(ζ i )
.

2. The parameters ζ maximizing J τ,ν,ζ,θ,ψ (Y o , R) when others are held fixed are

ζi = a 2 + b 2 D 2 i + 2ab Di , ∀i ∈ N .
3. The optimal ν in J τ,ν,ζ,θ,ψ (Y o , R) when all other parameters are held fixed verify

νij = logistic (q, )∈Q 2 τ iq τ j log π q 1 -π q -b + 2h(ζ i ) ab + b 2 (1 + D-j i ) + 2h(ζ j ) ab + b 2 (1 + D-i j )
. ( 14)

Moreover, λ iq = 1 ∀(i, q) ∈ N × Q for optimization of τ in Proposition 3.b).
Model selection. In NMAR cases, ICL can be useful not only to select the appropriate number of blocks but also for selecting the most appropriate sampling design when it is unknown. Contrary to the MAR case, ICL is no longer a straightforward generalization of [START_REF] Daudin | A mixture model for random graphs[END_REF]. Indeed, the complete likelihood and thus the penalization needs to take into account the sampling design. Let us consider a model with Q blocks and a sampling design with K parameters (i.e. the dimension of ψ). The ICL criterion is a Laplace approximation of the complete likelihood p(Y o , Y m , R, Z|Q, K) with p(θ, ψ|Q, K) the prior distributions on the parameters such that

p(Y o , Y m , R, Z|Q, K) = Θ×Ψ p θ,ψ (Y o , Y m , R, Z|Q, K)p(θ, ψ|Q, K)dθdψ.
Proposition 7. For a model with Q blocks, a sampling design with a vector of parameters ψ ∈ R K and ( θ, ψ) = arg max (θ,ψ) log p θ,ψ (Y o , Y m , R, Z), then

ICL(Q) = -2E pτ,ν ; θ, ψ log p θ, ψ(Y o , Y m , R, Z|Q, K) + pen ICL (Q), pen ICL =    K + Q(Q+1) 2 log n(n-1) 2 + (Q -1) log(n) for dyad-centered sampling Q(Q+1) 2 log n(n-1) 2 + (K + Q -1) log(n) for node-centered sampling
Note that an ICL criterion for MAR sampling designs can be constructed in the same fashion for the purpose of comparison with NMAR sampling designs.

MAR condition

Algorithm 1 for MAR samplings is tested on affiliation networks with 3 blocks. The number of blocks is assumed to be known. For this topology the probability of connection within a block is η and is ten times stronger than the probability of connections between nodes from different blocks. We generate networks with n = 200 nodes and marginal probabilities of belonging to blocks α = (1/3, 1/3, 1/3). The sampling design is chosen as a random-dyad sampling with a varying ρ. The difficulty is controlled by two parameters: the sampling effort ρ and the overall connectivity in matrix π, defined by c = q α q α π q , which is directly related to the choice of η: the lower the η, the sparser the network and the harder the inference. The simulation is repeated 500 times for each configuration (c, ρ). Figure 2 displays the results in terms of estimation of π and of classification recovery, for varying connectivity c and sampling effort ρ. Our method achieves good performances even with a low sampling effort provided that the connectivity is not too low. The adjacency matrix Y is generated under random-dyad sampling strategy for various connectivity c = q α q α π q .

π -π F / π F Adjusted Rand Index

NMAR condition

Under NMAR conditions we conduct an extensive simulation study by considering various network topologies (namely affiliation, star and bipartite), the connectivity matrix of which are given in Figure 3. We use a common tuning parameter to control the connectivity of the networks in each topology: the lower the , the more contrasted the topology. Among the three schemes developed in Section 2.3.2, we investigate thoroughly the double standard sampling, for which we exhibit a large panel of situations where the gap is large between the performances of the algorithms designed for MAR or NMAR cases. Other sampling designs are explored in the supplementary materials.

Simulated networks have n = 100 nodes, with varying in {0.05, 0.15, 0.25}. Prior probabilities α are chosen specifically for affiliation, star and bipartite topologies, respec- (1/3, 1/3, 1/3), (1/6, 1/3, 1/6, 1/3) and(1/4, 1/4, 1/4, 1/4). The exploration of the sampling parameters ψ = (ρ 0 , ρ 1 ) is done on a grid [0.1, 0.9] × [0.1, 0.9] discretized by steps of 0.1. Algorithm 2 is initialized with several random initializations and spectral clustering.

  1 - 1 - 1 -   (a) affiliation     1 - 1 - 0 0 1 - 0 0 0 1 - 1 - 0 0 1 - 0     (b) star     1 - 1 - 1 - 1 -     (c) bipartite
In Figure 4, the estimation error is represented as a function of the difference between the sampling design parameters (ρ 0 , ρ 1 ): the closer this difference to zero, the closer to the MAR case. As expected, Algorithm 1 designed for MAR only performs well when ρ 1 -ρ 0 ≈ 0. Algorithm 2 designed for NMAR double-standard sampling shows relatively flat curves which means that its performances are roughly constant no matter the sampling condition.

Figure 5 reports estimation accuracy for the sampling parameters ρ 0 and ρ 1 . Results show a good ability of the VEM to estimate these parameters. As expected, performances deteriorate for uncontrasted topologies with low sampling rate.

Model selection. Simulations are also conducted to study the performances of ICL. We compare results for the different topologies described in Figure 3 for = 0.05. Rates of correct answers for selecting the number of blocks Q under a double standard sampling with different sampling rates are displayed in Table 1. The ARI is also provided. The ICL shows a satisfactory ability to select the true Q even if the selection task obviously needs a larger sampling effort than the estimation task. It is worth mentioning that a whole block may not be sampled, which leads the ICL to select a lower number of blocks. In such a case the ARI is a meaningful additional information to demonstrate that the clustering remains coherent with the true one.

sampling rate affiliation bipartite star (0.154, 0.405] 0.58/0.96 0.46/0.84 0.45/0.84 (0.405, 0.656] 0.95/0.99 0.87/0.98 0.90/0.98 (0.656, 0.908] 1/1 0.99/1 0.99/1

Table 1: Performance of the ICL criterion: rates of correct answers when choosing the number of blocks and Adjusted Rand Indexes. Tested configurations are different sampling rates and three topologies (affiliation, bipartite and star) under a double standard sampling. Each configuration is simulated 500 times.

In Table 2, results concern the rates of correct selections of the sampling design when the two designs in competition are the random-dyad and the double standard samplings. As expected, the rate of correct answers increases with the sampling rate. 5 Importance of accouting for missing values in real networks

Seed exchange network in the region of Mount Kenya

In a context of subsistence farming, studies which investigate the relationships between crop genetic diversity and human cultural diversity patterns have shown that seed exchanges are embedded in farmers' social organization. Data on seed exchanges of sorghum in the region of Mount Kenya were collected and analyzed in [START_REF] Labeyrie | Seed exchange networks, ethnicity, and sorghum diversity[END_REF][START_REF] Labeyrie | Influence of ethnolinguistic diversity on the sorghum genetic patterns in subsistence farming systems in eastern kenya[END_REF]. The sampling is node-centered since the exchanges are documented by interviewing farmers who are asked to declare to whom they gave seeds and from whom they receive seeds. Since an interview is time consuming, the sampling is not exhaustive. A limited space area was defined where all the farmers were interviewed. The network is collected with missing dyads since information on the potential links between two farmers who were cited but not interviewed is missing. With the courtesy of Vanesse Labeyrie, we analyzed the Mount Kenya seed exchange network involving 568 farmers among which 155 were interviewed. Although other farmers in this region might be connected to non-interviewed farmers, we focus on this closed network of 568 nodes.

Since we only know that the sampling is node-centered, we fit SBM under the three node-centered sampling designs presented in Section 2.2 (star (MAR), class and star degree sampling). The ICL criterion is minimal for 10 blocks under the star degree sampling and for 11 blocks under the class degree sampling. The clusterings between the SBMs obtained with either class or star degree sampling remain close from each other (ARI: 0.6) and both unravel a strong community structure. The model selected by ICL for MAR sampling is composed by 11 blocks. The ARIs between MAR clustering and the two other clusterings are lower (around 0.4). Finally, note that interviewed and non-interviewed farmers are mixed up in the blocks of the three selected models. The ICL criteria computed for the three sampling designs are a slightly in favor of the MAR sampling.

On top of network data, categorical variables are available for discriminating the farmers such as the ntora2 they belong to (10 main ntoras plus 1 grouping all the others) and the dialect they speak (4 dialects). In Figure 6, we compute ARIs between the ntoras (left panel), the dialects (right panel) and the clusterings obtained with the SBM under the three node-centered sampling designs for a varying number of blocks. Even though the ARIs remain low, the clusterings from class or star degree sampling seem to catch a non negligible fraction of the social organization, larger than the one caught by the clustering from the MAR sampling. These two categorical variables, reflecting some aspects of the social organization, could partially explain the structure of the exchange network.

ARI q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.05 0.10 0.15 4 8 12 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.05 0.10 4 8 12 number of blocks sampling q q q Class Degree MAR

Figure 6: ARIs computed between the clusterings given by an SBM under class, star degree and MAR samplings with a varying number of blocks Q and ntora of farmers (left-hand-side) or dialect spoken by farmers (right-hand-side)

ER (ESR1) Protein-Protein Interaction network in breast cancer

Estrogen receptor 1 (ESR1) is a gene that encodes an estrogen receptor protein (ER), a central actor in breast cancer. Uncovering its relations with other proteins is essential for a better understanding of the disease. To this end, various bioinformatics tools are available to centralize knowledge about possible relations between proteins into networks known as Protein-Protein Interaction (PPI) networks. The platform string [START_REF] Szklarczyk | String v10: protein-protein interaction networks, integrated over the tree of life[END_REF] accessible via http://www.string-db.org is one of the most popular tools for this task.

Given a set of one (or several) initial protein(s) provided by the user, it is possible to recover a valued network between all proteins connected to the initial set. The value of an edge in this network corresponds to a score obtained by aggregating different types of knowledge (wet-lab experiments, textmining, co-expression data, etc. . . ), reflecting a level of confidence. Thus, it is possible for a given protein -we choose ER here -to recover the PPI network between all proteins involved. Our ambition is to rely on a SBM with missing data to finely analyze such networks: we rather describe a dyad as missing (thus not choosing between 0 or 1) if its level of confidence is too low. The PPI network in the neighborhood of ER is composed by 741 proteins connected by edges with values in (0, 1]. We remove ER from this set of proteins, as well as the zinc finger protein 44. Indeed, they were both connected to most of the other proteins and would thus only blur the underlying clustering structure. We denote ω ij the weight associated with dyad (i, j). By means of a tuning parameter γ reflecting the level of confidence, the adjacency matrix is defined as follows:

A γ = (A γ ) ij =    1 if ω ij > 1 -γ, NA if γ ≤ ω ij ≤ 1 -γ, 0 if ω ij < γ. (15)
In order to analyze the ER-centered network, Algorithm 1 (random-dyad MAR sampling) and Algorithm 2 (double-standard NMAR sampling) were applied on A γ for γ varying in {.15, .25, .35}, hence taking the uncertainties on the missing dyads into account with various thresholds. The ICL criterion in Figure 7 systematically chooses the NMAR modeling against the MAR modeling, whatever the value of γ. We study the best MAR and NMAR models associated with γ = 0.35, which value exhibits a clearer choice of the ICL than for γ = {.15, .25} for both MAR and NMAR modelings. The two corresponding SBMs have 11 clusters for MAR sampling and 13 clusters for NMAR sampling. The ARI between the two clusterings is around 0.39: this is mainly due to a large block in the random-dyad MAR clustering which contains much more nodes than any of the blocks in the NMAR clustering. The latter dispatches many of these nodes in four blocks (see the supplementary materials for a more detailed exposition of results). To prove that this finest clustering of the nodes is more relevant from the biological point of view, we propose a validation based on external biological knowledge. To this end, we rely on the Gene Ontology (GO) annotation [START_REF] Ashburner | Gene ontology: tool for the unification of biology[END_REF] which provides a DAG of ontologies to which genes are annotated if the proteins encoded by these genes are involved in a known biological process. Here, we use GO to perform enrichment analysis (that is to say identifying classes of genes that are over-represented in a large set of genes, via a simple hypergeometric test) on genes corresponding to the proteins present in the large block for MAR, and the corresponding four blocks for NMAR. Interestingly, at a significance level of 1%, we find a single significant biological process for MAR modeling while 13 are found significant in the NMAR case. We check that it is not due to a simple threshold effect by looking at the ranks of the p-values of the 13 NMAR significant processes in the 100 first most significant terms found in the MAR model: only 5 of the NMAR processes are found, with high ranks (24, 33, 39, 56 and 77) far from the smallest MAR p-values.

Conclusion

This paper shows how to deal with missing data on dyads in the SBM. We study MAR and NMAR sampling designs motivated by network data and propose variational approaches to perform inference under this series of designs, accompanied with model selection criteria. Relevance of the method is illustrated on numerical experiments both on simulated and realworld networks. An R-package is available at https://github.com/jchiquet/missSBM. . This work focuses on undirected binary networks. However, it can be adapted to other SBMs, in particular those developed in [START_REF] Mariadassou | Uncovering latent structure in valued graphs: A variational approach[END_REF] for (un)directed valued networks with a distribution of weights belonging to the exponential family. It could also be adapted to the degree-corrected SBM [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF], where the sampling design would depend on the degree correction parameters. This should lead to a design close to the star degree sampling. In future works, we plan to investigate the consistency of the variational estimators of SBM under missing data conditions, looking for similar results as the ones obtained in [START_REF] Bickel | Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels[END_REF] for fully observed networks. Another path of research is to consider missing data where we cannot distinguish between a missing dyad and the absence of an edge like in [START_REF] Priebe | Statistical inference on errorfully observed graphs[END_REF] and [START_REF] Balachandran | On the propagation of low-rate measurement error to subgraph counts in large networks[END_REF].

Figure 2 :

 2 Figure2: Estimation error of π and Adjusted Rand Index averaged over 500 simulations in the MAR setting. The adjacency matrix Y is generated under random-dyad sampling strategy for various connectivity c = q α q α π q .

Figure 3 :

 3 Figure 3: Matrix π in different topologies with inter/intra block probabilities.

Figure 4 :

 4 Figure 4: Double standard setting: estimation error of π and adjusted Rand index averaged over 500 simulations for affiliation, bipartite and star topologies.

|ρ 0 -Figure 5 :

 05 Figure 5: Double standard setting: estimation error of ρ 0 and ρ 1 averaged over 500 simulations for affiliation, bipartite and star topologies.

Figure 7 :

 7 Figure 7: ICL criteria for SBMs with random-dyad MAR sampling and double-standard NMAR sampling in the thresholded ER network.

Table 2 :

 2 Rates of correct answers of the ICL criterion when choosing between Random-dyad sampling (MAR) and double standard sampling (NMAR) in each of the 18 configurations. A configuration is the combination of a topology (affiliation, bipartite and star), a sampling rate and a sampling design. Each configuration is simulated 500 times.

	sampling rate sampling affiliation bipartite star
	(0.096, 0.367]	MAR	0.73	0.67	0.63
		NMAR	0.72	0.75	0.75
	(0.367, 0.638]	MAR	1	1	1
		NMAR	0.91	0.78	0.82
	(0.638, 0.909]	MAR	1	1	1
		NMAR	0.91	0.8	0.95

More complex sampling schemes -for instance adversarial strategies -are thus not handled

Simulation studyIn this section, we illustrate the relevance of our approaches on network data simulated under the SBM and sampled under MAR and NMAR conditions. The quality of the inference is assessed by computing the distance between the estimated and the true connectivity matrices π in terms of Frobenius norm. The quality of the clustering recovery is measured with the adjusted Rand index (ARI,[START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF] between the true classification and the clustering obtained by maximum posterior probabilities for each τ i .

The ntora is a small village or a group of neighborhoods

Acknowledgment. The authors thank Sophie Donnet (INRA-MIA, AgroParisTech) and Mahendra Mariadassou (INRA-MaIAGE, Jouy-en-Josas) for their helpful remarks and suggestions. We also thank all members of MIRES for fruitful discussions on network sampling designs and for providing the original problems from social science. In particular, we thank Vanesse Labeyrie (CIRAD-Green) for sharing the seed exchange data and for related discussions on the analysis. This work is supported by the INRA MetaProgram « GloFoods » through the project « SEARS » and by two public grants overseen by the French National research Agency (ANR): first as part of the « Investissement d'Avenir » program, through the « IDI 2017 » project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02, and second by the « EcoNet » project.


Anna Carla Araujo
email: araujo@insa-toulouse.fr

Rossana Mara Da Silva
email: rossana@metalmat.ufrj.br

Moreira Thiré

Damien Coutouly
email: damiencoutouly@insa-toulouse.fr

Anas Rais
email: anasrais@insa-toulouse.fr

Léo Sudriès
email: sudries@insa-toulouse.fr

Tom Gaboriaud
email: gaboriau@insa-toulouse.fr

Giovanna Siqueira
email: giovannassa@poli.ufrj.br

Sardela Andrade

Raphael Értola
email: raphaelertola@poli.ufrj.br

Pereira De

Vitor Deus Santos

Amadeu Corrêa

COLLABORATIVE PROJECT AS A DIDACTIC PLATFORM OF VIRTUAL MOBILITY FOR ADDITIVE MANUFACTURING COURSE

Keywords: Industry 4.0, Internacional Virtual Mobility, Additive Manufacturing, FDM, Composite filament.)

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The additive manufacturing (AM) is a recent manufacturing process that has been growing up the last 10 years and it is getting to production scale [START_REF] Campos | The influence of additive manufacturing on the micromilling machinability of ti6al4v: A comparison of slm and commercial workpieces[END_REF]. Its inclusion in Engineering courses is mandatory due to the technological contribution that this concept causes to workpiece design, to the development of new materials [START_REF] Coelho | Manufacturing of gypsum-sisal fiber composites using binder jetting[END_REF] and to adapt the industrial chain as it is one of the nine technological pillars of Industry 4.0 [START_REF] Motyl | How will change the future engineers' skills in the industry 4.0 framework? a questionnaire survey[END_REF].

There is a needed for the development of Industrial Internet of Things in the Engineering course [START_REF] Louw | Teaching industry 4.0 technologies in a learning factory through problem-based learning: case study of a semi-automated robotic cell design[END_REF] with a strong basis of manufacturing processes. "Modern manufacturing education plays a key role in preparing qualifies workforce for the manufacturing industry", as claimed [START_REF] Ngaile | Challenges in teaching modern manufacturing technologies[END_REF] detailing that they need new innovation, technologies and customer-driven products. [START_REF] Motyl | Trends in engineering education for additive manufacturing in the industry 4.0 era: a systematic literature review[END_REF] discussed about the educational strategies adopted at university level that involve AM. They made a cartography of 30 articles in journals and conferences that studied this topic and they claimed that "what emerges is a strong interest in improving the education of engineering students through the introduction of courses or workshops dedicated to Industry 4.0 subjects, in particular to manufacturing improvement and to the use and deepening of AM techniques for production. (...) and the need to introduce new skills."

The classical way for inclusion of a new process in the curriculum is to add a topic in the current course of "Introduction to Manufacturing Processes" [START_REF] Groover | Introdução aos Processos de Fabricação[END_REF], to create optional courses of the new thematic and to promote the thematic in the final engineering project subjects [START_REF] Coelho | Análise Dimensional de Peças Produzidas por Manufactura Aditiva[END_REF][START_REF] Gutnik | Fabricação de Protótipo para Reabilitação Nervosa por Manufatura Aditiva com Base Líquida[END_REF]. Although, the changes that this process causes to the way of designing parts and assembly mechanisms, should include the interdisciplinary in Industrial Engineering, as mechanical design and metrology in order to impact the development of a specialized profile [START_REF] Ngaile | Challenges in teaching modern manufacturing technologies[END_REF].

The new curriculum should include the classical FDM (Fused Deposition Modeling) process, the most economic AM process, but also the other AM strategies, as powder bed (SLS, SLM), selective powder deposition and photo-curing polymers (SLA) process [START_REF] Magaldi | Comparison of the porosity of scaffolds manufactured by two additive manufacturing technologies: Sla and fdm[END_REF]. It should provide discussing about what process will be more adapted to each project and material, as shown by [START_REF] Diegel | Teaching Design for Additive Manufacturing Through Problem-Based Learning[END_REF][START_REF] Flores Ituarte | Additive Manufacturing Validation Methods, Technology Transfer Based on Case Studies[END_REF]. [START_REF] Havenga | Preparing first-year engineering students for cooperation in real-world projects[END_REF] proposed an problem-based (active learning) course in the first year of Engineering for cooperation in real word projects. This article claims that the student works together solving multidisciplinary Engineering problems developing inquiry, critical dialog, creative thinking and active participation. In INSA Toulouse, the first year has a common curriculum for all Engineering courses with an specific module that includes AM in first year: the students should design a small piece (printing time under 25 min) in CAD software and send it to the 3D printers. This activity was chosen to fit on the available time of the students in one of the two practical manufacturing classes. So, adding this activity to the last year of Engineering could use the knowledge that they developed in the five years of Engineering and explore more possibilities of collaboration inside the group. On the contrary, students of Materials engineering and Metallurgical Engineering courses in UFRJ only take contact with additive manufacturing techniques when they decided to develop final projects or training activities within specific research groups of the Department that work with these techniques. Up to now, there is no specific (elective or compulsory) subject within their academic curriculum related to this topic.

In order to overcome these academic gaps, we conducted a virtual project course in additive manufacturing for a group composed of three students from the fifth-year of mechanical engineering students in a French Institution (INSA-Toulouse) and four fourth-year students from Metallurgical Engineering and Materials Engineering in a Brazilian University (UFRJ). This original educational approach aimed to bring together students from different establishments in the form of virtual meetings to build a mechanical project with an online manufacturing workshop. The target of the project was partially defined by the professors who guided through deadlines and technical suggestions.

The objective of this article is to share the experience of this course and motivate professors to introduce actual tools that prepare our Engineers for the international connected industries developing soft skills of: integration from different scientific fields, multicultural communication and contributing to the practice of the technical English language in a professional situation.

The structure of this article includes: the virtual mobility proposition to the partners, the organization before the course, the description of the students internal organization inside the group, a brief overview of the technical development and its results, and, in the last section, the pedagogical results and conclusions.

COURSE ORGANIZATION WITH PARTNERS

The program of the standard "Advanced additive manufacturing" module at INSA-Toulouse is currently presented in the format: a lecture introducing the main additive manufacturing processes and a review of the process basics, three courses in "conference" format presented by two manufacturers and by a researcher in the field and, then, 7 weeks of experimental and supervised projects in groups of 5 students. The last session will be devoted to project presentations by all groups, including the presentations related to this project that will be done by video-conference for class.

The inclusion of the virtual mobility in this course was gradual : in the first year of implementation, the classes were developed from October 2020 to January 2021 between INSA-Toulouse and Universitá di Trento supported by a project entitled "3Di Project -Teaching the industry of the future through international monitoring of the FDM additive manufacturing process". The objective was to create a hybrid course with teaching mobility when manufacturing could be done in both countries. The aim of the project was to reproduce a small structure of a multinational industry with two plants and one common project.

This project was submitted in 2019 (before COVID-19 crisis) for financial support inside the call "IDEFI DEFI Diversités 2020"1 from Toulouse University to adapt teaching practices to the diversity of student populations. The common actions carried out by this project aim to to develop the skills expected by the socioeconomic world.

The financial support for this project was approved in 2020, when the physical part of the collaboration could not be developed due to the mobility limitation due to the lock-down. The first exercise was entirely virtual for the development of the mechanical design, although the final manufacturing could be completed in each country, producing two products: one in Italy and in France.

In this new 2021 version, the objective was to include a third partner: one university in South America without physical student mobility, enlarging the collaboration and repeating the experience. In 2021, two groups from the French class should interact with Italian students and one group with the Brazilian students, activity described in this article. In the case of French-Brazilian group, an adaptation of pedagogic method was necessary, since due to the restrictions imposed by the Covid-19 pandemic, the Brazilian students were not able to carry out the course activities in person at UFRJ, such as Italian students did in the first version.

The academic and industrial lectures were available on an open Moodle platform (Open Insa2 ) or it was developed with synchronous meetings using Zoom Software so that it could be open to other establishments. The main objective provided for this group was: "Develop challenge activity with a final target product that uses an extruder to produce your own polymer filament".

The students should carry out the project in 7 sections:

• Section 1: Presentation videoconference and Brainstorming between international students

• Section 2: Definition of the goals of the project and discovery of the new equipment

• Section 3: Presentation of extrusion results for the supervisors

• Section 4: Improvement of the experimental part (corrections and progresses)

• Section 5: Printing of the final product

• Section 6: Elaboration of the presentation of the project results

• Section 7: Presentation of the project online to the three institutions, including the students that do not participate in the project in INSA-Toulouse.

The expected benefits of the project were related to the learning of students in the field of manufacturing and remote monitoring, the international interaction of trainers and students and the dissemination of partnership to other students who are not in mobility.

INTERNAL ORGANIZATION

Internal meeting and organization of the responsibilities

The students defined the title "FDM analysis on the change of mechanical properties of PLA", the demand (or the specific challenge) as "3D-printed bicycle storage support" and the technical objective to achieve : "Increase the Young's modulus of the material and increase the maximum stresses" comparing the regular filament with a new proposed one, manufactured by their own experience. As both institutions explained the target, students from both institutions understood the same objective, although each side had specific methodology for practical courses.

In the beginning, even after defining their own way of communication and registering the information, the main fail on the synchronization of progress was related to the difference between the core of the activities: the French students having the practical difficulties onsite and the Brazilian students trying to be useful and define their responsibilities, but online, far from the problem. At this point it was very important to insist that they would have to define each student's responsibility and that the grade will depend on the homogeneity of the work. The mechanical engineering students had practical extrusion problems that the Brazilian students had the answers and the skill to be able to guide to the solutions. The professors tried to guide and identify the mechanics of the group collaboration and after this point, they were able to develop a structured project with more information flow. The students will be faced with a similar situation in an international company with technological transfer.

The Brazilian group described in their report: "Furthermore, as a secondary objective, but no less important, we should develop interpersonal skills, which would encompass both the normal difficulties of a multidisciplinary project -in this case, materials, metallurgical and mechanical engineering -as well as the problems of online interaction between students from different cultures in different time zones, communicating in a language other than their mother tongue. (. . . ) we realized throughout the course that we could merge our knowledge in order to seek the best development of the project, even in topics that were a little outside our area."

Analysis of the materials and equipments

The students were free to choose the material for manufacturing the filament. They could use recycled pieces, commercial filaments or polymer pellets. In the beginning of the project, they started evaluating the possibility of producing a filament with magnetic properties. This idea was a cutting edge proposition but they realised that there were practical experimental problems related to a new good idea, what gave them the insight of what could be done in a RD department.

Finally they chose to produce a composite filament with theoretical higher rigidity in order to produce a bicycle support that needs a higher strength. A composite filament with PLA matrix, widely used in 3D printing and available on stock, reinforced by fibers. We had two fibers available on laboratory stock, as we work with aerospace research: glass and carbon fibers. Both materials have a high Young modulus but they decided that carbon fiber offered fewer risks in terms of safety while handling it. This shows that they understood that several components for decision making in Industry Industry should be considered.

The available equipment for the project included an extrusion filament equipment (Filafab 3D3) with Spool Winder (FilaFab Mk3) and a dedicated FDM printer (Ender 3), so there was no impact on other projects. The equipment can be seen in Figure 1.

The produced filament was used to print samples for tensile tests with dog-bone shape type V (ASTM D638-10 Standard). Five workpiece of each condition were printed. Table 1 shows some experimental parameters used to produce the PLA/carbon fibers filaments and printed parts. Table 1: Parameters used to produce filaments and printed parts

Theoretical material behavior

The theoretical background of the students is linked to elastic behavior of homogeneous material to simulate the mechanical behavior of composites. They did not consider other parameters, for example, friction forces and temperature in the extrusion head and its effects.

The students the elasticity modulus using the Young modulus of each composite material. The students used the properties values from the literature (approximate values). According to [START_REF] Zweben | Comprehensive Composite Materials II[END_REF], the theoretical Young modulus of carbon fiber is between about 35-1000 GPa, depending on fiber orientation (along the axial axis or layer plane). They used the average of 517 GPa for fibers and between 3-4 GPa for PLA, based on a comparative study [START_REF] Santana | A comparative study between petg and pla for 3d printing through thermal, chemical and mechanical characterization[END_REF]. Regarding the mass ratio, they chosen 1 wt% and 5 wt% fiber after consulting some commercial websites.

They chose two models for calculating the composite rigidity:

1. Krenchel model: based on the mixing rule, including an experimental factor η related to the orientation of the fibers inside the matrix. They used η = 3/8, because it is the value for random fibers and homogeneously distributed in the lamina.

The composite elastic module E c based on Krenchel model is dependant on the elastic modulus of the matrix E M and the fibers E F and also its volumetric ratio V M and V F :

E c = E M V M + ηE F V F (1)
2. Halpin-Tsai Model: It is used to predict elastic properties of the material based on an empirical model. It presents an efficiency factor ξ related to the length of the fibers l F and its diameter d F : xi = 2l F d F And a factor η that relates the behavior of the elastic properties of the matrix and the fiber, as the previous model.

Eq. ( 2) presents the composite elastic module E c calculation using this model:

E c = E M 1 + ξηV F 1 -ηV F (2)
The theorical values obtained, respectively, are 5.44 GPa and 5.73 GPa, considering 1 wt% of carbon fiber. It is important to highlight that these are theorical results, which means that every mathematical model has approximations. For instance, the Krenchel model presuppose that the fibers are homogeneously and randomly distributed. So, in this model, the fiber orientation and distribution play an important role when the elastic property is calculated. On the other hand, the Halpin-Tsai is an empirical model that can foresee satisfactorily if the fibers are well oriented in the main loading direction. Once again, the fibers orientation plays decisive role. The students faced some experimental problems:

Manufacturing workpieces and mechanical tests

• The extruder was quite small, without any technical after sales service. So, the students in France faced a lot of problems to regulate the flow and to achieve a constant diameter of the filament. Therefore, the single screw extruder does not reach high values of pressure and shear rate, which impaired the production of a composite. As a solution, the students proposed a low volume fraction of the reinforcement fibers and/or to increase the extrusion temperature to improve the fluidity of the polymer, although this increasing could led to the degradation of PLA molecules.

• In the beginning, they observed that the pure PLA filament produced using polymer pellets was getting thin and brittle. They decided to use cut pieces of a PLA commercial filament as raw material for extrusion process instead of polymer granules in order to benefit of the additives the industry includes in their products. . The cooling process of the extruded filament was also evaluated. They proposed to change the cooling based on the theoretic basis on polymer behavior. The filament was cooled in air.

• Adding carbon fiber to the PLA, they had the difficulty to separate the fibers and mixing it in a homogeneous way inside the container of the extruder, before the hot channel. To overcome this problem, they reduced the cut pieces of FDM filament to improve mixture capability, and then, mixed with carbon fibers before feeding the extruder.

• For the printing stage, they realize that the printing parts with home produced filament should take into account the diameter they achieved in extrusion: space between layers should compensate the filaments diameter (smaller than expected). The 3D printer was setup for filaments with 1.75 mm in diameter. However, thinner filaments (1.4 mm in diameter) were produced. Thus, students realized that an adjustment in machine was necessary to avoid under-extrusion during printing process. They adapted the FDM parameters to be able to produce the additive manufacturing pieces (bed temperature, nozzle temperature, speed and printing trajectory). Since a large difference between the maximum stresses values was observed, the printing speed were changed to improve the adhesion between the layers.

Figure 3 shows the stress-strain curves of the printed workpieces produced using extruded filaments with and without fibers. According to these results, the incorporation of fibers did not increase the Young modulus. Moreover, composite parts presented lower ductility. The average value of Young Modulus for the composite samples was 2.2 GPa, which is lower than those of the pure PLA. The results obtained by the theoretical models are not comparable to those obtained experimentally, , mainly because the carbon fibers probably agglomerated throughout the filament extruded. One justification could be that the maximum extruder's shear rate was not enough to promote a proper distribution of the filler, which could weaken the samples by creating regions of tension concentration.

Final presentation and report

In the end of the semester, an online presentation for the other students of this course (other groups not in virtual mobility) and for the professors from France, Brazil and Italy, as it can be observed in Fig. 4. Also, they produced one poster for the evaluation in INSA Toulouse (Fig. 4b) and one report in Portuguese for their evaluation in Brazil (Fig. 4c).

Figure 3: Comparison of mechanical behavior between some samples with and without carbon fibers.

CONCLUSIONS

Clearly, the experimental processing should be improved to obtain better technical results. However, it should be highlighted that the objective of this work was not to explore the technical results, but to show that in course hours the students were able to overcome several difficulties and were successful in producing pure PLA and PLA/carbon fiber filaments by extrusion, 3D printing workpieces from these filaments and analyzing samples by tensile tests. Definitely, in a research exercise, the process should be improved with more regular filaments and understanding why the composite samples was less rigid. For further studies, the student also suggested to study of the effect of printing parameters on the adhesion between the layers and to evaluate the influence of the infill density on the produced composite 3D-printed parts. The Additive Manufacturing module responds perfectly to the needs of Industry 4.0, by training engineers capable of improving manufacturing performance thanks to their multidisciplinary skills in industrial engineering, materials sciences, online connectivity and data analysis, and by reinforcing their transversal skills such as intercultural communication, teamwork and creativity. Thanks to the pedagogical approach (problem-based learning) and the hybrid format, students will be directly confronted with problem solving in real situations in manufacturing and worked in an international multidisciplinary team.

Figure 1 :

 1 Figure 1: Equipment available exclusively for the group (Extruder, winder and filament stocking) and FDM printer. Processes Parameters Extrusion Barrel temperature = 187 o C Screw speed = 3 rpm Winder-extruder distance = 15 cm 3D printing Printing speed = 35 mm/s Nozzle temperature = 200 o C Bed temperature = 60 o C Infill density = 100% Raster angle = ± 45 o

Figure 2

 2 Figure2presents some images of the three experimental parts of their work: the filament production (with and without fibers), the specimen manufactured by FDM and the tensile tests.

Figure 2 :

 2 Figure 2: Manufacturing and Test Sequence.

Figure 4 :

 4 Figure 4: Final presentation and reports

https://www.univ-toulouse.fr/des-formations-pour-tous/idefi-defi-diversites

https://open.insa-toulouse.fr/

ACKNOWLEDGMENTS

The authors thank to the financial support of IDEFI DEFI Diversités 2020, the technical support of Philippe Seitier (Mechanical Engineering Department of INSA Toulouse) and to the International Relation Departments of INSA Toulouse (DRI) and Escola Politécnica (DARI-POLI-UFRJ).