Pierre Barbillon

Julien Chiquet

Tabouy Timothée

missSBM: An R Package for Handling Missing Values in the Stochastic Block Model

Keywords: Network, Missing data, Stochastic Block Model

The Stochastic Block Model (SBM) is a popular probabilistic model for random graphs. It is commonly used for clustering network data by aggregating nodes that share similar connectivity patterns into blocks. When fitting an SBM to a network which is partially observed, it is important to take into account the underlying process that generates the missing values, otherwise the inference may be biased. This paper introduces missSBM, an R-package fitting the SBM when the network is partially observed, i.e., the adjacency matrix contains not only 1's or 0's encoding presence or absence of edges but also NA's encoding missing information between pairs of nodes. This package implements a set of algorithms for fitting the binary SBM, possibly in the presence of external covariates, by performing variational inference adpated to several observation processes. Our implementation automatically explores different block numbers to select the most relevant model according to the Integrated Classification Likelihood (ICL) criterion. The ICL criterion can also help determine which observation process better corresponds to a given dataset. Finally, missSBM can be used to perform imputation of missing entries in the adjacency matrix. We illustrate the package on a network data set consisting of interactions between political blogs sampled during the French presidential election in 2007.

Introduction

In many fields of science, networks are a natural way to represent interaction data. To cite a few examples, a network may represent social interactions such as friendship or collaboration between people in a social network, regulation between genes and their products in a gene regulatory network, or predation between animals in a food web. In this paper, we only consider networks which can be represented by graphs composed of binary edges connecting pairs of nodes (also referred to as dyads in the following).

To this day, there exist many pieces of software performing network-related analyses. Un-arXiv:1906.12201v3 [stat.CO] 27 May 2021 surprisingly, the R community is extremely active in this area. Indeed, the R programming language is especially well-designed for performing data manipulation and visualization, and is thus appropriate for handling network data. Among the many available R packages related to networks, we suggest a classification into three groups1 : i) Packages for representation, manipulation or visualization tasks, and packages computing descriptive statistics. We mention non-exhaustively the following top representatives: igraph [START_REF] Csardi | The igraph software package for complex network research[END_REF], network and sna (Butts 2008a,b).

ii) Packages learning the structure of a network from an external source of data, such as huge [START_REF] Zhao | The huge package for highdimensional undirected graph estimation in R[END_REF], glasso [START_REF] Friedman | glasso: Graphical Lasso: Estimation of Gaussian Graphical Models[END_REF], bnlearn [START_REF] Scutari | Learning Bayesian networks with the bnlearn R package[END_REF] or bnstruct [START_REF] Franzin | bnstruct: an R package for Bayesian Network Structure Learning in the Presence of Missing Data[END_REF]. These packages generally rely on a specific graphical modeling of the data (e.g., Gaussian graphical models (Lauritzen 1996) in huge and glasso, or Bayesian networks [START_REF] Pearl | Bayesian networks[END_REF]) in bnlearn and bnstruct).

iii) Packages fitting (probabilistic) models on network data. The ergm package (Hunter, Handcock, Butts, Goodreau, and Morris 2008) fits the family of exponential random graph models (ERGM) introduced in Hunter and [START_REF] Hunter | Inference in curved exponential family models for networks[END_REF]: it is part of the collection of tools around ERGM regrouped in the statnet metapackage (Handcock, Hunter, Butts, Goodreau, and Morris 2008); latentnet [START_REF] Krivitsky | Fitting latent cluster models for networks with latentnet[END_REF] implements the latent space approach of [START_REF] Hoff | Latent space approaches to social network analysis[END_REF]; mixer (Ambroise, Grasseau, Hoebeke, Latouche, Miele, and Picard 2015) and blockmodels [START_REF] Léger | Blockmodels: A R-package for estimating in latent block model and stochastic block model, with various probability functions, with or without covariates[END_REF] fit the Stochastic Block Model (SBM) when the distribution of the edges belongs to the exponential family [START_REF] Snijders | Estimation and prediction for stochastic blockmodels for graphs with latent block structure[END_REF][START_REF] Nowicki | Estimation and Prediction for Stochastic Blockstructures[END_REF]. Other R packages related to the SBM and its extensions include sbm [START_REF] Chiquet | sbm: Stochastic Blockmodels. R package version 0.4.0[END_REF], sbmr [START_REF] Strayer | sbmr: Fit and investigate Stochastic Block Models in R. R package version 0.0.2[END_REF], dynSBM [START_REF] Matias | dynsbm: Dynamic Stochastic Block Models. R package version 0.2. Mejean A[END_REF], blockmodeling (Žiberna 2020), dBlockmodeling [START_REF] Brusco | dBlockmodeling: Deterministic Blockmodeling of Signed, One-Mode and Two-Mode Networks[END_REF], expSBM [START_REF] Rastelli | expSBM: An Exponential Stochastic Block Model for Interaction Lengths[END_REF], MLVSBM (Chabert-Liddell 2021), greed (Étienne Côme and Jouvin 2021), sbmSDP (Amini 2015), hergm [START_REF] Schweinberger | HERGM: Hierarchical exponential-family random graph models[END_REF], lda [START_REF] Chang | lda: Collapsed Gibbs Sampling Methods for Topic Models[END_REF], graphon [START_REF] You | graphon: A Collection of Graphon Estimation Methods[END_REF], GREMLINS [START_REF] Donnet | GREMLINS: Generalized Multipartite Networks. R package version 0[END_REF] and noisySBM [START_REF] Rebafka | noisySBM: Noisy Stochastic Block Mode: Graph Inference by Multiple Testing[END_REF]. Some of these packages, as well as some implementations in other programming languages, are presented in the following.

The missSBM package which we introduce here belongs to the third category, that is, software that fits a specific probabilistic model on network data. More specifically, missSBM is dedicated to the estimation of the Stochastic Block Model (SBM), a mixture of Erdős-Rényi random graphs [START_REF] Erdős | On random graphs[END_REF] offering a high degree of heterogeneity in connectivity profiles (see Abbe 2017, for a recent review). The SBM generally fits well real-world network data while keeping the advantage of being a probabilistic generative model (contrary to mechanistic approaches such as the Barabási-Albert model [START_REF] Albert | Statistical mechanics of complex networks[END_REF], defined by a preferential attachment algorithm). The main outcome of an SBM fit is a clustering of the nodes -or "blocks" -so that the nodes share the same properties within the same block. To our knowledge, the reference package for fitting the SBM with the R programming system is blockmodels. It includes efficient implementations of variational algorithms to fit different flavors of the SBM, adapted to binary network data and valued networks, with optional covariates on the edges. Two other important extensions of the SBM are available as R packages: the degree-corrected Stochastic Block Model in randnet [START_REF] Li | randnet: Random Network Model Selection and Parameter Tuning[END_REF] and a dynamic version of the Stochastic Block Model in dynsbm [START_REF] Matias | dynsbm: Dynamic Stochastic Block Models. R package version 0.2. Mejean A[END_REF].

Beyond the R framework, there also exist Python packages and C++ libraries providing efficient codes for some particular SBM: the Python packages CommunityDetection (Mejean and Maison 2017) and BipartiteSBM [START_REF] Yen | Blockmodeling on a Bipartite Network with Bipartite Prior[END_REF] are dedicated to the estimation of special network structures using various heuristics and network models, among which the SBM. Beyond variational approaches, MCMC methods exist for inferring the SBM, solving the exact problem but being generally more computationally demanding: the Python library graph-tool [START_REF] Peixoto | The graph-tool python library[END_REF]) includes an MCMC sampler to fit the binary SBM and its degreecorrected variant; C++ libraries sbm_canonical_mcmc (Young, Desrosiers, Hébert-Dufresne, Laurence, and Dubé 2017) and bipartiteSBM-MCMC (Yen and Larremore 2019) respectively implement a MCMC sampler for the SBM and the bipartite SBM. Finally MODE-NET (Decelle, Krzakala, and Zhang 2019) implements the belief propagation algorithm for inferring the degree-corrected SBM.

Despite their high quality, an important limitation of the aforementioned software is to require a network that is fully observed, that is, no missing value is supported. The main feature of missSBM is to deal with cases where the network data is only partially observed. More precisely, we consider situations where the adjacency matrix of the network data contains not only 1's or 0's for presence or absence of an edge, but also NA's encoding missing information for some dyads. Note that this situation is different from the case considered in noisySBM: there, a similarity matrix is fully observed between all pairs of nodes, and the goal is to separate the 'true' interactions from noise by means of a dedicated SBM.

When inferring the SBM from network data with missing values, it is important to take into account the underlying process that generates these missing values in the estimation of the model parameters, otherwise it may be biased. More specifically, one has to identify whether the values are Missing at Random or not (MAR and MNAR, see [START_REF] Little | Statistical analysis with missing data[END_REF]. This issue has been studied in the context of network data by [START_REF] Handcock | Modeling social networks from sampled data[END_REF] for the ERGM and in our methodological paper [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF] for the SBM. missSBM is an implementation of the methodology developed therein. It also considers new sampling designs and the inclusion of covariates simultaneously in the SBM and in the observation process, which was not studied by [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF]. Specifically, missSBM implements variational algorithms in the vein of [START_REF] Daudin | A mixture model for random graphs[END_REF] and [START_REF] Léger | Blockmodels: A R-package for estimating in latent block model and stochastic block model, with various probability functions, with or without covariates[END_REF] for estimating the SBM, with or without covariates, under various missing data mechanisms. This includes cases of incomplete data where the inference can be made only on the observed part of the data (MAR), or cases where it is necessary to take the sampling design into account in the inference (MNAR).

Some frameworks deal with missing data but rather from the cross-validation perspective than the sampling perspective. Cross-validation is used to perform model selection for networks such as the choice of the number of blocks or communities [START_REF] Li | Network cross-validation by edge sampling[END_REF][START_REF] Chen | Network cross-validation for determining the number of communities in network data[END_REF] or the choice of the latent structure [START_REF] Hoff | Modeling homophily and stochastic equivalence in symmetric relational data[END_REF]). Hence, these frameworks are quite different from ours since cross-validation is done under a MAR sampling while our main goal is to be able to infer an SBM under several MNAR sampling mechanisms.

The paper is organized as follows: Section 2 introduces the statistical framework of the binary SBM, with or without covariates, and summarizes the key points of its inference under missing data conditions. Section 3 provides basic user guidelines for the main functions and classes of objects. We finally detail in Section 4 a case study which analyzes a network data set describing the French blogosphere during the period preceding the 2007 French presidential election, illustrating the most striking features of the package.

Statistical Framework

Binary Stochastic Block Model (SBM)

In an SBM, nodes from a set N {1, . . . , n} are distributed among a set Q {1, . . . , Q} of hidden blocks which model the latent structure of the graph. The group membership is described by independent categorical variables (Z i , i ∈ N ) with multinomial distribution M(1, α = (α 1 , ..., α Q )). The probability of having an edge between any pair of nodes (or dyad) 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 of the other edges conditionally on the latent blocks:

Y ij |Z i , Z j ∼ ind B(π Z i Z j ), for all (i, j) ∈ D, (1)
where B stands for the Bernoulli distribution and D the set of dyads. This set may be either equal to {(i, j) ∈ N2 ; i = j} if the network is directed or to {(i, j) ∈ N 2 ; i < j}, otherwise 2 . In the following, we denote by π

= (π q ) (q, )∈Q 2 ∈ [0, 1] Q 2 the connectivity matrix, α ∈ D Q = {(α 1 , . . . , α Q ) ∈ [0, 1] Q ; α 1 + . . . α Q = 1} the block proportions, Z = (Z 1 , ..., Z n ) T the n × Q membership matrix and Y = (Y ij ) (i,j
)∈D the n × n adjacency matrix. This matrix is binary, with a diagonal filled with NA's and is symmetric if and only if the network is undirected. The vector encompassing all the unknown model parameters is θ = (α, π). A schematic representation of the binary SBM in the undirected case is given in Figure 1, where we highlight the latent clustering.

Accounting for External Covariates

On top of information about connections between nodes, it is common for network data to be accompanied with additional information on nodes or dyads, that we call covariates: for instance in social networks, nodes may belong to different categories (gender, occupation, nationality). Covariates on dyads usually represent similarity or dissimilarity between nodes: for example in a context of spatial data where nodes correspond to entities with explicit geographic locations, dyad covariates may be the distances between the nodes.

Depending on the analysis, we may want to detect a connectivity pattern beyond the covariate effect. To do so, we implemented in missSBM a variant of Model (1) for including covariates.

Let X ij ∈ R m denote the vector of m covariates for dyad (i, j). If the covariates correspond to the nodes, i.e., X i ∈ R N is associated with node i for all i ∈ N , they are transferred onto the dyad level through a symmetric "similarity" function φ(•,

•) : R N × R N → R m : X ij φ(X i , X j ). In the following, X [X ij ] i,j∈N ∈ (R m ) n×n denotes the array of covariates. A1 A2 A3 π •• B1 B2 B3 B4 B5 π •• C1 C2 π •• π •• π •• π •• • Q = {•, •, •} blocks • α • = P(i ∈ •), • ∈ Q, i = 1, . . . , n • π •• = P(Y ij = 1|i ∈ •, j ∈ •)
Figure 1: Schematic representation of an undirected network following the stochastic block model with 3 blocks. Colors are blocks in which nodes are dispatched with probabilities α and the distribution of the dyads depends on colors of nodes with probabilities π.

An SBM including the effect of these covariates is

Z i ∼ iid M(1, α), for all i ∈ N , Y ij | Z i , Z j , X ij ∼ ind B(g(γ z i z j + β X ij ), for all (i, j) ∈ D, (2)
where

γ q ∈ R, β ∈ R m , α = (α 1 , ..., α Q ), g(x) = (1 + e -x
) -1 . The vector of unknown parameters is now defined by θ = (γ, β, α). Note the connection with logistic regression: Model (2) assumes a logistic link between the presence of an edge Y ij and the corresponding covariates. The intercept term (γ q ) q depends on the blocks of the nodes and describes the heterogeneity of the connections that is not explained by the regression term β T X ij .

Connections with Similar Models. The SBM was originally introduced by [START_REF] Nowicki | Estimation and Prediction for Stochastic Blockstructures[END_REF]. Many extensions have been proposed since, including other distributions on dyads on top of the incorporation of covariates [START_REF] Mariadassou | Uncovering latent structure in valued graphs: A variational approach[END_REF]. Contrary to the latent space model of [START_REF] Hoff | Latent space approaches to social network analysis[END_REF] where the latent space is continuous, the latent variables in the SBM lie in a discrete space. When covariates are included as in Equation ( 2), their effect is removed and the blocks shall then explain the structure in the network beyond the covariates. This approach is similar to [START_REF] Vu | Model-based clustering of large networks[END_REF] and opposite to the one used by [START_REF] Tallberg | A Bayesian approach to modeling stochastic blockstructures with covariates[END_REF] or [START_REF] Binkiewicz | Covariate-assisted spectral clustering[END_REF] where the covariates help learn the underlying clustering of the nodes since their distribution is assumed to depend on the same latent variables as the SBM.

Missing Data and SBM

The main feature of missSBM is to deal with some processes that generate missing values in order to provide more accurate estimates of the parameters underlying an incompletely observed network. The number of nodes n is assumed to be known and the missing information only concerns the dyads. Hence the sampled data can be encoded in an adjacency matrix Y where missing information -dyads whose value is unobserved -is encoded by NA's. We also define the n × n observation matrix R such as R ij = 1 if the dyad Y ij is observed and R ij = 0 otherwise. For convenience we define Y o = {Y ij : R ij = 1} and Y m = {Y ij : R ij = 0} the respective sets of observed and unobserved dyads.

In our framework, an observation process -or sampling design -is a stochastic process that generates R. We then rely on the standard missing data theory of [START_REF] Little | Statistical analysis with missing data[END_REF] to classify those designs either into Missing Completely At Random (MCAR), Missing At Random (MAR) or Missing Not At Random (MNAR) cases. This framework has to be extended to handle the latent variables Z in the SBM as we did in [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF]:

Sampling design for SBM is

       MCAR if R |= (Y, Z) | X, MAR if R |= (Y m , Z) | (Y o , X), MNAR otherwise.
(3)

The notation |= stands for independence between random variables. Note that MCAR missingness is a particular case of MAR missingness. This definition provides the general case when covariates X are available. Otherwise, the definition remains valid just by removing X. Denoting by ψ the set of parameters associated with the distribution that generates the sampling matrix R, we assume that ψ and θ live in a product space, so that we can derive the following proposition:

Proposition 1. From (3), if the sampling is MAR or MCAR then maximizing p θ,ψ (Y o , R, X) or p θ (Y o , X) in θ is equivalent for θ.
This proposition was proven in [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF] in the absence of covariate. The generalization to handle covariates is straightforward. In words, the inference can be conducted on the observed part of the network data when the sampling is M(C)AR without incurring any bias. In these cases, adaptating the existing algorithms for SBM inference is simple. MNAR sampling designs require, however, more refined inference strategies since the observation process has to be included in the inference.

Examples of Sampling Designs for Networks

This section reviews a set of stochastic processes defining sampling designs available in missSBM. These sampling designs may depend either on i) the values of the dyads in the network; ii) the latent clustering of the nodes; or iii) some covariates, via the vector of parameters ψ. All examples detailed in the following assume that the observations are independent conditionally on Y, Z and X (either observations of the dyads for dyad-centered sampling designs, or observations of the nodes for node-centered sampling designs). From a practical viewpoint, the sampling designs implemented in missSBM allow the user to either i) generate new data by partially observing an existing network according to a predefined sampling design, with user-defined parameters ψ, or ii) to fit an SBM model under missing data condition, assuming that the missing entries arise from a given type of sampling design for which the unknown parameters ψ are estimated jointly with the SBM parameters θ.

Dyad-Centered Sampling Designs

• Dyad Sampling (MCAR): each dyad (i, j) ∈ D has the same probability P(R ij = 1) ψ ∈ [0, 1] to be observed.

• Double Standard Sampling (MNAR): let ψ (ρ 1 , ρ 0 ) ∈ [0, 1] 2 . 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 .
The probability for sampling a dyad thus intrinsically depends on the presence/absence of the corresponding edge. This double standard sampling is especially likely in real world applications, if it is easier to observe an existing connection than the absence of connection. For instance, in protein-protein networks, the sampling effort is more important to determine the absence of a link than its existence.

• Block-Dyad Sampling (MNAR): this sampling consists in observing all dyads with probabilities ψ (ψ q ) (q, )∈Q 2 ∈ [0, 1] Q 2 depending on the underlying clustering of the network:

ψ q = P(R ij = 1 | Z iq = 1, Z j = 1). • Covar-Dyad Sampling (MAR): let us define ψ (α, κ) ∈ R × R m .
Here the probability for observing a dyad is driven by the effect of a given covariate:

P(R ij = 1|X ij ) = g(α + κ t X ij ),
where we recall that g(x) = (1 + e -x ) -1 . Under this sampling, the external covariates may impact both a connection and the ability to observe it. In this case, the sampling remains MAR provided that the covariates are available.

Node-Centered Sampling Designs

A node-centered sampling consists in observing some nodes sampled with probabilities given by the sampling design. Observing a node means observing all the dyads involving that node. For all i ∈ N , we denote by V i the indicator variable for observing node i. Hence if V i = 1 we have R ij = R ji = 1 for all j ∈ N . Node-centered sampling designs are likely in social sciences since a network is sampled through direct interviews. During an interview, individuals (nodes) indicate to whom they are connected. Some individuals may then indicate a connection with an individual not available for an interview. The resulting missing dyads concern dyads between individuals who were not interviewed. Even if the connection is oriented (directed network), we assume that an individual, when interviewed, provides its ingoing and outgoing connections.

• Node Sampling (MCAR): the probabilities for observing nodes are uniform:

P(V i = 1) = ψ ∈ [0, 1] for all i ∈ N .
• Snowball sampling (MAR): a first batch of nodes is sampled as in node sampling. Then, a second batch is composed of the neighbors of the first batch (the set of nodes linked to at least a node of the first batch). Other batches can then be obtained through several sampling steps which are called waves. These successive waves are then MAR and not MCAR since they are built on the basis of the previously observed part of Y .

• Degree Sampling (MNAR): for all node i ∈ N , P(

V i = 1) = ρ i where (ρ 1 , . . . , ρ n ) ∈ [0, 1] n are such that ρ i = g(a+bD i ) for all i ∈ N where ψ (a, b) ∈ R 2 and D i = j Y ij .
This sampling may be the consequence of a situation where popular individuals are more likely to be interviewed.

• Block-Node Sampling (MNAR): this sampling consists in observing all dyads corresponding to nodes selected with probabilities ψ

(ψ 1 , . . . , ψ Q ) ∈ [0, 1] Q such that ψ q = P(V i = 1 | Z iq = 1
) for all (i, q) ∈ N × Q. This sampling may happen if some communities that shape the connections in the network are not equally reachable.

• Covar-Node Sampling (MAR): let ψ (ν, η) ∈ R × R N . The probability to observe a node is

P(V i = 1|X i ) = g(ν + η t X i ).
In this sampling, some external information shapes the sampling process of the nodes. Even if the covariates also impact the probabilities of connection, as in the Covar-dyad sampling, the sampling design is MAR provided that the covariates are available.

Estimation Procedure: a Variational EM

The SBM is a latent state space model which can be seen as a mixture model for random graphs. Therefore, the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is the natural choice for the inference since it generally proves very useful for inferring various types of mixture models. It is based on the evaluation of the expectation of the complete log-likelihood of the model, with respect to the conditional distribution of the latent variables given the data. However, this expectation is intractable in the SBM due to the structure of dependency between the latent variables Z and the network Y. In fact, it would require to sum over all possible clusterings for all pairs of nodes, which is out of reach even for a moderate number of nodes or blocks. To address this shortcoming when the network Y is fully observed, [START_REF] Daudin | A mixture model for random graphs[END_REF] introduced a variational EM, based on the variational principles of [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF]. The idea is to maximize a lower bound of the loglikelihood based on an approximation of the true conditional distribution of the latent variable Z. In the case of an SBM with missing data, the level of difficulty is higher since the set of latent variables encompasses both Z (the latent blocks) and Y m (the missing dyads). We propose here a variational distribution of the conditional distribution p θ,ψ (Z, Y m |Y o ) where complete independence is forced on Z and Y m , using a multinomial, respectively a Bernoulli distribution for Z and Y m . We denote by m(•) and b(•) the probability density functions of, respectively, the multinomial and the Bernoulli distributions which gives the following expression of the variational distribution:

pτ,ν (Z, Y m ) = pτ (Z) pν (Y m ) = i∈N m(Z i ; τ i ) (i,j)∈D m b(Y ij ; ν ij ),
where

τ = {τ i = (τ i1 , . . . , τ iQ ) ∈ [0, 1] Q : Q q=1 τ iq = 1, i ∈ N } and ν = {ν ij ∈ [0, 1], (i, j) ∈ D m }
are the two sets of variational parameters respectively associated with Z and Y m . Interestingly, τ 's are proxies for the posterior probabilities of the group memberships for all nodes, and ν's correspond to the imputed values of the missing dyads in the network data. This variational distribution was chosen in order to replace the intractable E-step of the EM algorithm with a tractable variational E-step. This approximation leads to the following lower bound J of the log-likelihood, where KL is the Kullback-Leibler divergence between the true conditional distribution and its variational approximation:

log p θ,ψ (Y o , R) ≥ J τ ,ν,θ,ψ (Y o , R) log p θ,ψ (Y o , R) -KL pτ,ν (Z, Y m ) p θ (Z, Y m Y o ) , = E pτ ,ν [log p θ,ψ (Y o , R, Y m , Z)] -E pτ ,ν [log pτ,ν (R, Y m )] .
If we choose p = p θ,ψ (Z, Y m |Y o ), the true conditional distribution of the latent variables Z, Y m , we retrieve the standard EM algorithm, requiring the evaluation of the intractable quantity

E p θ,ψ (Z,Y m |Y o ) [log p θ,ψ (Y o , R, Y m , Z)].
Note that we alleviated the notations above by not explicitly writing the possible conditioning on covariates in the log-likelihoods.

Based on this approximation, the variational EM algorithm consists in alternating updates of the variational parameters {τ , ν} (the VE-step) with updates of the model parameters θ, ψ (the M-step) maximizing J.

Steps VE and M are iterated until convergence like in a standard EM. The algorithm converges to a local maximum of the lower bound of the log-likelihood. This variational is translated into a collection of algorithms for handling missing data with all sampling designs introduced in Section 2.4. When an algorithm reaches convergence, we obtain estimates of the parameters involved in the SBM (θ), in the sampling process (ψ), and also estimates of the variational parameters which bring information on the clustering (τ ) and on the missing dyads (ν). The variational estimator in the SBM is proven to be asymptotically normal in [START_REF] Bickel | Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels[END_REF] when the network is fully observed. The extension to the MCAR case is proven in [START_REF] Mariadassou | Consistency and asymptotic normality of stochastic block models estimators from sampled data[END_REF].

Initialization. It is well known that EM-like algorithms have a great sensitivity to the initialization step, which therefore requires a special attention. In missSBM, the initial clustering is obtained by applying the popular Absolute Eigenvalues Spectral Clustering (detailed in Rohe, Chatterjee, and Yu 2011) to the adjacency matrix where the NA's are replaced by zero.

The initial clustering can also be provided by the user. As specified in the next paragraph, the exploration of the number of blocks also provides several other relevant initializations.

Selection of the Number of Blocks.

A main difficulty met when conducting SBM inference lies in the estimation of the number of blocks, generally unknown to the user. To remedy this problem, we use the Integrated Classification Likelihood (ICL) criterion defined in [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] and routinely used in the framework of mixture models. Note that Saldana, Yu, and Feng ( 2017

ICL(Q) = -2E pτ ,ν log p θ, ψ(Y o , Y m , R, Z|Q, K) + pen ICL (Q, K), 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.

We also implemented in missSBM an exploration procedure designed to avoid getting stuck in local minima by producing a convex and robust ICL curve. This exploration procedure is divided into two steps, forward and backward: the forward step creates new initializations for each number Q of block considered, by splitting blocks obtained from estimations with Q -1 blocks. On the other hand, the backward step tries new initializations for each Q by merging groups of the model with Q + 1 groups. The best model in terms of ICL is always retained. The procedure can be iterated until a satisfying shape of the ICL curve is met.

Implementation Details. missSBM adopts an oriented-object programming spirit for representing most models by means of R6-classes and the R6 package of [START_REF] Chang | R6: Classes with Reference Semantics[END_REF], which paves the way for future extensions. This approach is non-visible to the user, who essentially only has to deal with classical R functions. The time consuming pieces of code are written in C++ using the armadillo library for linear algebra [START_REF] Sanderson | Armadillo: a template-based C++ library for linear algebra[END_REF], in conjunction with the Rcpp and RcppArmadillo packages [START_REF] Eddelbuettel | Rcpp: Seamless R and C++ integration[END_REF][START_REF] Eddelbuettel | RcppArmadillo: Accelerating R with high-performance C++ linear algebra[END_REF]) to interface C++ with R. The numerical performance of our implementation is of the same order as existing variational implementations of the binary SBM (like blockmodels). It allows to deal with networks with up to a couple of thousands of nodes. To give an idea to the reader, Figure 2 reports the timings for adjusting an SBM to the network data considered in Section 4 (the 2007 French political blogosphere, 194 nodes), for a varying number of blocks with a single Intel Core i9-9900 CPU at 3.10GHz, replicated 50 times per block-size. We report in the discussion some interesting tracks for improving the speed in the future and eventually tackle larger networks.

Guidelines for Users

This section gives an overview of the basic usage of the package. In particular, it describes inputs and outputs of the main functions at play in the procedure that fits an SBM from partially observed network data. Along the section, we also describe the classes of objects included in the package to ease the manipulation of the resulting fitted models. As a complement to this section and to the usual R package documentation, a full documentation including a vignette is available as a pkgdown website at URL http://grosssbm.github.io/missSBM.

Parameter Estimation, Prediction and Clustering

Estimation of an SBM from a partially observed network is done by means of the function estimateMissSBM, with the following usage: missSBM:::available_samplings ## [1] "dyad" "covar-dyad" "node" "covar-node" ## [5] "block-node" "block-dyad" "double-standard" "degree" ## [9] "snowball" Argument covariates is an optional list with as many entries as covariates. If the covariates are nodal, each entry must be a size-n vector; if the covariates are defined between pairs of nodes, each entry must be an n × n matrix.

Advanced Tuning. The argument control is a list to control the estimation and finely tune the variational EM algorithm, with the following entries: i) useCov: Logical indicating whether the covariates should be incorporated within the SBM (or just in the sampling);

ii) clusterInit: Initial method for clustering. Either a character ("spectral") or a list with length(vBlocks) vectors, each with size ncol(adjacencyMatrix), providing a userdefined clustering. Default is "spectral", for absolute spectral clustering;

iii) similarity: An R × R → R function to compute similarities between nodal covariates (default is missSBM:::l1_similarity, that is, (x, y) → -|x -y|); iv) threshold: Optimization stops when a V-EM step changes the objective function or the connection parameters by less than threshold (default is 1.10 -2 ); v) maxIter: Optimization stops when the number of iteration exceeds maxIter (default is 50); vi) fixPointIter: Number of iterations in the fixed-point algorithm used to solve the variational E step (default is 3); vii) exploration: Character indicating what kind of exploration should be used among "forward", "backward", "both" or "none" (default is "both"); viii) iterates: Integer for the number of iterations during the exploration process. Only relevant when exploration is different from "none" (default is 1); ix) trace: Logical for verbosity (default is TRUE).

Return Value: Classes missSBM_collection and missSBM_fit. An output of the function estimateMissSBM is an instance of an R6 object with class missSBM_collection, which can be handled as a standard list. Its fields are described in Table 1.

Field Description models

a list of models with class missSBM_fit ICL the vector of ICL associated to the models in the collection bestModel best model according to the ICL (a missSBM_fit object) optimizationStatus a data.frame summarizing the optimization process for all models Method Description plot(object, type) plot either "icl", "elbo" or "monitoring" Table 1: Structure of missSBM_collection (output of estimateMissSBM).

Among the fields of missSBM_collection, models is a list with as many elements as in vBlocks. These elements are R6 objects with class missSBM_fit, the fields of which are detailed in Table 2. It gives access to the results of the inference for a fixed number of blocks.

We finally give additional details on fields fittedSBM and fittedSampling in Table 3. Note that fittedSBM enjoys some standard plot, print, coef, predict and fitted methods, most of them inherited from the class simpleSBM_fit of the package sbm.

Models Exploration.

At the end of the estimation process, it is common that the algorithm gets stuck in some local minima for some values of Q, the current number of blocks. The consequence is "non-smooth" ICL curve which is theoretically convex and rather smooth. This may lead the user to choose a sub-optimal number of blocks. Thus, the ICL criterion is automatically "smoothed" after a first pass across all the models. The idea is to apply a split and merge strategy to the path of models stored in missSBM_collection in order to find a better initialization for each value of Q in the V-EM algorithm, so that it converges to the global minimum. This model exploration can be tuned with the control argument, see paragraph Advanced tuning. The default goes back and forth a single time. Method Description plot(object, type) plot either "imputed", "expected", "meso" or "monitoring" print(object) a print method recalling important fields and methods available coef(object, type) model parameters (type="mixture", "connectivity", "covariates" or "sampling" ) predict(object) estimated adjacency matrix (with imputation)

fitted(object)

expected value of the SBM Parallel Computing. Some internal components of estimateMissSBM (initialization, estimation, exploration) use the future framework [START_REF] Bengtsson | A Unifying Framework for Parallel and Distributed Processing in R using Futures[END_REF] to speed-up the whole process with parallel computing. The future::plan must be set by the user (default is sequential with no parallelism). For instance, in order to run missSBM on 10 cores with forking on Unix systems, the following command should be used before calling estimateMissSBM().

future::plan("multicore", workers = 10)

Partial Observation ("sampling") of Some Network Data

The function observeNetwork generates missing entries in a fully observed adjacency matrix according to a given network sampling design. The usage is the following:

observeNetwork( # the original adjacency matrix of the network to be partially observed adjacencyMatrix, # a character for the sampling design which defines the observation process sampling, # a set of parameters associated to the sampling design parameters, # an optional clustering membership vector required for block samplings clusters = NULL, # an optional list of covariates that may influence the observation process covariates = list(), # an optional function to compute similarities between node covariates similarity = l1_similarity, # an optional intercept term added in case of the presence of covariates intercept = 0)

Note that the dimension (or length) of parameters depends on the sampling design selected, as described in Section 2.4. Argument clusters only needs to be specified for "block-dyad" and "block-node" sampling designs. Arguments covariates is by-default list(), covarSimilarity is set to the l1_similarity function defined by (x,

y) ∈ R d ×R d → -|x -y| 1 ∈ R d
, and finally intercept is set to 0. These last three arguments only need to be specified in the case of an SBM with covariate(s). Note that the intercept is not included in parameters and must be specified independently. The output of observeNetwork is a matrix encoding the adjacency-matrix, with NA values for dyads not observed by the observation process, which can be provided as input to estimateMissSBM.

Illustration: the 2007 French political blogosphere

This section illustrates the features of missSBM by conducting the analysis of a real-world network data. It is a sub-network of the French political blogosphere, extracted from a snapshot of over 1100 blogs collected during a period preceding the 2007 French presidential election, and manually classified by the "Observatoire Présidentielle project" (see [START_REF] Zanghi | Fast online graph clustering via Erdős-Rényi mixture[END_REF]. The network is composed of 194 blogs representing nodes in the network and 1432 edges indicating that at least one of the two blogs references the other. On top of missSBM, our analysis relies on aricode [START_REF] Chiquet | aricode: Efficient Computations of Standard Clustering Comparison Measures[END_REF] for computing clustering comparison measures and tidyverse for performing data manipulations and producing graphical outputs. The igraph package, imported by missSBM, is needed for basic graph manipulations. The future package is used for parallel computing. We also fix the seed for reproducibility: library(missSBM) library(aricode) library(tidyverse); theme_set(theme_bw()) library(igraph) library(future) set.seed(03052008)

We set our future::plan to multicore with 10 workers3 . future::plan("multicore", workers = 10)

The frenchblog2007 data set is provided with missSBM4 as an igraph object. We extract the adjacency matrix corresponding to the blog network after removing the isolated nodes.

We also extract the political party of each blog from the vertex attribute party, which gives us a natural classification of the nodes that could be used as a reference in our analyses.

data("frenchblog2007", package = "missSBM") frenchblog2007 <-delete_vertices(frenchblog2007, which(degree(frenchblog2007) == 0)) blog <-as_adj(frenchblog2007) party <vertex.attributes(frenchblog2007)$party For this network, we explore models with a number of blocks varying from 1 to 18: blocks <-1:18 Standard SBM Estimation. At this stage, the data set has no missing entry: every dyad and every node is observed. The adjacency matrix Y of the fully-observed network is stored in the variable blog. We first perform a standard SBM estimation on the fully observed network, including smoothing of the ICL.

sbm_full <-estimateMissSBM(blog, blocks, "node")

We inspect the result of the optimization process by plotting the ELBO (variational lower bound of the log-likelihood) against the number of iterations in the V-EM algorithm, accumulated along all the numbers of blocks considered (Figure 3). The ELBO is typically expected to increase with the number of blocks at the end of the successive optimizations, which is indeed the case here: plot(sbm_full, type = "monitoring") or plot diverse outputs, for instance a matrix view of the original network with columns and rows reordered according to the block memberships of the nodes, see Figure 4).

plot(sbm_full$bestModel, dimLabels = list(row = "blogs", col = "blogs")) For illustrative purposes, we sample a sub-graph from the blog network to mimic missing data and create a new adjacency matrix with missing entries. Since data consists in interactions between blogs (who references who), it is natural to sample the network with a node-centered sampling design: the following code generates a realization of a 194 × 194 sampling matrix according to the block-node sampling, where the blocks correspond to the clustering estimated by the SBM fitted on the full data set. The sampling rate is either low (0.2) in small blocks or high (0.8) in large blocks.

samplingParameters <ifelse(sbm_full$bestModel$fittedSBM$blockProp < .1, 0.2, .8) blog_obs <-observeNetwork( adjacencyMatrix = blog, sampling = "block-node", parameters = samplingParameters, clusters = sbm_full$bestModel$fittedSBM$memberships ) Estimation of a Partially Observed Network. We now perform SBM inference under missing data conditions by fitting two types of model: first, the SBM under the MNAR blocknode sampling design, i.e., under the design that truly generated the missing entries; second, the SBM under the MCAR node sampling design which basically performs inference only on the observed part of the network, neglecting the process that generates the missing values. The estimation is run on both models with the same settings as for the fully observed data. The ICL curve is smoothed with five iterations of forward-backward exploration since the presence of missing values typically increases the chances for falling into local minima. sbm_block <-estimateMissSBM(blog_obs, blocks, "block-node", control = list(iterates=5)) sbm_node <-estimateMissSBM(blog_obs, blocks, "node", control = list(iterates=5))

Sampling Design Comparison. We plot on Figure 5 the ICL of the three models ("fully observed", "block-node" sampling and "node" sampling) to show how it can be compared to select which sampling design fits at best the data: rbind( tibble(Q = blocks, ICL = sbm_node$ICL , sampling ="node"), tibble(Q = blocks, ICL = sbm_block$ICL, sampling = "block-node"), tibble(Q = blocks, ICL = sbm_full$ICL , sampling = "fully observed") ) %>% ggplot(aes(x = Q, y = ICL, color = sampling)) + geom_line() + geom_point() + ggtitle("Model Selection") + labs(x = "# blocks", y = "Integrated Classification Likelihood")

Note that the curves associated with the block-node sampling and the node sampling are quite close for small numbers of blocks (less than 10) and then depart from each other: the choice of the sampling design is a tough question for the network data at play. We also represent the ICL curve for the collection of SBM estimated on the fully observed network: although the values of ICL cannot be compared between this model and the two obtained on the partially observed network (data sets are not the same), we underline that the numbers of blocks selected in the different cases remain comparable. A model with 10 blocks is selected with the fully observed network, while a model with only 9 blocks is chosen for the block-node sampling SBM. Indeed, due to the partial sampling, some blocks are less well represented than others, and it seems more likely to gather some blocks together considering the information available. Regarding the clustering obtained by the three variants, we compare them with the Adjusted Rand Index (ARI, Rand 1971) computed with the aricode package [START_REF] Chiquet | aricode: Efficient Computations of Standard Clustering Comparison Measures[END_REF]. We use the classification of the fully-observed SBM as a reference, since its clustering was used to sample the network with the block-node sampling design. We typically expect Figure 5: ICL for models with fully observed data, block-node sampling and node sampling that a model relying on a better modeling of the missing values shall lead to a clustering closer to the reference. This is indeed the case when looking at the next piece of code, where it is shown that the ARI with the reference clustering is higher for the block-sampling SBM than for the MCAR node sampling SBM: ARI(sbm_block$bestModel$fittedSBM$memberships, sbm_full$bestModel$fittedSBM$memberships)

[1] 0.6570555 ARI(sbm_node$bestModel$fittedSBM$memberships , sbm_full$bestModel$fittedSBM$memberships)

[1] 0.5436008 With myModel, we now have at hand a tool for analyzing the clustering of the French political blogosphere. The first output is the connectivity matrix of the network, which puts into light the community structure of the blogosphere. Indeed, it is revealed by a diagonal filled with high probabilities and off-diagonal with low probabilities. Thus, nodes (blogs) in blocks connect with a high probability with other nodes in the same block and with a low probability with nodes in other blocks. Such a network concentrates most of its edges between nodes of the same blocks. This can be seen by displaying the probability of connection predicted by the SBM at the whole network scale (see Figure 6): plot(myModel, type = "expected", dimLabels = list(row="blogs", col="blogs"))

For validation, we suggest comparing the clustering of the model with the node attribute corresponding to the political parties to which the blogs belong. First, we remark that the SBM fitted on missing entries carries a little bit less information regarding the political party than the SBM adjusted on the fully observed network: ARI(party, myModel$fittedSBM$memberships)

[1] 0.4126328 ARI(party, sbm_full$bestModel$fittedSBM$memberships)

[1] 0.463709 A more detailed comparison between blocks inferred by the SBM and political parties is reported in Figure 7 with an alluvial diagram.

Also remember that missSBM performs imputation of the missing dyads in the adjacency matrix. Thus, we can compare the imputed values with the values of the dyad in the fully observed network to validate the performance of our approach. Using the R package pROC [START_REF] Robin | pROC: an Open-Source Package for R and S+ to Analyze and Compare ROC Curves[END_REF], we check the quality of the imputation. We replicate this experiment 500 times with a sampling rate varying between ≈ 0.4 and ≈ 0.9 (always with block-sampling design), fixing the number of blocks to the best one found on the fully observed network. The Area Under the Curve (AUC) is plotted in Figure 8 against the sampling rate, showing the robustness and the good performance of the imputation method. Dealing with Covariates. In order to illustrate how we can deal with covariates in the observation process, we now consider a sampling which depends on the political parties of the blog. For illustrative purposes, we extract a sub-graph that only contains the nodes whose political party is either left or right: blog_subgraph <-frenchblog2007 %>% igraph::induced_subgraph(V(frenchblog2007)$party %in% c( "right", "left")) blog_subgraph <-delete_vertices(blog_subgraph, which(degree(blog_subgraph) == 0))

The sub-graph is given in Figure 9, generated with the following piece of code: We then build a simple binary covariate on nodes indicating its party (1/0 = left/right).

dummy_party <-(V(blog_subgraph)$party == "left") * 1

Now, the observation process of the network is assumed to depend on this covariate so that a node belonging to the left party is more likely to be observed than a node belonging to the right party:

blog_subgraph_obs <-blog_subgraph %>% as_adj() %>% observeNetwork(sampling="covar-node", parameters = 10, covariates = list(dummy_party)) blocks <-2:8 future::plan("multicore", workers = 10) Since the covariate is nodal, a nodal observation process is considered. When it comes to take into account the effect of this covariate on the probability of connection between two nodes in the SBM as in Equation ( 2), the covariate has to be transferred at the dyad level. This consists in computing an n × n matrix whose elements are equal to one if the corresponding two nodes belong to the same political party, zero otherwise. This matrix computation is directly handled in missSBM, when the option useCov is TRUE.

From the observed network blog_subgraph_obs, four modeling choices (labeled i to iv) are possible whether the covariate is taken into account in the SBM or not, and whether the sampling is set as depending on the covariate or not. For the latter choice, we know that the sampling which depends on the covariate is more appropriate but this information is generally not available a priori. We then run the estimation in these four scenarios below and print the estimated parameters related to the SBM and the sampling. i) Take the covariate into account in both the sampling and the SBM: For models i) and ii), we notice that the estimates of the sampling parameters are quite accurate. For models iii) and iv), the estimation of the unique sampling parameter boils down to the proportion of observed nodes. For models i) and iii), the estimate of the parameter β in Equation ( 2) is positive and quite high. Thus, we conclude that the probability of connection between two nodes is strengthened by the fact that the nodes belong to the same party.

Figure 10 shows clearly that the "covar-node" sampling is uncovered from the data since the ICL criterion favors this sampling no matter if the covariate is taken into account in the SBM or not. The ICL criterion also points out that the models including the covariate effect in the SBM are better. rbind( tibble(Q=blocks, ICL=sbm_covar1$ICL, sampling="covar-node", useCov= true ), tibble(Q=blocks, ICL=sbm_covar2$ICL, sampling="covar-node", useCov= false ), tibble(Q=blocks, ICL=sbm_covar3$ICL, sampling="node" , useCov= true ), tibble(Q=blocks, ICL=sbm_covar4$ICL, sampling="node" , useCov= false ) ) %>% ggplot(aes(x = Q, y = ICL, color = sampling, shape = useCov)) + geom_line() + geom_point() + labs(x = "#blocks", y = "ICLs")

Finally, we compare the clustering obtained with the four models together by computing pairwise ARIs. We also compare them with the clustering obtained by the SBM fitted on the fully observed network, adjusted as follows: sbm_covar_full <-as_adj(blog_subgraph) %>% estimateMissSBM(blocks, "node", covariates = list(dummy_party))

For the sake of clarity, the ARIs are displayed in Table 4. We obtain the same clustering with model i) and iii), which is the closest to the one obtained with the fully observed network. This is expected since these models also take into account the effect of the covariate.

Discussion

The R package missSBM enables the estimation of an SBM from partially observed binary networks even for some observation processes which generate MNAR data.

Although version 1.0.0 of missSBM only deals with binary networks, we deploy a structure that lets the possibility to easily include other variants in the future. In particular, extending missSBM to weighted SBM where the distribution on the edges belongs to the exponential family (e.g. Poisson distribution or Gaussian distribution, see [START_REF] Mariadassou | Uncovering latent structure in valued graphs: A variational approach[END_REF]) should be straightforward in the MAR case. To pave the way of such a generalization, missSBM relies on the package sbm which is designed to offer a collection of methods and algorithms to manipulate more general SBM, yet in the absence of missing data, or at least, with imputed data. Hence, missSBM could take advantage in the future of improvements and new advances in sbm, while focusing on the handling of missing data specifically, by dealing only with the modeling of the observation process and the imputation of the missing entries. The modular object-oriented coding of missSBM also allows the developer to make it easily evolving in the way it can take into account the missing data: new observation processes could be incorporated by providing new sampling design which should contain the functions corresponding to the sampling specific steps in the VEM algorithm. Other dependency structures between the SBM and the covariates could also be modeled and integrated. For example, the latent variables Z could directly depend on the covariates X instead of the adjacency matrix Y.

Some recent work is concerned with boosting the speed of the variational inference (see e.g. [START_REF] Blei | Variational inference: A review for statisticians[END_REF]. This work could be adapted in the SBM context with or without missing data in order to enable our package to deal with larger networks. Resorting to minorization-maximization algorithms within the Variational E-step as in [START_REF] Vu | Model-based clustering of large networks[END_REF] could also be interesting for speeding up the algorithm. Finally, a common drawback of the variational inference is to provide too narrow standard errors for the estimated parameters [START_REF] Westling | Beyond prediction: A framework for inference with variational approximations in mixture models[END_REF]. Moreover, there is no uncertainty measure on the stability of the node clustering. We consider it as future work to provide the user with confidence intervals and stability measures of the clustering based on resampling techniques such as the bootstrap.

  ),[START_REF] Wang | Likelihood-based model selection for stochastic block models[END_REF],[START_REF] Hu | Using maximum entry-wise deviation to test the goodness of fit for stochastic block models[END_REF] and[START_REF] Côme | Hierarchical clustering with discrete latent variable models and the integrated classification likelihood[END_REF] propose alternative methods for selecting the number of blocks. The ICL criterion was adapted to the SBM under missing data condition in[START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF], and is recalled here: for a model with Q blocks, a sampling design parametrized by K parameters (size of ψ) and ( θ, ψ) = arg max (θ,ψ) log p θ,ψ (Y o , Y m , R, Z), then

Figure 2 :

 2 Figure 2: Timings for adjusting binary SBM with missSBM on the French political blogosphere with a single core for a varying number of blocks (50 replicated runs per # blocks).

  network data: matrix with NA entries vBlocks , # vector of block numbers for model exploration sampling = "dyad", # string describing the sampling design covariates = list(), # optional list of covariates (dyad or nodal) control = list()) # optional list for tuning the algorithm Standard Usage. This function takes as a first argument a square base::matrix or a sparsely encoded Matrix::dgCMatrix, possibly with NA entries corresponding to missing (unobserved) values of the network. The second argument vBlocks contains the successive explored values for the number of blocks, this number being generally unknown. The third argument -sampling -specifies how NA entries are taken into account in the estimation. It must be one of the following character strings, corresponding to one of the sampling designs (or observation processes) depicted in Section 2.4:

  variational bound J τ ,θ at each step double vExpec value of Ep τ [log(p θ (Y, Z))] at each step double penalty penalty of the model with Q blocks double

Figure 3 :

 3 Figure 3: Evolution of the ELBO during the optimization for the successive numbers of blocks

Figure 4 :

 4 Figure 4: Network data reorganized by the estimated block memberships

Figure 6 :Figure 7 :

 67 Figure 6: Probabilities of connection predicted by the SBM with block-node sampling

  $bestModel$fittedSBM$memberships nBlocks <-sbm_full$bestModel$fittedSBM$nbBlocks future::plan("sequential") res_auc <mclapply(1:500, function(i) { subGraph <-observeNetwork(blog, "block-node", runif(nBlocks), cl0) missing <which(as.matrix(is.na(subGraph))) true_dyads <-blog[missing] sbm_block <-estimateMissSBM(subGraph, nBlocks, "block-node", control = list(cores = 1, trace = 0)) imputed_dyads <-sbm_block$bestModel$imputedNetwork[missing] c(rate = 1length(missing)/length(blog), auc = auc(true_dyads, imputed_dyads, quiet = TRUE)) }, mc.cores = 10) purrr::reduce(res_auc, rbind) %>% as.data.frame() %>% ggplot() + aes(x = rate, y = auc) + geom_point(size = 0.25) + geom_smooth(method = "lm", formula = y ~x, se = FALSE) + labs(x = "sampling rate", y = "Area under the ROC curve")

Figure 8 :Figure 9 :

 89 Figure 8: Area Under the Curve (AUC) of the imputation as a function of the sampling rate.

  sbm_covar1 <-estimateMissSBM(blog_subgraph_obs, blocks, "covar-node", covariates = list(dummy_party), control = list(useCov = TRUE, iterates = 2)) sbm_covar1$bestModel$fittedSampling$parameters # sampling parameters sbm_covar1$bestModel$fittedSBM$covarParam # the covariate into account in the sampling only: sbm_covar2 <-estimateMissSBM(blog_subgraph_obs, blocks, "covar-node", covariates = list(dummy_party), control = list(useCov = FALSE, iterates = 2)) sbm_covar2$bestModel$fittedSampling$parameters # sampling parameters ## [1] -0.3629055 2.9469030 iii) Take the covariate into account in the SBM only: sbm_covar3 <-estimateMissSBM(blog_subgraph_obs, blocks, "node", covariates = list(dummy_party), control = list(useCov = TRUE, iterates = 2)) sbm_covar3$bestModel$fittedSampling$parameters # sampling parameter sbm_covar3$bestModel$fittedSBM$covarParam # regression parameter ## [1] 0.71875 ## [1] 4.945801 iv) Ignore the covariate in both the sampling and the SBM: sbm_covar4 <-estimateMissSBM(blog_subgraph_obs, blocks, "node", control = list(useCov = FALSE, iterates = 2)) sbm_covar4$bestModel$fittedSampling$parameters # sampling parameter ## [1] 0.71875

Figure 10 :

 10 Figure 10: ICL criterion for different numbers of blocks under the four models which make different use of the covariate political party.

Table 2 :

 2 Selection of fields in object missSBM_fit with descriptions and types.

	R6 object	Field	Description	Correspondence
		probMemberships estimated probability of block belonging	{τ i}i∈N
		connectParam	connectivity matrix	π
	fittedSBM	expectation covarParam	imputed adjacency matrix regression parameter	Y o ∪ {νij} (i,j)∈D m β
		memberships	vector of blocks memberships	(which.max(τ i))i∈N
		blockProp	block proportions	α
	fittedSampling	parameters	sampling parameter	ψ

Table 3 :

 3 Selection of fields in fittedSBM, fittedSampling and mathematical counterparts.

Extraction of the SBM with Block-Sampling Design.

  The model that we finally retain is thus a block-sampling with 9 blocks.

	=====================================================================
	Dimension = ( 194 ) -( 9 ) blocks and no covariate(s).
	=====================================================================
	* Useful fields
	$nbNodes, $modelName, $dimLabels, $nbBlocks, $nbCovariates, $nbDyads
	$blockProp, $connectParam, $covarParam, $covarList, $covarEffect
	$expectation, $indMemberships, $memberships
	* R6 and S3 methods
	$rNetwork, $rMemberships, $rEdges, plot, print, coef
	myModel$fittedSampling
	block-node-model for network sampling
	==================================================================
	Structure for handling network sampling in missSBM.
	==================================================================
	* Useful fields
	$type, $parameters, $df
	$penalty, $vExpec
	Representation and Validation.
	myModel <-sbm_block$bestModel
	myModel is an object with class missSBM_fit with two special elements used for storing
	the results of the estimation of both the SBM (field fittedSBM) and the sampling design
	(fittedSampling). These two elements are special objects themselves with dedicated fields
	and methods which are recalled to the user thanks to the print/show methods:
	myModel$fittedSBM
	Simple Stochastic Block Model --bernoulli variant

Table 4 :

 4 Clustering comparison with adjusted Rand indices (ARIs) between models taking the covariate into account (models i), ii), iii), iv) and model with fully observed network).

		i	ii	iii	iv
	Full 0.64 0.48 0.64 0.44
	i	NA 0.42 1.00 0.39
	ii	NA NA 0.42 0.95
	iii	NA NA NA 0.39

In addition to this brief typology, the interested reader may consult the CRAN task view on the related topic of graphical modeling(Hojsgaard

2019).

Although self-edges (Yii) could be defined in the SBM, they are not considered in missSBM since they are scarce in real data.

Use multisession if you work on Windows or from Rstudio

Earlier versions of this data set were available in packages mixer and sand.

Acknowledgments

The authors thank all members of MIRES group for fruitful discussions on network sampling designs. This work is supported by public grants overseen by the French National research Agency (ANR) 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, ANR-18-CE02-0010.

Affiliation:

Pierre Barbillon, Julien Chiquet & Timothée Tabouy MIA Paris, Université Paris-Saclay, AgroParisTech, INRAE E-mails: pierre.barbillon@agroparistech.fr, julien.chiquet@inrae.fr, timothee.tabouy@gmail.com


Manon Blaise
email: manon.blaise@etu.unistra.fr

Franz Chouly
email: franz.chouly@u-bourgogne.fr

Pierre-Yves Rohan
email: pierre-yves.rohan@ensam.eu

A few remarks about the computer implementation and the verification of some hyperelastic constitutive laws and an illustration with the mechanical response of an artery

In this short technical note, we are interested in the constitutive equations used to model macroscopically the mechanical function of soft tissues under mechanical loading. Soft tissues have the ability to undergo large elastic reversible deformations under quasi-static loading and are usually modelled using hyperelastic constitutive laws. Several constitutive equations have been defined in the literature by means of strain components or strain invariants isotropic hyperelastic models. This short technical note recalls how to derive the weak form for hyperelasticity equations, and how this weak form can be linearized. It also presents a semi-analytical solution for a simplified geometry representing the mechanical response of an artery modelled under pressure loading.

Introduction

As recalled in the book "Biomechanics of Living Organs: Hyperelastic Constitutive Laws for Finite Element Modeling" by Yohan Payan and Jacques Ohayon [START_REF] Payan | Biomechanics of Living Organs -Hyperelastic Constitutive Laws for Finite Element Modeling[END_REF], the mechanics of human soft tissues has been an emerging research field since the publication, in 1981, of the book "Biomechanics: Mechanical Properties of Living Tissues" by Yuan-Cheng Fung [START_REF] Fung | Biomechanics: mechanical properties of living tissues[END_REF]. Since that date, many groups in the world have proposed biomechanical models of soft organs to investigate both the underlying mechanisms that drive normal physiology of biological soft tissues and the mechanical factors that contribute to the onset and development of diseases such as tumour growth [START_REF] Yankeelov | Clinically relevant modeling of tumor growth and treatment response[END_REF], atherosclerosis or aneurysms [START_REF] Romo | In vitro analysis of localized aneurysm rupture[END_REF], or multilevel lumbar disc degenerative diseases [START_REF] Schmidt | Effect of multilevel lumbar disc arthroplasty on spine kinematics and facet joint loads in flexion and extension: a finite element analysis[END_REF], to name a few. Such organs indeed deform under physiological conditions (such as muscle activations or interactions with other tissues) or because of the mechanical interaction with the surgical gesture.

Assuming we know the governing ordinary and partial differential equations, finite element models can predict the behavior of the system from given initial and boundary conditions measured at a few selected points. This approach is incredibly powerful, but, to close the system of equations, we need constitutive equations that characterize the behavior of the system, which we need to calibrate with experimental data, for instance.

In the last decades, significant theoretical results, many confirmed by experiments, have projected considerable light on the mechanical behaviour of elastic materials subjected to large deformations. It has consistently been shown that hyperelasticity can be used to model macroscopically the mechanical function of soft tissues under mechanical loading.

In this document we focus on general hyperelasticity formulations, and detail the derivation of the weak form, and then of the tangent weak problem, for the deformation of an hyperelastic body. These are the two fundamental ingredients for an implementation within a finite element library. Nowadays modern finite element libraries such as FEniCS [START_REF] Logg | Automated solution of differential equations by the finite element method: The FEniCS book[END_REF] or GetFEM++ [START_REF] Renard | GetFEM: automated FE modeling of multiphysics problems based on a generic weak form language[END_REF] make use of high level assembling languages and even of symbolic automatic differentiation, and are capable of solving an hyperelastic problem by specifying only the complete expression of the hyperelastic potential. Nevertheless, even in this case, it can be interesting to know how to carry out by hand this derivation, for instance to gain more control on the solving procedure. Furthermore, the linearized hyperelasticity problem is involved in adjoint-based solvers, that are useful in many applications, such as parameter calibration, shape optimization or dual-based error estimation, see, e.g., [START_REF] Allaire | A review of adjoint methods for sensitivity analysis, uncertainty quantification and optimization in numerical codes[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF]. This first topic that concerns the weak form and its linearization is of course covered in many documents about hyperelasticity and the presentation here is mostly inspired from [START_REF] Bonet | Nonlinear continuum mechanics for finite element analysis[END_REF][START_REF] Fortin | Les éléments finis: de la théorie à la pratique[END_REF].

The second part of this document details a semi-analytical solution for the hyperelastic deformation of a ring, that may represent an idealized artery. This can be helpful to verify numerical solutions with finite elements, in the same spirit as manufactured solutions [START_REF] Chamberland | Comparison of the performance of some finite element discretizations for large deformation elasticity problems[END_REF]. This is a first step towards semi-analytical solutions with active constitutive laws such as presented in [START_REF] Payan | Biomechanics of Living Organs -Hyperelastic Constitutive Laws for Finite Element Modeling[END_REF].

Hyperelasticity

We first recall that, when a material is purely elastic, all the energy that is supplied to it during loading can be regained after unloading. So, during the whole cycle of deformation, no energy is dissipated physically. It can be shown that the stress in such a material can be derived from a scalar functional representing the specific elastic energy density (elastic energy per unit undeformed volume).

Notations

Let Ω ⊂ R d be an open bounded set that denotes the reference configuration of one deformable solid in a space of dimension d = 2 or 3. A deformed configuration Ω t can be defined through a transformation ϕ which maps any point X of the reference configuration to a point x of the deformed one:

ϕ : Ω -→ R d X -→ x = ϕ(X).
We define the displacement u relatively to the reference configuration as:

u(X) := ϕ(X) -X.
The gradient of a quantity in the deformed (resp. reference) configuration will be noted ∇ (resp. ∇ X ). To describe the deformation, we introduce as usual the identity matrix of size d×d, denoted by I and the deformation gradient F = I + ∇ X u. We introduce as well

(u) := sym (∇u) = 1 2 ∇u + ∇u T = 1 2 (∇ X u)F -1 + F -T (∇ X u) T = sym ((∇ X u)F -1 )
where sym (•) denotes the symmetric part of a second order tensor. The jacobian of ϕ is denoted by J := det F. We introduce also the Cauchy-Green tensor C = F T F, and the Green-Lagrange tensor E = 1 2 (C -I). We will note σ the Cauchy stress tensor, Π = JσF -T the first Piola-Kirchhoff stress tensor and Σ = JF -1 σF -T the second Piola-Kirchhoff stress tensor.

Hyperelastic constitutive law

We consider a general hyperelastic constitutive law, derived from a potential W that depends on the deformation through E (or C) (see, e.g., [START_REF] Bonet | Nonlinear continuum mechanics for finite element analysis[END_REF][START_REF] Fortin | Les éléments finis: de la théorie à la pratique[END_REF]), so that the second Piola-Kirchhoff stress is

Σ = ∂W ∂E = 2 ∂W ∂C ,
with corresponding fourth-order elasticity tensor

C := ∂Σ ∂E = ∂ 2 W ∂E∂E .
We will need as well the isotropic tensor

I := 1 2 (e i ⊗ e j ⊗ e i ⊗ e j + e i ⊗ e j ⊗ e j ⊗ e i ),
where ⊗ denotes the tensor product of two vectors, (e i ) i=1,...,d is the canonical basis of R d and where Einstein's summation convention is used. The tensor I has the property I : A = A for any symmetric second-order tensor A (: denotes the double-dot product between two tensors).

We will denote the global potential energy by J (•). For example, if considering simple equilibrium under a gravity force, the potential energy is

J (u) = Ω W(E(u)) dX - Ω ρ g • u dX, ( 1
)
where ρ is the density in the reference configuration and g is the gravity acceleration vector. Of course, additional terms such as boundary loads, can be considered as well. Dirichlet conditions can also be prescribed, but, to simplify the formulation, the treatment of Dirichlet conditions will be omitted in the following. Moreover, we consider the following notations to simplify the mathematical presentation. The directional derivative of a quantity A with respect to the displacement u in direction δu will be denoted by A (u; δu). This directional derivative is defined as

A (u; δu) := lim ε→0 A(u + ε δu) -A(u) ε = d dε A(u + ε δu)| ε=0
when this limit exists.

Having some fun with differential calculus

We need first the following results on the directional derivatives of deformation tensors, that are obtained after simple computations (see, e.g. [6, Chapter 14]):

Proposition 2.1. The directional derivatives of F and E are:

F (u; δu) = ∇ X (δu), E (u; δu) = sym (F T (u)∇ X (δu)) = F T (u) (δu) F(u).
Proof: The expression of F (u; δu) is a direct consequence of the definition of F. Then let us compute

E (u; δu) = 1 2 C (u; δu) = 1 2 (F T ) (u; δu)F(u) + F T (u)F (u; δu) = 1 2 (∇ T X ) (δu)F(u) + F T (u)∇ X (δu)
where we used first the definition of E and C, then the product rule and finally the expression of F . As a result we recover the first expression. To get the second expression we can do as follows

E (u; δu) = sym (F T (u)∇ X (δu)) = sym (F T (u)∇(δu)F(u)) = F T (u) sym (∇(δu))F(u)
where we used ∇ X (δu) = ∇(δu)F(u) and the properties of the symmetric part.

The computation of the directional derivatives of stress tensors Π and Σ is more involved and we recall their expression below (see as well [START_REF] Fortin | Les éléments finis: de la théorie à la pratique[END_REF]Chapter 14]): Proposition 2.2. The directional derivatives of Σ, resp. Π, are:

Σ (u; δu) = C(u) : F T (u)∇ X (δu) = C(u) : F T (u) (δu) F(u).
(2)

Π (u; δu) = ∇ X (δu) Σ(u) + F(u)(C(u) : F T (u)∇ X (δu)). (3)
Proof: First, for an hyperelastic law, there holds in fact Σ(u) = Σ(E(u)) and we apply the chain rule:

Σ (u; δu) = Σ (E(u); E (u; δu)) = ∂Σ ∂E (E(u)) : E (u; δu).
Since C = ∂Σ/∂E and using Proposition 2.1 we get:

Σ (u; δu) = C(u) : sym (F T (u)∇ X (δu)).
Then ( 2) is obtained with the above formula and the symmetry properties of C. Using the relationship Π = FΣ and applying the product rule yield:

Π (u; δu) = F (u; δu)Σ(u) + F(u)Σ (u; δu).
We use (2) and once again Proposition 2.1 to obtain (3).

Weak formulation

Now we are ready to find the expression of the weak form for an hyperelastic material. This weak form is the first order optimality condition associated to the minimization of (1) and reads, for any virtual displacement δu:

0 = J (u; δu) = Ω (W • E) (u; δu) dX - Ω ρ g • δu dX.
We compute, using the chain rule and Proposition 2.1:

(W • E) (u; δu) = W (E(u); E (u; δu)) = ∂W ∂E (E(u)) : E (u; δu) = Σ(u) : sym (F T (u)∇ X (δu)).
Using the symmetry property of Σ and the relationship Π = FΣ, we obtain:

(W • E) (u; δu) = Σ(u) : F T (u)∇ X (δu) = F(u)Σ(u) : ∇ X (δu) = Π(u) : ∇ X (δu).
Finally we define (substituting v to δu)

a(u; v) = Ω Π(u) : ∇ X v dX, L(v) = Ω ρ g • v dX,
and the weak form for an hyperelastic material reads, for all v:

a(u; v) = L(v). (4)
This is the kind of expression that can be implemented in GetFEM++ or FEniCS, as an alternative to the expression of W.

Linearization for hyperelasticity

We want to linearize the weak form (4). This linearization can be helpful for different purposes. Particularly it is involved in iterative procedures for numerical solving, such as the Newton method, that needs the exact jacobian of the weak form. It can also be useful for the theoretical study of linearized systems around some states or prestress calculations. It has also an interest for adjoint-based solvers when doing dual based a posteriori error estimation, sensitivity analysis, parameter calibration or optimization. As already mentioned, some computer software has now the capability of carrying out this linearization automatically, but it can be necessary sometimes, or just more amusing, to do this by hand. We can proceed as follows.

Details of the linearization procedure

Since L(•) is linear, there only remains to compute the directional derivative of a(•; •). We fix u and v and consider a direction of differentiation δu, and first note that:

a (u; δu, v) = Ω Π (u; δu) : ∇ X v dX.
There just remains to apply (3):

a (u; δu, v) = Ω ∇ X (δu) Σ(u) + F(u)(C(u) : F T (u)∇ X (δu)) : ∇ X v dX.
Splitting the sum, there holds:

a (u; δu, v) = Ω ∇ X (δu) Σ(u) : ∇ X v dX + Ω F(u)(C(u) : F T (u)∇ X (δu)) : ∇ X v dX.
This can be rewritten as (using symmetries of C):

a (u; δu, v) = Ω Σ(u) : ∇ X (δu) T ∇ X v dX + Ω F T (u)∇ X (δu) : C(u) : F(u) T ∇ X v dX. ( 5
)
Remark that for a Newton-Raphson iteration, at each step, and for a given u, we will find an increment δu solution to the tangent system:

a (u; δu, v) = L(v) -a(u; v).

Re-interpretation "à la Schotté" of the tangent stiffness weak form

We can provide a beatiful decomposition of the tangent stiffness system, following [START_REF] Schotté | Various modelling levels to represent internal liquid behaviour in the vibration analysis of complex structures[END_REF]. First we use the relationship F = I + ∇ X , and reformulate the term:

Ω F T (u)∇ X (δu) : C(u) : F(u) T ∇ X v dX = Ω (I + ∇ T X (u))∇ X (δu) : C(u) : (I + ∇ T X (u))∇ X v dX = Ω ∇ X (δu) : C(u) : ∇ X v dX + Ω ∇ X (δu) : C(u) : ∇ T X (u)∇ X v dX + Ω ∇ T X (u)∇ X (δu) : C(u) : ∇ X v dX + Ω ∇ T X (u)∇ X (δu) : C(u) : ∇ T X (u)∇ X v dX.
With the above calculation, we can rewrite the tangent weak form as:

a (u; δu, v) = a G (u; δu, v) + a L (u; δu, v) + a D (u; δu, v), where a G (u; δu, v) := Ω Σ(u) : ∇ X (δu) T ∇ X v dX is the geometric stiffness contribution, a L (u; δu, v) := Ω ∇ X (δu) : C(u) : ∇ X v dX
is the linear stiffness contribution and

a D (u; δu, v) := Ω ∇ X (δu) : C(u) : ∇ T X (u)∇ X v dX + Ω ∇ T X (u)∇ X (δu) : C(u) : ∇ X v dX + Ω ∇ T X (u)∇ X (δu) : C(u) : ∇ T X (u)∇ X v dX.
is the predeformation stiffness contribution. For the terminology and more explanations, see [START_REF] Schotté | Various modelling levels to represent internal liquid behaviour in the vibration analysis of complex structures[END_REF].

Saint-Venant-Kirchhoff material

As an example, and as in [START_REF] Mlika | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF], suppose that the constitutive law is those of a Saint-Venant-Kirchhoff material, i.e., that:

W(E) = λ 2 (tr (E)) 2 + µ tr (E 2 ),
where λ and µ are material parameters, see, e.g., [START_REF] Bonet | Nonlinear continuum mechanics for finite element analysis[END_REF]Chapter 5]. The associated second Piola-Kirchhoff stress tensor and elasticity tensor are:

Σ = λtr (E)I + 2µE, C = λ I ⊗ I + 2µ I.
Let us detail the expression

Π (u; δu) = ∇ X (δu) Σ(u) + F(u)(C(u) : F T (u)∇ X (δu)) = ∇ X (δu) (λtr (E(u))I + 2µE(u)) + F(u)((λ I ⊗ I + 2µ I) : F T (u)∇ X (δu)).
We compute separately

I ⊗ I : F T (u)∇ X (δu) = (I : F T (u)∇ X (δu))I = tr (F T (u)∇ X (δu))I,
and

I : F T (u)∇ X (δu) = sym (F T (u)∇ X (δu)).
This yields

Π (u; δu) = ∇ X (δu) (λtr (E(u))I + 2µE(u)) + F(u)(λ tr (F T (u)∇ X (δu))I + 2µ sym (F T (u)∇ X (δu))) = λ tr (E(u))∇ X (δu) + tr (F T (u)∇ X (δu))F(u) + 2µ ∇ X (δu)E(u) + F(u) sym (F T (u)∇ X (δu)) .

Inflation of an artery modelled as a hyperelastic cylinder : semi-analytical solution

The prediction of the stress distribution within arteries has been a major focus of cardiovascular biomechanics. The complexity of the physical system presents a formidable modelling challenge. The walls of blood vessels are mixtures of solid and fluid constituents, with a very sizable fraction of fluid content. The solid constituent of soft tissues is comprised of extra-cellular matrix containing various types of collagen and elastin, and cells which also contain fluids. Blood vessels are anisotropic and the anisotropy changes from point to point, and complicating matters further, blood vessels are inhomogeneous. To make some initial progress, idealisations of various degrees of realism have necessarily been made in the litterature. The larger arteries are commonly idealised as being cylindrical in shape in cardiovascular biomechanics and assumed to be both incompressible and non-linearly elastic. The complex loadings experienced by arteries in situ are typically modelled using the approximation of uniform internal inflation. It is commonly assumed that the artery is cylindrically isotropic.

All of these idealisations are adopted in this section. The artery is modelled using a strain-energy function per unit undeformed volume. We derive here a semi-analytical solution for a simplified geometry (a ring).

Geometry and boundary conditions

We consider a perfect ring that represents the transverse section of an idealized arterial wall, see Figure 1. This configuration is established as the reference one Ω. The notations Γ i and Γ e define internal and external boundaries, respectively. A predefined pressure p ext is imposed on the external boundary Γ e (where x = r e e r ):

σ n = p ext n. (6)
On Γ i (where x = r i e r ) we can impose either a Neumann (pressure) or Dirichlet (displacement) boundary condition, as:

σ n = p int n, (7)
u = u int n, (8)
where p int and u int are fixed internal pressure and displacement, and n normal vector pointing outward to the boundary, which is either e r for Γ e or -e r for Γ i .

Incompressible neo-hookean material

In biomechanical studies, arterial tissue is commonly assumed to be incompressible. A challenge of this study is also to determine relevance of this incompressibility condition. Indeed, although ex-vivo experiences confirm it, in-vivo tests -more invasive -could be in favour of a quasi-incompressible model (see, e.g. [START_REF] Bouten | Identification des propriétés mécaniques des tissus constitutifs du mollet pour l´étude mécanique de la contention[END_REF]). Considering an incompressible neo-hookean material, the associated constitutive law yields:

W(C) = µ 2 (tr (C) -d).
Second Piola-Kirchhoff tensor is computed as follows:

Σ = 2 ∂W ∂C -pC -1 ,
where p is the hydrostatic pressure defined in each point of Ω. Plus, we introduce the following condition for incompressibility (conservation of volume):

J = det(F) = 1. (9)

Equilibrium equation

Supposing the absence of body forces on the whole domain Ω, the fundamental principle of statics is:

-div (σ) = 0. (10)

First calculations in cylindrical coordinates

The gradient of the displacement u is given in cylindrical coordinates by:

∇u =         ∂u r ∂r 1 r ∂u r ∂θ -u θ ∂u r ∂z ∂u θ ∂r 1 r ∂u θ ∂θ + u r ∂u θ ∂z ∂u z ∂r 1 r ∂u z ∂θ ∂u z ∂z         .
In the present case, the displacement u is unidirectional and depends only on r. Hence as sketched in Fig. 1, we assume that:

∀M (r, θ, z), u(M ) = u r (r) e r .
The gradient of the displacement can be directly simplified as follows:

∇u =     du r dr 0 0 0 u r r 0 0 0 0     .
The tensors of deformation F = I + ∇u, Cauchy C and Green-Lagrange E are derived in Eqs. ( 11) to ( 13) respectively.

F =   λ r 0 0 0 λ θ 0 0 0 1   (11)
C =   λ 2 r 0 0 0 λ 2 θ 0 0 0 1   (12)
E = 1 2        2 du r dr + du r dr 2 0 0 0 2 u r r + u r r 2 0 0 0 0        =        du r dr + 1 2 du r dr 2 0 0 0 u r r + 1 2 u r r 2 0 0 0 0        . ( 13
)
where λ r = 1 + du r dr and λ θ = 1 + u r r .

Expression for the displacement

First, the displacement u r is obtained by solving the differential equation of incompressibility:

0 = det(F) -1 = λ r λ θ -1 = du r dr u r r + du r dr + u r r .
We assume that a solution is

u r (r) = c + r 2 -r,
where c is a constant which will be noted c = r 2 o . Let us proceed with intermediate calculations:

du r dr = 2r 2 √ c + r 2 -1 det(F) -1 = r √ c + r 2 -1 √ c + r 2 r -1 + r √ c + r 2 -1 + √ c + r 2 r -1 = - r √ c + r 2 + 1 - √ c + r 2 r + 1 + r √ c + r 2 + √ c + r 2 r -2 = 0
So we checked that the expression for u r is correct. Then the constant c is determined using the boundary conditions.

• In the case where u r (r i ) is known, we consider Eq. ( 8), so that:

u r (r i ) = -u int = 2 o + r 2 i -r i .
Then we find

r 2 o = u int (u int -2r i ).
Finally, response of the displacement is defined for r i r r e and given by:

u r (r) = u int (u int -2r i ) + r 2 -r. (14)
• In the case where u r (r i ) is unknown, r 2 o is determined by pressure conditions Eqs. ( 6) and ( 7). Hence, u r is noted: u r (r) = r 2 + r 2 o -r.

Stress and equilibrium

We compute the Cauchy stress tensor σ under the incompressibility hypothesis Eq. ( 9):

σ = 1 J FΣF T = FµIF T -pFC -1 F T = µC -pFF -1 F -T F T = µC -pI, (15)
i.e.:

σ =   µλ 2 r -p 0 0 0 µλ 2 θ -p 0 0 0 µ -p   . (16)
Then remember that the divergence of a second order tensor A in cylindrical coordinates is given by:

-→ div(A) =           ∂A rr ∂r + 1 r ∂A rθ ∂θ + ∂A rz ∂z + A rr -A θθ r ∂A θr ∂r + 1 r ∂A θθ ∂θ + ∂A θz ∂z + A rθ + A θr r ∂A zr ∂r + 1 r ∂A zθ ∂θ + ∂A zz ∂z + A zr r           .
For the Cauchy stress tensor, invariant in θ and z:

-→ div(σ) =        µ dλ 2 r dr - dp dr + 1 r µ λ 2 r -λ 2 θ 0 0        .
Applying the fundamental principle of statics Eq. ( 10) and incompressible equation Eq. ( 9), assuming that

λ r = 1 λ θ : dp dr = µ dλ 2 r dr + 1 r µ λ 2 r -λ 2 θ = µ d dr r 2 r 2 + r 2 o + 1 r r 2 r 2 + r 2 o -1 + r 2 o r 2 = µ d dr r 2 r 2 + r 2 o + 2r 2(r 2 + r 2 o ) - 1 r - r 2 o r 3 .
By integration, we obtain then:

p r (r) = µ r 2 r 2 + r 2 o + 1 2 ln(r 2 + r 2 o ) -ln(r) + r 2 o 2r 2 + p o = µ 2r 4 + r 2 o (r 2 + r 2 o ) 2r 2 (r 2 + r 2 o ) + ln r 2 + r 2 o r + p o , (17)
with p o the constant of integration. From Eqs. ( 6) and ( 7) on boundary conditions, with an external pressure p ext = 0:

-(µλ 2 ri -p(r i )) = p int , µλ 2
re -p(r e ) = 0. Thus

p int = -µ r 2 i r 2 i + r 2 o - 2r 4 i + r 2 o (r 2 i + r 2 o ) 2r 2 i (r 2 i + r 2 o ) -ln r 2 i + r 2 o r i + p o (18)
with

p o = µ
We sum Eq. (18) + Eq. ( 19) to get the expression of p int in function of r i , r e and r o :

p int = -µ r 2 i r 2 i + r 2 o - r 2 e r 2 e + r 2 o + µ 2 2r 4 i + r 2 o (r 2 i + r 2 o ) r 2 i (r 2 i + r 2 o ) - 2r 4 e + r 2 o (r 2 e + r 2 o ) r 2 e (r 2 e + r 2 o ) + µ ln r e r i r 2 i + r 2 o r 2 e + r 2 o .
Finally, from:

p r (r) = µ 2r 4 + r 2 o (r 2 + r 2 o ) 2r 2 (r 2 + r 2 o ) + ln r 2 + r 2 o r + p o ,
the final expression of p r is given by:

p r (r) = µ 2r 4 + r 2 o (r 2 + r 2 o ) 2r 2 (r 2 + r 2 o ) + r 2 e r 2 e + r 2 o - 2r 4 e + r 2 o (r 2 e + r 2 o ) 2r 2 e (r 2 e + r 2 o ) + ln r 2 + r 2 o r -ln r 2 e + r 2 o r e .
Hence, we obtain the expressions of the principal radial stress tensor (cf Eq. ( 15)), and therefore derived expressions such as the Von Mises stress, etc. Wonderful, isn't it ?

In order to find a numerical value of r o :

We introduce the following function Φ, such as Φ(r o ) = 0:

Φ(x) = -p int -µ r 2 i r 2 i + x 2 - r 2 e r 2 e + x 2 + µ 2 2r 4 i + x 2 (r 2 i + x 2 ) r 2 i (r 2 i + x 2 ) - 2r 4 e + x 2 (r 2 e + x 2 ) r 2 e (r 2 e + x 2 ) +µ ln r e r i r 2 i + x 2 r 2 e + x 2
. By dichotomy, a numerical value of r o can be obtained for r e = 2.25 mm and r i = 1.93 mm: r o = 0.01707...

Numerical application

Let us take a very simple example of an expansion movement. The value for the parameters are: r i = 4.0 mm, r e = 7.0 mm, and u int = -2.0 mm. See configuration Fig. 2: For the next figure 3, we use following values: r i = 1.93 mm, r e = 2.25 mm, and µ = 27.9 kP a, see [START_REF] Payan | Biomechanics of Living Organs -Hyperelastic Constitutive Laws for Finite Element Modeling[END_REF].

Saint-Venant-Kirchhoff material

As an other example, suppose that the constitutive law is those of a Saint-Venant-Kirchhoff material. According to the constitutive law adapted to this case Σ = λtr (E)I + 2µE, Second Piola-Kirchhoff stress is written:

Σ =     λtr (E) + 2µ du r dr + µ( du r dr ) 2 0 0 0 λtr (E) + 2µ u r r + µ( u r r ) 2 0 0 0 λtr (E)    
and allows to find the first Piola-Kirchhoff stress tensor Π = FΣ:

FΣ =     (1 + du r dr )[λtr (E) + 2µ du r dr + µ( du r dr ) 2 ] 0 0 0 (1 + u r r )λtr (E) + 2µ u r r + µ( u r r ) 2 0 0 0 λtr (E)     . Thus -→ div(FΣ) =    ∂FΣ 11 ∂r + 1 r (FΣ 11 -FΣ 22 ) 0 0    and tr (E) = du r dr + 1 2 ( du r dr ) 2 + u r r + 1 2 ( u r r ) 2 .
Term by term, calculations show that:

FΣ
+ (λ + 3µ) ( dur dr ) 2 -( ur r ) 2 + ( λ 2 + µ) ( dur dr ) 3 -( ur r ) 3 .
Furthermore

FΣ 11 = (1 + du r dr ) λ du r dr + 1 2 ( du r dr ) 2 + u r r + 1 2 ( u r r ) 2 + µ(2 du r dr + ( du r dr ) 2 ) = λ u r r + 1 2 ( u r r ) 2 + du r dr + 3 2 ( du r dr ) 2 + 1 2 ( du r dr ) 3 + u r r du r dr + 1 2 ( u r r ) 2 du r dr + µ 2 du r dr + 3( du r dr ) 2 + ( du r dr ) 3 .
Some intermediates calculations could be helpful:

d dr ( u r r ) = 1 r ( du r dr - u r r ), d dr ( u r r ) 2 = 2 u r r 2 ( du r dr - u r r ), d dr du r dr u r r = u r r d 2 u r dr 2 + 1 r ( du r dr ) 2 -( du r dr u r r ) , d dr ( u r r ) 2 du r dr = ( u r r ) 2 d 2 u r dr 2 + 2 u r r 2 ( du r dr - u r r )
and

dFΣ 11 dr = λ r 1 + u r r du r dr - u r r + (2µ + λ) d 2 u r dr 2 + (6µ + 3λ) du r dr d 2 u r dr 2 + 3( λ 2 + µ)( du r dr ) 2 d 2 u r dr 2 + λ u r r d 2 u r dr 2 + 1 r ( du r dr ) 2 - u r r du r dr + 1 2 ( u r r ) 2 d 2 u r dr 2 + 1 r u r r du r dr ( du r dr - u r r ) .
We finally obtain:

-→ div(FΣ)   1 0 0   = -( 2µ + λ r ) u r r + ( 2µ + λ r ) du r dr + (2µ + λ)( du r dr ) 2 -( 2λ + 3µ r )( u r r ) 2 - 1 r ( λ 2 + µ)( u r r ) 3 + ( 2λ + 3µ r )( du r dr ) 2 + 1 r ( λ 2 + µ)( du r dr ) 3 - 1 r (1 - λ 2 ) ( u r r ) 2 du r dr + 1 r (1 - λ 2 ) u r r ( du r dr ) 2 + λ u r r d 2 u r dr 2 + λ 2 ( u r r ) 2 d 2 u r dr 2 + (6µ + 3λ) du r dr d 2 u r dr 2 + 1 2 (6µ + 3λ)( du r dr ) 2 d 2 u r dr 2 (20)
Then, using the above expression and the equilibrium equation Eq. ( 10), an approximated numerical solution could be obtained.

Conclusion

In this short technical note, after a brief review of the general hyperelasticity framework and the detail of the derivation of the weak form, the theory was applied to a the case of the mechanical response of an artery modelled under pressure loading. This can be helpful to verify numerical solutions with finite elements, in the same spirit as manufactured solutions [START_REF] Chamberland | Comparison of the performance of some finite element discretizations for large deformation elasticity problems[END_REF]. This is a first step towards semi-analytical solutions with active constitutive laws such as presented in [START_REF] Payan | Biomechanics of Living Organs -Hyperelastic Constitutive Laws for Finite Element Modeling[END_REF].

On Γ i , where x = r i e r , we apply either:

Π n( x) = -p int n (A.2) or u r (r i ) = u int n, (A.3)
where p int and u int are fixed internal pressure and displacement, and n normal vector pointing outward to the boundary, which is either e r for Γ e or -e r for Γ i .

Appendix A.2. Equilibrium We start from

-→ div(Π) = -→ div(FΣ) =    ∂FΣ 11 ∂r + 1 r (FΣ 11 -FΣ 22 ) 0 0    (A.4)
and compute second Piola-Kirchhoff tensor as follows:

Σ = 2 ∂W ∂C -pC -1 = 2 ∂ ∂C µ 2 (tr(C) -d -pC -1 = µ ∂(C kk -d) ∂C -pC -1 = µI -pC -1 C -1 =         1 (1 + du r dr ) 2 0 0 0 1 (1 + u r r ) 2 0 0 0 1         Σ =        µ - p (1 + du r dr ) 2 0 0 0 µ - p (1 + u r r ) 2 0 0 0 µ -p       
in order to obtain:

FΣ =         µ(1 + dur dr ) - p (1 + dur dr ) 0 0 0 µ(1 + ur r ) - p (1 + ur r ) 0 0 0 µ -p         dFΣ11 dr + 1 r (FΣ11 -FΣ22) = µ d 2 ur dr 2 + p d 2 ur dr 2 1 1 + dur dr 2 - dp dr 1 1 + dur dr + 1 r     µ dur dr - ur r -p     1 (1 + dur dr ) - 1 (1 + ur r )         (A.5)
It results from Equilibrium equation Eq. ( 10) combined with Eqs. (A.4) and (A.5):

     µ + p 1 + du r dr 2      d 2 u r dr 2 - dp dr 1 1 + du r dr + µ r du r dr - u r r - p r     1 1 + du r dr - 1 1 + u r r     = 0 (A.6)
This equation will give the differential equation of u r , as a function of hydrostatic pressure p. Let us finally take the result of the fundamental principle of statics given by Eq. ( 10):

p      d 2 u r dr 2 1 + du r dr 2 + 1 r    1 1 + u r r - 1 1 + du r dr         - dp dr 1 1 + du r dr = µ r u r r - du r dr -µ d 2 u r dr 2 .

Appendix B. A few lines of Python

We provide here a Python script for the numerical application. + mu * l o g ( s q r t ( Ri+Ro ) / r i ) + po # v e r i f i c a t i o n : a l p h a = p i n t return ( mu -p r ( r ) ) def sigma VM ( r ) :

return ( 1 / ( s q r t ( 2 ) ) * s q r t ( ( s i g m a r r ( r )-s i g m a t t ( r ) ) * * 2 + ( s i g m a t t ( r )-s i g m a z z ( r ) ) * * 2 + ( s i g m a z z ( r )-s i g m a r r ( r ) ) * * 2 ) )

# C o n t r a i n t e r a d i a l e Si g m a r r ( P r i n c i p a l r a d i a l s t r e s s ) p l t . f i g u r e (

Figure 1 :

 1 Figure 1: Reference configuration Ω0

Figure 2 :

 2 Figure 2: Visualisation of displacement from the reference to the "current" configuration under pressure.

Figure 3 :

 3 Figure 3: Hydrostatic pressure at 20, 40 and 60 mmHg

Listing 1 :

 1 ro po import m a t p l o t l i b . p y p l o t a s p l t import math import numpy a s np from p y l a b import * # I l s ' a g i t de d e t e r m i n e r d e s v a l e u r s numé r i q u e s approch é e s d e s c o n s t a n t e s ro e t po # Ceci p e r m e t t a n t e n s u i t e de c a l c u l e r en t o u t p o i n t du domaine : # l e dé placement , l a p r e s s i o n h y d r o s t a t i q u e , # e t d ' a u t r e s paramè t r e s t e l s que l e s c o n t r a i n t e s p r i n c i p a l e s , # c o n t r a i n t e de Von Mises , dé f o r m a t i o n . . . # ---Paramè t r e s ( u n i t é SI ) r i = 1 . 9 3 e-3 # en m. Rayon i n t e r n e du c y l i n d r e p a r f a i t r e = 2 . 2 5 e-3 # rayon e x t e r n e Ri = r i * r i Re = r e * r e r 1 = 1 . 9 9 e-3 # p o i n t s de c a l c u l du modele numerique r e a l i t e , c e t t e v a l e u r de 2 7 . 9 kPa # c o r r e s p o n d a une donnee pour l ' intima , e t non l a media # C o n d i t i o n s de b o r d s en p r e s s i o n . # on e t u d i e une d i f f e r e n c e de p r e s s i o n : 100-80 = 20 mmHg p i n t = 2 . 6 6 6 e3 #1 3 . 3 3 e3 = 100 mmHg p e x t = 0 #1 0 . 6 6 4 e3 = 80 mmHg # ---Fonction pour d e t e r m i n e r ro ---( Phi ( ro ) = 0) def Phi ( x ) : #t q Phi ( ro )=0 terme mu = -mu * ( ( Ri / ( Ri+x * * 2 ) ) -( Re / ( Re+x * * 2 ) ) ) ; terme mu12 = (mu/ 2 ) * ( ( 2 * Ri * Ri + x * * 2 * ( Ri+x * * 2 ) ) / ( Ri * ( Ri+x * * 2 ) ) -( 2 * Re * Re + x * * 2 * ( Re+x * * 2 ) ) / ( Re * ( Re+x * * 2 ) ) ) ; t e r m e l n = mu * l o g ( s q r t ( ( Re * ( Ri+x * * 2 ) ) / ( Ri * ( Re+x * * 2 ) ) ) ) ; return ( -p i n t -p e x t + terme mu + terme mu12 + t e r m e l n ) ; # Fonction d i c h o t o m i e pour t r o u v e r l e z e r o d ' une f o n c t i o n f def d i c h o ( f , a , b , p r e c ) : while b-a>p r e c : c = ( a+b ) / 2 i f f ( a ) * f ( c ) <= 0 : b = c e l s e : a = c return a # Determinatio n d e s c o n s t a n t e s --ro-e t --po-r o = d i c h o ( Phi , 0 , 0 . 1 , 1 e -15) #a t t e n t i o n , ro e s t une c o n s t a n t e en m Ro = r o * * 2 po = mu * ( Re / ( Re+Ro ) ) -(mu/ 2 ) * ( 2 * Re * * 2 + Ro * ( Re+Ro ) ) / ( Re * ( Re+Ro ) ) -mu * l o g ( s q r t ( Re+Ro ) / r e ) + p e x t a l p h a= -mu * ( Ri / ( Ri+Ro ) ) + (mu/ 2 ) * ( 2 * Ri * * 2 + Ro * ( Ri+Ro ) ) / ( Ri * ( Ri+Ro ) )

4 )

 4 #p l t . s u b p l o t ( 1 2 1 ) p l t . g r i d( True ) x d e f = np . l i n s p a c e ( r i * 1 e3 , r e * 1 e3 , 1 0 0 0 ) y d e f = [ -s i g m a r r ( i * 1 e -3) * 1 e-3 f o r i in x d e f ] #y 3 d e f = [ s i g m a t t ( i * 1 e -3) f o r i i n x d e f ] #y 4 d e f = [ s i g m a z z ( i * 1 e -3) f o r i i n x d e f ] p l t . p l o t ( x d e f , y d e f , l a b e l = " A n a l y t i c a l S11 " ) #p l t . p l o t ( x d e f , y 3 d e f , l a b e l = " s i g m a t t ") #p l t . p l o t ( x d e f , y 4 d e f , l a b e l = " s i g m a z z ") p l t . x l i m ( r i * 1 e3 -0 . 0 5 , r e * 1 e3 + 0 . 0 5 ) x = np . a r r a y ( [ 1 . 9 3 , 1 . 9 9 , 2 . 0 6 , 2 . 1 2 , 2 . 1 9 , 2 . 2 5 ] ) # n u m e r i c a l model y = np . a r r a y ( [ 0 . 0 . p l o t ( x , y , ' r o ' , l a b e l =" N u m e r i c a l S11 " ) p l t i t l e ( " P r i n c i p a l R a d i a l S t r e s s ( kPa ) " ) p l t . l e g e n d ( l o c = ' b e s t ' ) #met l a l e g e n d e au m e i l l e u r e n d r o i t x l a b e l ( " r (mm) " ) y l a b e l ( " s i g m a 1 1 ( kPa ) " ) p l t . show ( ) # C o n t r a i n t e p r i n c i p a l e s e l o n t h e t a # # p l t . f i g u r e ( 4 ) # p l t . s u b p l o t ( 1 2 2 ) # p l t . g r i d ( True ) # x d e f = np . l i n s p a c e ( r i * 1 e3 , r e * 1 e3 , 1 0 0 0 ) # y 3 d e f = [ s i g m a t t ( i * 1 e -3) f o r i i n x d e f ] # p l t . p l o t ( x d e f , y 3 d e f , l a b e l = " s i g m a t t ") # #p l t . x l i m ( r i * 1 e3 -0.05 , r e * 1 e3 +0.05) # #p l t . y l i m ( s i g m a r r ( r i ) -1000 , s i g m a z z ( r e )+1000) # p l t . t i t l e (" s t r e s s ( r ) " ) # p l t . l e g e n d ( l o c = ' b e s t ' ) # x l a b e l (" r (mm) " ) # y l a b e l (" s i g m a t t (Pa ) " ) # p l t . show ( ) # A f f i c h a g e c o n t r a i n t e Von Mises en MPa p l t . f i g u r e ( 5 ) p l t . g r i d ( True ) x d e f = np . l i n s p a c e ( r i * 1 e3 , r e * 1 e3 , 1 0 0 0 ) y 5 d e f = [ ( sigma VM ( i * 1 e -3)) * 1 e-3 f o r i in x d e f ] p l t . p l o t ( x d e f , y 5 d e f , l a b e l = " 20 mmHg" ) p l t . x l i m ( 1 . 9 , 2 . 2 7 ) x = np . a r r a y ( [ 1 . 9 3 , 1 . 9 9 , 2 . 0 6 , 2 . 1 2 , 2 . 1 9 , 2 . 2 5 ] ) # n u m e r i c a l model y = np . a r r a y ( [ 0 . 0 . p l o t ( x , y , ' r o ' ) p l t . t i t l e ( "Von M i s e s S t r e s s ( kPa ) " ) p l t . l e g e n d ( l o c = ' b e s t ' ) x l a b e l ( " r (mm) " ) y l a b e l ( "sigma VM ( kPa ) " ) p l t . show ( ) Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 Hyperelastic constitutive law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3 Having some fun with differential calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4 Weak formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Appendix A.1. Boundary conditions A pressure p ext is imposed on Γ e , where x = r e e r :

Π n( x) = p ext , (A.1)

# gr aphe p e r m e t t a n t de s i t u e r ro , ave c Phi ( ro )=0 p l t . f i g u r e [START_REF] Allaire | A review of adjoint methods for sensitivity analysis, uncertainty quantification and optimization in numerical codes[END_REF] x d e f = np . l i n s p a c e ( -0 . 4 , 0 . 4 , 5 0 0 ) y d e f = [ Phi ( i ) f o r i in x d e f ] p l t . p l o t ( x d e f , y d e f , marker="+" ) #p l t . x l i m ( -0 . 1 , 0 . 1 ) #p l t . y l i m ( -10000 ,10000) p l t . t i t l e ( " Phi ( x ) | Phi ( r o ) = 0 " ) p l t . l e g e n d ( ) x l a b e l ( "x" ) y l a b e l ( " Phi " ) p l t . g r i d ( True ) p l t . show ( ) # p l t . p l o t ( x , y , ' r o ' ) p l t . x l i m ( 1 . 9 , 2 . 3 ) p l t . t i t l e ( " D i s p l a c e m e n t u r ( r ) " ) p l t . l e g e n d ( ) x l a b e l ( " r (mm) " ) y l a b e l ( "u (mm) " ) p l t . g r i d ( True ) p l t . show ( ) # ------------------------